--- /dev/null
+//---------------------------- step-10.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- step-10.cc ---------------------------
+
+
+// a hp-ified version of step-10
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-10/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/convergence_table.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <grid/tria.h>
+#include <grid/grid_out.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <fe/fe_q.h>
+#include <fe/hp_fe_values.h>
+
+#include <fe/mapping_q.h>
+
+#include <iostream>
+#include <fstream>
+#include <cmath>
+
+const long double pi = 3.141592653589793238462643;
+
+
+
+template <int dim>
+void gnuplot_output()
+{
+ deallog << "Output of grids into gnuplot files:" << std::endl
+ << "===================================" << std::endl;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_ball (triangulation);
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+
+ for (unsigned int refinement=0; refinement<2;
+ ++refinement, triangulation.refine_global(1))
+ {
+ deallog << "Refinement level: " << refinement << std::endl;
+
+ std::string filename_base = "ball";
+ filename_base += '0'+refinement;
+
+ for (unsigned int degree=1; degree<4; ++degree)
+ {
+ deallog << "Degree = " << degree << std::endl;
+
+ const MappingQ<dim> mapping (degree);
+
+
+ GridOut grid_out;
+ GridOutFlags::Gnuplot gnuplot_flags(false, 30);
+ grid_out.set_flags(gnuplot_flags);
+
+ std::string filename = filename_base+"_mapping_q";
+ filename += ('0'+degree);
+ filename += ".dat";
+ std::ofstream gnuplot_file (filename.c_str());
+
+ grid_out.write_gnuplot (triangulation, gnuplot_file, &mapping);
+ }
+ deallog << std::endl;
+ }
+}
+
+template <int dim>
+void compute_pi_by_area ()
+{
+ deallog << "Computation of Pi by the area:" << std::endl
+ << "==============================" << std::endl;
+
+ const hp::QCollection<dim> quadrature(QGauss<dim>(4));
+
+ for (unsigned int degree=1; degree<5; ++degree)
+ {
+ deallog << "Degree = " << degree << std::endl;
+
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_ball (triangulation);
+
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+
+ MappingQ<dim> m(degree);
+ const hp::MappingCollection<dim> mapping (m);
+
+ const hp::FECollection<dim> dummy_fe (FE_Q<dim>(1));
+
+ hp::DoFHandler<dim> dof_handler (triangulation);
+
+ hp::FEValues<dim> x_fe_values (mapping, dummy_fe, quadrature,
+ update_JxW_values);
+
+ ConvergenceTable table;
+
+ for (unsigned int refinement=0; refinement<6;
+ ++refinement, triangulation.refine_global (1))
+ {
+ table.add_value("cells", triangulation.n_active_cells());
+
+ dof_handler.distribute_dofs (dummy_fe);
+
+ long double area = 0;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+ for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i)
+ area += fe_values.JxW (i);
+ };
+
+ table.add_value("eval.pi", static_cast<double> (area));
+ table.add_value("error", static_cast<double> (std::fabs(area-pi)));
+ };
+
+ table.omit_column_from_convergence_rate_evaluation("cells");
+ table.omit_column_from_convergence_rate_evaluation("eval.pi");
+ table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
+
+ table.set_precision("eval.pi", 16);
+ table.set_scientific("error", true);
+
+ table.write_text(deallog.get_file_stream());
+
+ deallog << std::endl;
+ };
+}
+
+
+template <int dim>
+void compute_pi_by_perimeter ()
+{
+ deallog << "Computation of Pi by the perimeter:" << std::endl
+ << "===================================" << std::endl;
+
+ const hp::QCollection<dim-1> quadrature(4);
+
+ for (unsigned int degree=1; degree<5; ++degree)
+ {
+ deallog << "Degree = " << degree << std::endl;
+ Triangulation<dim> triangulation;
+ GridGenerator::hyper_ball (triangulation);
+
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+
+ const MappingQ<dim> m (degree);
+ const hp::MappingCollection<dim> mapping(m);
+ const hp::FECollection<dim> fe (FE_Q<dim>(1));
+
+ hp::DoFHandler<dim> dof_handler (triangulation);
+
+ hp::FEFaceValues<dim> x_fe_face_values (mapping, fe, quadrature,
+ update_JxW_values);
+ ConvergenceTable table;
+
+ for (unsigned int refinement=0; refinement<6;
+ ++refinement, triangulation.refine_global (1))
+ {
+ table.add_value("cells", triangulation.n_active_cells());
+
+ dof_handler.distribute_dofs (fe);
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ long double perimeter = 0;
+ for (; cell!=endc; ++cell)
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ if (cell->face(face_no)->at_boundary())
+ {
+ x_fe_face_values.reinit (cell, face_no);
+ const FEFaceValues<dim> &fe_face_values
+ = x_fe_face_values.get_present_fe_values ();
+
+ for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i)
+ perimeter += fe_face_values.JxW (i);
+ };
+ table.add_value("eval.pi", static_cast<double> (perimeter/2.));
+ table.add_value("error", static_cast<double> (std::fabs(perimeter/2.-pi)));
+ };
+
+ table.omit_column_from_convergence_rate_evaluation("cells");
+ table.omit_column_from_convergence_rate_evaluation("eval.pi");
+ table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2);
+
+ table.set_precision("eval.pi", 16);
+ table.set_scientific("error", true);
+
+ table.write_text(deallog.get_file_stream());
+
+ deallog << std::endl;
+ };
+}
+
+
+int main ()
+{
+ logfile.precision(2);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ gnuplot_output<2>();
+
+ compute_pi_by_area<2> ();
+ compute_pi_by_perimeter<2> ();
+
+ return 0;
+}
--- /dev/null
+//---------------------------- step-6.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- step-6.cc ---------------------------
+
+
+// a hp-ified version of step-6
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-6/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <dofs/hp_dof_handler.h>
+#include <grid/grid_generator.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+
+#include <fstream>
+#include <iostream>
+
+#include <fe/fe_q.h>
+#include <grid/grid_out.h>
+
+
+#include <dofs/dof_constraints.h>
+
+#include <grid/grid_refinement.h>
+
+#include <numerics/error_estimator.h>
+
+
+
+
+template <int dim>
+class LaplaceProblem
+{
+ public:
+ LaplaceProblem ();
+ ~LaplaceProblem ();
+
+ void run ();
+
+ private:
+ void setup_system ();
+ void assemble_system ();
+ void solve ();
+ void refine_grid ();
+ void output_results (const unsigned int cycle) const;
+
+ Triangulation<dim> triangulation;
+
+ hp::DoFHandler<dim> dof_handler;
+ hp::FECollection<dim> fe;
+
+ ConstraintMatrix hanging_node_constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+
+
+
+template <int dim>
+class Coefficient : public Function<dim>
+{
+ public:
+ Coefficient () : Function<dim>() {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+double Coefficient<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+{
+ if (p.square() < 0.5*0.5)
+ return 20;
+ else
+ return 1;
+}
+
+
+
+template <int dim>
+void Coefficient<dim>::value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component) const
+{
+ const unsigned int n_points = points.size();
+
+ Assert (values.size() == n_points,
+ ExcDimensionMismatch (values.size(), n_points));
+
+ Assert (component == 0,
+ ExcIndexRange (component, 0, 1));
+
+ for (unsigned int i=0; i<n_points; ++i)
+ {
+ if (points[i].square() < 0.5*0.5)
+ values[i] = 20;
+ else
+ values[i] = 1;
+ }
+}
+
+
+
+
+template <int dim>
+LaplaceProblem<dim>::LaplaceProblem () :
+ dof_handler (triangulation),
+ fe (FE_Q<dim>(2))
+{}
+
+
+
+template <int dim>
+LaplaceProblem<dim>::~LaplaceProblem ()
+{
+ dof_handler.clear ();
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+
+
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+
+ hanging_node_constraints.close ();
+
+ hanging_node_constraints.condense (sparsity_pattern);
+
+ sparsity_pattern.compress();
+
+ system_matrix.reinit (sparsity_pattern);
+}
+
+
+template <int dim>
+void LaplaceProblem<dim>::assemble_system ()
+{
+ const hp::QCollection<dim> quadrature_formula(3);
+
+ hp::FEValues<dim> x_fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_q_points | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe[0].dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points;
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ const Coefficient<dim> coefficient;
+ std::vector<double> coefficient_values (n_q_points);
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+
+ coefficient.value_list (fe_values.get_quadrature_points(),
+ coefficient_values);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += (coefficient_values[q_point] *
+ fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ 1.0 *
+ fe_values.JxW(q_point));
+ }
+
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+
+ hanging_node_constraints.condense (system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+
+ std::map<unsigned int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(),
+ boundary_values);
+ MatrixTools::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+}
+
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ SolverCG<> cg (solver_control);
+
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.2);
+
+ cg.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+
+ hanging_node_constraints.distribute (solution);
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::refine_grid ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ QGauss<dim-1>(3),
+ typename FunctionMap<dim>::type(),
+ solution,
+ estimated_error_per_cell);
+
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.03);
+
+ triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
+{
+ Assert (cycle < 10, ExcNotImplemented());
+
+ std::string filename = "grid-";
+ filename += ('0' + cycle);
+ filename += ".eps";
+
+ std::ofstream output (filename.c_str());
+
+ GridOut grid_out;
+ grid_out.write_eps (triangulation, output);
+}
+
+
+
+
+template <int dim>
+void LaplaceProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<8; ++cycle)
+ {
+ deallog << "Cycle " << cycle << ':' << std::endl;
+
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_ball (triangulation);
+
+ static const HyperBallBoundary<dim> boundary;
+ triangulation.set_boundary (0, boundary);
+
+ triangulation.refine_global (1);
+ }
+ else
+ refine_grid ();
+
+
+ deallog << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl;
+
+ setup_system ();
+
+ deallog << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ assemble_system ();
+ solve ();
+ output_results (cycle);
+ }
+
+ DataOutBase::EpsFlags eps_flags;
+ eps_flags.z_scaling = 4;
+
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
+ data_out.set_flags (eps_flags);
+
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+ data_out.build_patches ();
+
+ std::ofstream output ("final-solution.eps");
+ data_out.write_eps (output);
+}
+
+
+
+int main ()
+{
+ logfile.precision(2);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ deallog.depth_console (0);
+
+ LaplaceProblem<2> laplace_problem_2d;
+ laplace_problem_2d.run ();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
--- /dev/null
+//---------------------------- step-7.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- step-7.cc ---------------------------
+
+
+// a hp-ified version of step-7
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-7/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_constraints.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/fe_q.h>
+#include <numerics/matrices.h>
+#include <numerics/error_estimator.h>
+#include <numerics/data_out.h>
+
+#include <dofs/dof_renumbering.h>
+#include <base/smartpointer.h>
+#include <numerics/vectors.h>
+#include <base/convergence_table.h>
+#include <fe/hp_fe_values.h>
+
+#include <typeinfo>
+#include <fstream>
+#include <iostream>
+
+
+
+template <int dim>
+class SolutionBase
+{
+ protected:
+ static const unsigned int n_source_centers = 3;
+ static const Point<dim> source_centers[n_source_centers];
+ static const double width;
+};
+
+
+template <>
+const Point<1>
+SolutionBase<1>::source_centers[SolutionBase<1>::n_source_centers]
+= { Point<1>(-1.0 / 3.0),
+ Point<1>(0.0),
+ Point<1>(+1.0 / 3.0) };
+
+template <>
+const Point<2>
+SolutionBase<2>::source_centers[SolutionBase<2>::n_source_centers]
+= { Point<2>(-0.5, +0.5),
+ Point<2>(-0.5, -0.5),
+ Point<2>(+0.5, -0.5) };
+
+template <int dim>
+const double SolutionBase<dim>::width = 1./3.;
+
+
+
+template <int dim>
+class Solution : public Function<dim>,
+ protected SolutionBase<dim>
+{
+ public:
+ Solution () : Function<dim>() {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual Tensor<1,dim> gradient (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+
+template <int dim>
+double Solution<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+{
+ double return_value = 0;
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> x_minus_xi = p - this->source_centers[i];
+ return_value += std::exp(-x_minus_xi.square() /
+ (this->width * this->width));
+ }
+
+ return return_value;
+}
+
+
+template <int dim>
+Tensor<1,dim> Solution<dim>::gradient (const Point<dim> &p,
+ const unsigned int) const
+{
+ Tensor<1,dim> return_value;
+
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> x_minus_xi = p - this->source_centers[i];
+
+ return_value += (-2 / (this->width * this->width) *
+ std::exp(-x_minus_xi.square() /
+ (this->width * this->width)) *
+ x_minus_xi);
+ }
+
+ return return_value;
+}
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>,
+ protected SolutionBase<dim>
+{
+ public:
+ RightHandSide () : Function<dim>() {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+};
+
+
+template <int dim>
+double RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int) const
+{
+ double return_value = 0;
+ for (unsigned int i=0; i<this->n_source_centers; ++i)
+ {
+ const Point<dim> x_minus_xi = p - this->source_centers[i];
+
+ return_value += ((2*dim - 4*x_minus_xi.square()/
+ (this->width * this->width)) /
+ (this->width * this->width) *
+ std::exp(-x_minus_xi.square() /
+ (this->width * this->width)));
+ return_value += std::exp(-x_minus_xi.square() /
+ (this->width * this->width));
+ }
+
+ return return_value;
+}
+
+
+
+template <int dim>
+class HelmholtzProblem
+{
+ public:
+ enum RefinementMode {
+ global_refinement, adaptive_refinement
+ };
+
+ HelmholtzProblem (const hp::FECollection<dim> &fe,
+ const RefinementMode refinement_mode);
+
+ ~HelmholtzProblem ();
+
+ void run ();
+
+ private:
+ void setup_system ();
+ void assemble_system ();
+ void solve ();
+ void refine_grid ();
+ void process_solution (const unsigned int cycle);
+
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler;
+
+ SmartPointer<const hp::FECollection<dim> > fe;
+
+ ConstraintMatrix hanging_node_constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ const RefinementMode refinement_mode;
+
+ ConvergenceTable convergence_table;
+};
+
+
+
+
+template <int dim>
+HelmholtzProblem<dim>::HelmholtzProblem (const hp::FECollection<dim> &fe,
+ const RefinementMode refinement_mode) :
+ dof_handler (triangulation),
+ fe (&fe),
+ refinement_mode (refinement_mode)
+{}
+
+
+
+template <int dim>
+HelmholtzProblem<dim>::~HelmholtzProblem ()
+{
+ dof_handler.clear ();
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (*fe);
+ DoFRenumbering::Cuthill_McKee (dof_handler);
+
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+ hanging_node_constraints.condense (sparsity_pattern);
+ sparsity_pattern.compress();
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::assemble_system ()
+{
+ hp::QCollection<dim> quadrature_formula(3);
+ hp::QCollection<dim-1> face_quadrature_formula(3);
+
+ const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points;
+ const unsigned int n_face_q_points = face_quadrature_formula[0].n_quadrature_points;
+
+ const unsigned int dofs_per_cell = (*fe)[0].dofs_per_cell;
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ hp::FEValues<dim> x_fe_values (*fe, quadrature_formula,
+ update_values | update_gradients |
+ update_q_points | update_JxW_values);
+
+ hp::FEFaceValues<dim> x_fe_face_values (*fe, face_quadrature_formula,
+ update_values | update_q_points |
+ update_normal_vectors | update_JxW_values);
+
+ const RightHandSide<dim> right_hand_side;
+ std::vector<double> rhs_values (n_q_points);
+
+ const Solution<dim> exact_solution;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator
+ cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+
+ right_hand_side.value_list (fe_values.get_quadrature_points(),
+ rhs_values);
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += ((fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point)
+ +
+ fe_values.shape_value(i,q_point) *
+ fe_values.shape_value(j,q_point)) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += (fe_values.shape_value(i,q_point) *
+ rhs_values [q_point] *
+ fe_values.JxW(q_point));
+ }
+
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->at_boundary()
+ &&
+ (cell->face(face)->boundary_indicator() == 1))
+ {
+ x_fe_face_values.reinit (cell, face);
+ const FEFaceValues<dim> &fe_face_values
+ = x_fe_face_values.get_present_fe_values ();
+
+ for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+ {
+ const double neumann_value
+ = (exact_solution.gradient (fe_face_values.quadrature_point(q_point)) *
+ fe_face_values.normal_vector(q_point));
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ cell_rhs(i) += (neumann_value *
+ fe_face_values.shape_value(i,q_point) *
+ fe_face_values.JxW(q_point));
+ }
+ }
+
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+
+ hanging_node_constraints.condense (system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+
+ std::map<unsigned int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ Solution<dim>(),
+ boundary_values);
+ MatrixTools::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ SolverCG<> cg (solver_control);
+
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.2);
+
+ cg.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+
+ hanging_node_constraints.distribute (solution);
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::refine_grid ()
+{
+ switch (refinement_mode)
+ {
+ case global_refinement:
+ {
+ triangulation.refine_global (1);
+ break;
+ }
+
+ case adaptive_refinement:
+ {
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ typename FunctionMap<dim>::type neumann_boundary;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ hp::QCollection<dim-1>(3),
+ neumann_boundary,
+ solution,
+ estimated_error_per_cell);
+
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.03);
+
+ triangulation.execute_coarsening_and_refinement ();
+
+ break;
+ }
+
+ default:
+ {
+ Assert (false, ExcNotImplemented());
+ }
+ }
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::process_solution (const unsigned int cycle)
+{
+ Vector<float> difference_per_cell (triangulation.n_active_cells());
+ VectorTools::integrate_difference (dof_handler,
+ solution,
+ Solution<dim>(),
+ difference_per_cell,
+ hp::QCollection<dim>(3),
+ VectorTools::L2_norm);
+ const double L2_error = difference_per_cell.l2_norm();
+
+ VectorTools::integrate_difference (dof_handler,
+ solution,
+ Solution<dim>(),
+ difference_per_cell,
+ hp::QCollection<dim>(3),
+ VectorTools::H1_seminorm);
+ const double H1_error = difference_per_cell.l2_norm();
+
+ const QTrapez<1> q_trapez;
+ const QIterated<dim> q_iterated (q_trapez, 5);
+ VectorTools::integrate_difference (dof_handler,
+ solution,
+ Solution<dim>(),
+ difference_per_cell,
+ q_iterated,
+ VectorTools::Linfty_norm);
+ const double Linfty_error = difference_per_cell.linfty_norm();
+
+ const unsigned int n_active_cells=triangulation.n_active_cells();
+ const unsigned int n_dofs=dof_handler.n_dofs();
+
+ deallog << "Cycle " << cycle << ':'
+ << std::endl
+ << " Number of active cells: "
+ << n_active_cells
+ << std::endl
+ << " Number of degrees of freedom: "
+ << n_dofs
+ << std::endl;
+
+ convergence_table.add_value("cycle", cycle);
+ convergence_table.add_value("cells", n_active_cells);
+ convergence_table.add_value("dofs", n_dofs);
+ convergence_table.add_value("L2", L2_error);
+ convergence_table.add_value("H1", H1_error);
+ convergence_table.add_value("Linfty", Linfty_error);
+}
+
+
+
+template <int dim>
+void HelmholtzProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<7; ++cycle)
+ {
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (1);
+
+ typename Triangulation<dim>::cell_iterator
+ cell = triangulation.begin (),
+ endc = triangulation.end();
+ for (; cell!=endc; ++cell)
+ for (unsigned int face=0;
+ face<GeometryInfo<dim>::faces_per_cell;
+ ++face)
+ if ((cell->face(face)->center()(0) == -1)
+ ||
+ (cell->face(face)->center()(1) == -1))
+ cell->face(face)->set_boundary_indicator (1);
+ }
+ else
+ refine_grid ();
+
+
+ setup_system ();
+
+ assemble_system ();
+ solve ();
+
+ process_solution (cycle);
+ }
+
+
+ std::string gmv_filename;
+ switch (refinement_mode)
+ {
+ case global_refinement:
+ gmv_filename = "solution-global";
+ break;
+ case adaptive_refinement:
+ gmv_filename = "solution-adaptive";
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ switch ((*fe)[0].degree)
+ {
+ case 1:
+ gmv_filename += "-q1";
+ break;
+ case 2:
+ gmv_filename += "-q2";
+ break;
+
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ gmv_filename += ".gmv";
+ std::ofstream output (gmv_filename.c_str());
+
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+
+ data_out.build_patches ((*fe)[0].degree);
+ data_out.write_gmv (output);
+
+
+
+ convergence_table.set_precision("L2", 3);
+ convergence_table.set_precision("H1", 3);
+ convergence_table.set_precision("Linfty", 3);
+
+ convergence_table.set_scientific("L2", true);
+ convergence_table.set_scientific("H1", true);
+ convergence_table.set_scientific("Linfty", true);
+
+ convergence_table.set_tex_caption("cells", "\\# cells");
+ convergence_table.set_tex_caption("dofs", "\\# dofs");
+ convergence_table.set_tex_caption("L2", "$L^2$-error");
+ convergence_table.set_tex_caption("H1", "$H^1$-error");
+ convergence_table.set_tex_caption("Linfty", "$L^\\infty$-error");
+
+ convergence_table.set_tex_format("cells", "r");
+ convergence_table.set_tex_format("dofs", "r");
+
+ deallog << std::endl;
+ convergence_table.write_text(deallog);
+
+ std::string error_filename = "error";
+ switch (refinement_mode)
+ {
+ case global_refinement:
+ error_filename += "-global";
+ break;
+ case adaptive_refinement:
+ error_filename += "-adaptive";
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ switch ((*fe)[0].degree)
+ {
+ case 1:
+ error_filename += "-q1";
+ break;
+ case 2:
+ error_filename += "-q2";
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ error_filename += ".tex";
+ std::ofstream error_table_file(error_filename.c_str());
+
+ convergence_table.write_tex(error_table_file);
+
+
+
+ if (refinement_mode==global_refinement)
+ {
+ convergence_table.add_column_to_supercolumn("cycle", "n cells");
+ convergence_table.add_column_to_supercolumn("cells", "n cells");
+
+ std::vector<std::string> new_order;
+ new_order.push_back("n cells");
+ new_order.push_back("H1");
+ new_order.push_back("L2");
+ convergence_table.set_column_order (new_order);
+
+ convergence_table
+ .evaluate_convergence_rates("L2", ConvergenceTable::reduction_rate);
+ convergence_table
+ .evaluate_convergence_rates("L2", ConvergenceTable::reduction_rate_log2);
+ convergence_table
+ .evaluate_convergence_rates("H1", ConvergenceTable::reduction_rate_log2);
+
+ deallog << std::endl;
+ convergence_table.write_text(deallog);
+
+ std::string conv_filename = "convergence";
+ switch (refinement_mode)
+ {
+ case global_refinement:
+ conv_filename += "-global";
+ break;
+ case adaptive_refinement:
+ conv_filename += "-adaptive";
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ switch ((*fe)[0].degree)
+ {
+ case 1:
+ conv_filename += "-q1";
+ break;
+ case 2:
+ conv_filename += "-q2";
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+ conv_filename += ".tex";
+
+ std::ofstream table_file(conv_filename.c_str());
+ convergence_table.write_tex(table_file);
+ }
+}
+
+
+int main ()
+{
+ logfile.precision(2);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ const unsigned int dim = 2;
+
+ try
+ {
+ deallog.depth_console (0);
+
+
+ {
+ deallog << "Solving with Q1 elements, adaptive refinement" << std::endl
+ << "=============================================" << std::endl
+ << std::endl;
+
+ hp::FECollection<dim> fe(FE_Q<dim>(1));
+ HelmholtzProblem<dim>
+ helmholtz_problem_2d (fe, HelmholtzProblem<dim>::adaptive_refinement);
+
+ helmholtz_problem_2d.run ();
+
+ deallog << std::endl;
+ }
+
+ {
+ deallog << "Solving with Q1 elements, global refinement" << std::endl
+ << "===========================================" << std::endl
+ << std::endl;
+
+ hp::FECollection<dim> fe(FE_Q<dim>(1));
+ HelmholtzProblem<dim>
+ helmholtz_problem_2d (fe, HelmholtzProblem<dim>::global_refinement);
+
+ helmholtz_problem_2d.run ();
+
+ deallog << std::endl;
+ }
+
+ {
+ deallog << "Solving with Q2 elements, global refinement" << std::endl
+ << "===========================================" << std::endl
+ << std::endl;
+
+ hp::FECollection<dim> fe(FE_Q<dim>(2));
+ HelmholtzProblem<dim>
+ helmholtz_problem_2d (fe, HelmholtzProblem<dim>::global_refinement);
+
+ helmholtz_problem_2d.run ();
+
+ deallog << std::endl;
+ }
+
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
+
+
+template const double SolutionBase<2>::width;
--- /dev/null
+//---------------------------- step-8.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- step-8.cc ---------------------------
+
+
+// a hp-ified version of step-8
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-8/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_cg.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <fe/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <dofs/dof_constraints.h>
+#include <numerics/error_estimator.h>
+
+#include <fe/fe_system.h>
+#include <fe/fe_q.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+template <int dim>
+class ElasticProblem
+{
+ public:
+ ElasticProblem ();
+ ~ElasticProblem ();
+ void run ();
+
+ private:
+ void setup_system ();
+ void assemble_system ();
+ void solve ();
+ void refine_grid ();
+ void output_results (const unsigned int cycle) const;
+
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler;
+
+ hp::FECollection<dim> fe;
+
+ ConstraintMatrix hanging_node_constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+};
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+ public:
+ RightHandSide ();
+
+ virtual void vector_value (const Point<dim> &p,
+ Vector<double> &values) const;
+
+ virtual void vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &value_list) const;
+};
+
+
+template <int dim>
+RightHandSide<dim>::RightHandSide ()
+ :
+ Function<dim> (dim)
+{}
+
+
+template <int dim>
+inline
+void RightHandSide<dim>::vector_value (const Point<dim> &p,
+ Vector<double> &values) const
+{
+ Assert (values.size() == dim,
+ ExcDimensionMismatch (values.size(), dim));
+ Assert (dim >= 2, ExcNotImplemented());
+
+ Point<dim> point_1, point_2;
+ point_1(0) = 0.5;
+ point_2(0) = -0.5;
+
+ if (((p-point_1).square() < 0.2*0.2) ||
+ ((p-point_2).square() < 0.2*0.2))
+ values(0) = 1;
+ else
+ values(0) = 0;
+
+ if (p.square() < 0.2*0.2)
+ values(1) = 1;
+ else
+ values(1) = 0;
+}
+
+
+
+template <int dim>
+void RightHandSide<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+ std::vector<Vector<double> > &value_list) const
+{
+ Assert (value_list.size() == points.size(),
+ ExcDimensionMismatch (value_list.size(), points.size()));
+
+ const unsigned int n_points = points.size();
+
+ for (unsigned int p=0; p<n_points; ++p)
+ RightHandSide<dim>::vector_value (points[p],
+ value_list[p]);
+}
+
+
+
+
+
+
+template <int dim>
+ElasticProblem<dim>::ElasticProblem ()
+ :
+ dof_handler (triangulation),
+ fe (FESystem<dim>(FE_Q<dim>(1), dim))
+{}
+
+
+
+
+template <int dim>
+ElasticProblem<dim>::~ElasticProblem ()
+{
+ dof_handler.clear ();
+}
+
+
+
+template <int dim>
+void ElasticProblem<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+ hanging_node_constraints.condense (sparsity_pattern);
+
+ sparsity_pattern.compress();
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void ElasticProblem<dim>::assemble_system ()
+{
+ hp::QCollection<dim> quadrature_formula(2);
+
+ hp::FEValues<dim> x_fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_q_points | update_JxW_values);
+
+ const unsigned int dofs_per_cell = fe[0].dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points;
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ std::vector<double> lambda_values (n_q_points);
+ std::vector<double> mu_values (n_q_points);
+
+ ConstantFunction<dim> lambda(1.), mu(1.);
+
+ RightHandSide<dim> right_hand_side;
+ std::vector<Vector<double> > rhs_values (n_q_points,
+ Vector<double>(dim));
+
+
+ typename hp::DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(),
+ endc = dof_handler.end();
+ for (; cell!=endc; ++cell)
+ {
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+
+ lambda.value_list (fe_values.get_quadrature_points(), lambda_values);
+ mu.value_list (fe_values.get_quadrature_points(), mu_values);
+
+ right_hand_side.vector_value_list (fe_values.get_quadrature_points(),
+ rhs_values);
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const unsigned int
+ component_i = fe[0].system_to_component_index(i).first;
+
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ {
+ const unsigned int
+ component_j = fe[0].system_to_component_index(j).first;
+
+ for (unsigned int q_point=0; q_point<n_q_points;
+ ++q_point)
+ {
+ cell_matrix(i,j)
+ +=
+ (
+ (fe_values.shape_grad(i,q_point)[component_i] *
+ fe_values.shape_grad(j,q_point)[component_j] *
+ lambda_values[q_point])
+ +
+ (fe_values.shape_grad(i,q_point)[component_j] *
+ fe_values.shape_grad(j,q_point)[component_i] *
+ mu_values[q_point])
+ +
+ ((component_i == component_j) ?
+ (fe_values.shape_grad(i,q_point) *
+ fe_values.shape_grad(j,q_point) *
+ mu_values[q_point]) :
+ 0)
+ )
+ *
+ fe_values.JxW(q_point);
+ }
+ }
+ }
+
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ const unsigned int
+ component_i = fe[0].system_to_component_index(i).first;
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ cell_rhs(i) += fe_values.shape_value(i,q_point) *
+ rhs_values[q_point](component_i) *
+ fe_values.JxW(q_point);
+ }
+
+ cell->get_dof_indices (local_dof_indices);
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ }
+ }
+
+ hanging_node_constraints.condense (system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+
+ std::map<unsigned int,double> boundary_values;
+ VectorTools::interpolate_boundary_values (dof_handler,
+ 0,
+ ZeroFunction<dim>(dim),
+ boundary_values);
+ MatrixTools::apply_boundary_values (boundary_values,
+ system_matrix,
+ solution,
+ system_rhs);
+}
+
+
+
+
+template <int dim>
+void ElasticProblem<dim>::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ SolverCG<> cg (solver_control);
+
+ PreconditionSSOR<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.2);
+
+ cg.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+
+ hanging_node_constraints.distribute (solution);
+}
+
+
+
+template <int dim>
+void ElasticProblem<dim>::refine_grid ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ typename FunctionMap<dim>::type neumann_boundary;
+ KellyErrorEstimator<dim>::estimate (dof_handler,
+ hp::QCollection<dim-1>(2),
+ neumann_boundary,
+ solution,
+ estimated_error_per_cell);
+
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.3, 0.03);
+
+ triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+template <int dim>
+void ElasticProblem<dim>::output_results (const unsigned int cycle) const
+{
+ std::string filename = "solution-";
+ filename += ('0' + cycle);
+ Assert (cycle < 10, ExcInternalError());
+
+ filename += ".gmv";
+ std::ofstream output (filename.c_str());
+
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
+ data_out.attach_dof_handler (dof_handler);
+
+
+
+ std::vector<std::string> solution_names;
+ switch (dim)
+ {
+ case 1:
+ solution_names.push_back ("displacement");
+ break;
+ case 2:
+ solution_names.push_back ("x_displacement");
+ solution_names.push_back ("y_displacement");
+ break;
+ case 3:
+ solution_names.push_back ("x_displacement");
+ solution_names.push_back ("y_displacement");
+ solution_names.push_back ("z_displacement");
+ break;
+ default:
+ Assert (false, ExcNotImplemented());
+ }
+
+ data_out.add_data_vector (solution, solution_names);
+ data_out.build_patches ();
+ data_out.write_gmv (output);
+}
+
+
+
+
+template <int dim>
+void ElasticProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<8; ++cycle)
+ {
+ deallog << "Cycle " << cycle << ':' << std::endl;
+
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (2);
+ }
+ else
+ refine_grid ();
+
+ deallog << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl;
+
+ setup_system ();
+
+ deallog << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ assemble_system ();
+ solve ();
+ output_results (cycle);
+ }
+}
+
+
+int main ()
+{
+ logfile.precision(2);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ deallog.depth_console (0);
+
+ ElasticProblem<2> elastic_problem_2d;
+ elastic_problem_2d.run ();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+
+ return 0;
+}
--- /dev/null
+//---------------------------- step-9.cc ---------------------------
+// $Id$
+// Version: $Name$
+//
+// Copyright (C) 2005, 2006 by the deal.II authors
+//
+// This file is subject to QPL and may not be distributed
+// without copyright and license information. Please refer
+// to the file deal.II/doc/license.html for the text and
+// further information on this license.
+//
+//---------------------------- step-9.cc ---------------------------
+
+
+// a hp-ified version of step-9
+
+
+#include <base/logstream.h>
+#include <fstream>
+std::ofstream logfile("step-9/output");
+
+
+#include <base/quadrature_lib.h>
+#include <base/function.h>
+#include <base/logstream.h>
+#include <lac/vector.h>
+#include <lac/full_matrix.h>
+#include <lac/sparse_matrix.h>
+#include <lac/solver_bicgstab.h>
+#include <lac/precondition.h>
+#include <grid/tria.h>
+#include <grid/grid_generator.h>
+#include <grid/grid_refinement.h>
+#include <grid/tria_accessor.h>
+#include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
+#include <dofs/hp_dof_handler.h>
+#include <dofs/dof_accessor.h>
+#include <dofs/dof_tools.h>
+#include <dofs/dof_constraints.h>
+#include <fe/hp_fe_values.h>
+#include <numerics/vectors.h>
+#include <numerics/matrices.h>
+#include <numerics/data_out.h>
+#include <fe/fe_q.h>
+#include <grid/grid_out.h>
+
+#include <base/thread_management.h>
+#include <base/multithread_info.h>
+
+#include <base/tensor_function.h>
+
+#include <numerics/error_estimator.h>
+
+#include <fstream>
+#include <iostream>
+
+
+
+#ifndef M_PI
+# define M_PI 3.14159265358979323846
+#endif
+
+
+template <int dim>
+class AdvectionProblem
+{
+ public:
+ AdvectionProblem ();
+ ~AdvectionProblem ();
+ void run ();
+
+ private:
+ void setup_system ();
+ void assemble_system ();
+ void assemble_system_interval (const typename hp::DoFHandler<dim>::active_cell_iterator &begin,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &end);
+
+ void solve ();
+ void refine_grid ();
+ void output_results (const unsigned int cycle) const;
+
+ Triangulation<dim> triangulation;
+ hp::DoFHandler<dim> dof_handler;
+
+ hp::FECollection<dim> fe;
+
+ ConstraintMatrix hanging_node_constraints;
+
+ SparsityPattern sparsity_pattern;
+ SparseMatrix<double> system_matrix;
+
+ Vector<double> solution;
+ Vector<double> system_rhs;
+
+ Threads::ThreadMutex assembler_lock;
+};
+
+
+
+
+template <int dim>
+class AdvectionField : public TensorFunction<1,dim>
+{
+ public:
+ AdvectionField () : TensorFunction<1,dim> () {};
+
+ virtual Tensor<1,dim> value (const Point<dim> &p) const;
+
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<Tensor<1,dim> > &values) const;
+
+ DeclException2 (ExcDimensionMismatch,
+ unsigned int, unsigned int,
+ << "The vector has size " << arg1 << " but should have "
+ << arg2 << " elements.");
+};
+
+
+
+template <int dim>
+Tensor<1,dim>
+AdvectionField<dim>::value (const Point<dim> &p) const
+{
+ Point<dim> value;
+ value[0] = 2;
+ for (unsigned int i=1; i<dim; ++i)
+ value[i] = 1+0.8*std::sin(8*M_PI*p[0]);
+
+ return value;
+}
+
+
+
+template <int dim>
+void
+AdvectionField<dim>::value_list (const std::vector<Point<dim> > &points,
+ std::vector<Tensor<1,dim> > &values) const
+{
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch (values.size(), points.size()));
+
+ for (unsigned int i=0; i<points.size(); ++i)
+ values[i] = AdvectionField<dim>::value (points[i]);
+}
+
+
+
+
+template <int dim>
+class RightHandSide : public Function<dim>
+{
+ public:
+ RightHandSide () : Function<dim>() {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+
+ private:
+ static const Point<dim> center_point;
+};
+
+
+template <>
+const Point<1> RightHandSide<1>::center_point = Point<1> (-0.75);
+
+template <>
+const Point<2> RightHandSide<2>::center_point = Point<2> (-0.75, -0.75);
+
+template <>
+const Point<3> RightHandSide<3>::center_point = Point<3> (-0.75, -0.75, -0.75);
+
+
+
+template <int dim>
+double
+RightHandSide<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+ const double diameter = 0.1;
+ return ( (p-center_point).square() < diameter*diameter ?
+ .1/std::pow(diameter,dim) :
+ 0);
+}
+
+
+
+template <int dim>
+void
+RightHandSide<dim>::value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component) const
+{
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch (values.size(), points.size()));
+
+ for (unsigned int i=0; i<points.size(); ++i)
+ values[i] = RightHandSide<dim>::value (points[i], component);
+}
+
+
+
+template <int dim>
+class BoundaryValues : public Function<dim>
+{
+ public:
+ BoundaryValues () : Function<dim>() {};
+
+ virtual double value (const Point<dim> &p,
+ const unsigned int component = 0) const;
+
+ virtual void value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component = 0) const;
+};
+
+
+
+template <int dim>
+double
+BoundaryValues<dim>::value (const Point<dim> &p,
+ const unsigned int component) const
+{
+ Assert (component == 0, ExcIndexRange (component, 0, 1));
+
+ const double sine_term = std::sin(16*M_PI*std::sqrt(p.square()));
+ const double weight = std::exp(-5*p.square()) / std::exp(-5.);
+ return sine_term * weight;
+}
+
+
+
+template <int dim>
+void
+BoundaryValues<dim>::value_list (const std::vector<Point<dim> > &points,
+ std::vector<double> &values,
+ const unsigned int component) const
+{
+ Assert (values.size() == points.size(),
+ ExcDimensionMismatch (values.size(), points.size()));
+
+ for (unsigned int i=0; i<points.size(); ++i)
+ values[i] = BoundaryValues<dim>::value (points[i], component);
+}
+
+
+
+
+class GradientEstimation
+{
+ public:
+ template <int dim>
+ static void estimate (const hp::DoFHandler<dim> &dof,
+ const Vector<double> &solution,
+ Vector<float> &error_per_cell);
+
+ DeclException2 (ExcInvalidVectorLength,
+ int, int,
+ << "Vector has length " << arg1 << ", but should have "
+ << arg2);
+ DeclException0 (ExcInsufficientDirections);
+
+ private:
+ typedef std::pair<unsigned int,unsigned int> IndexInterval;
+
+ template <int dim>
+ static void estimate_interval (const hp::DoFHandler<dim> &dof,
+ const Vector<double> &solution,
+ const IndexInterval &index_interval,
+ Vector<float> &error_per_cell);
+};
+
+
+
+
+
+template <int dim>
+AdvectionProblem<dim>::AdvectionProblem () :
+ dof_handler (triangulation),
+ fe(FE_Q<dim>(1))
+{}
+
+
+
+template <int dim>
+AdvectionProblem<dim>::~AdvectionProblem ()
+{
+ dof_handler.clear ();
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::setup_system ()
+{
+ dof_handler.distribute_dofs (fe);
+
+ hanging_node_constraints.clear ();
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ hanging_node_constraints);
+ hanging_node_constraints.close ();
+
+ sparsity_pattern.reinit (dof_handler.n_dofs(),
+ dof_handler.n_dofs(),
+ dof_handler.max_couplings_between_dofs());
+ DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
+
+ hanging_node_constraints.condense (sparsity_pattern);
+
+ sparsity_pattern.compress();
+
+ system_matrix.reinit (sparsity_pattern);
+
+ solution.reinit (dof_handler.n_dofs());
+ system_rhs.reinit (dof_handler.n_dofs());
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::assemble_system ()
+{
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ Threads::ThreadGroup<> threads;
+
+ typedef typename hp::DoFHandler<dim>::active_cell_iterator active_cell_iterator;
+ std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
+ thread_ranges
+ = Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
+ dof_handler.end (),
+ n_threads);
+
+ for (unsigned int thread=0; thread<n_threads; ++thread)
+ threads += Threads::spawn (*this, &AdvectionProblem<dim>::assemble_system_interval)
+ (thread_ranges[thread].first,
+ thread_ranges[thread].second);
+
+ threads.join_all ();
+
+
+ hanging_node_constraints.condense (system_matrix);
+ hanging_node_constraints.condense (system_rhs);
+}
+
+
+
+template <int dim>
+void
+AdvectionProblem<dim>::
+assemble_system_interval (const typename hp::DoFHandler<dim>::active_cell_iterator &begin,
+ const typename hp::DoFHandler<dim>::active_cell_iterator &end)
+{
+ const AdvectionField<dim> advection_field;
+ const RightHandSide<dim> right_hand_side;
+ const BoundaryValues<dim> boundary_values;
+
+ hp::QCollection<dim> quadrature_formula(QGauss<dim>(2));
+ hp::QCollection<dim-1> face_quadrature_formula(QGauss<dim-1>(2));
+
+ hp::FEValues<dim> x_fe_values (fe, quadrature_formula,
+ update_values | update_gradients |
+ update_q_points | update_JxW_values);
+ hp::FEFaceValues<dim> x_fe_face_values (fe, face_quadrature_formula,
+ update_values | update_q_points |
+ update_JxW_values | update_normal_vectors);
+
+ const unsigned int dofs_per_cell = fe[0].dofs_per_cell;
+ const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points;
+ const unsigned int n_face_q_points = face_quadrature_formula[0].n_quadrature_points;
+
+ FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+ Vector<double> cell_rhs (dofs_per_cell);
+
+ std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
+ std::vector<double> rhs_values (n_q_points);
+ std::vector<Tensor<1,dim> > advection_directions (n_q_points);
+ std::vector<double> face_boundary_values (n_face_q_points);
+ std::vector<Tensor<1,dim> > face_advection_directions (n_face_q_points);
+
+ typename hp::DoFHandler<dim>::active_cell_iterator cell;
+ for (cell=begin; cell!=end; ++cell)
+ {
+ cell_matrix = 0;
+ cell_rhs = 0;
+
+ x_fe_values.reinit (cell);
+ const FEValues<dim> &fe_values = x_fe_values.get_present_fe_values();
+
+ advection_field.value_list (fe_values.get_quadrature_points(),
+ advection_directions);
+ right_hand_side.value_list (fe_values.get_quadrature_points(),
+ rhs_values);
+
+ const double delta = 0.1 * cell->diameter ();
+
+ for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) += ((advection_directions[q_point] *
+ fe_values.shape_grad(j,q_point) *
+ (fe_values.shape_value(i,q_point) +
+ delta *
+ (advection_directions[q_point] *
+ fe_values.shape_grad(i,q_point)))) *
+ fe_values.JxW(q_point));
+
+ cell_rhs(i) += ((fe_values.shape_value(i,q_point) +
+ delta *
+ (advection_directions[q_point] *
+ fe_values.shape_grad(i,q_point)) ) *
+ rhs_values[i] *
+ fe_values.JxW (q_point));
+ };
+
+ for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+ if (cell->face(face)->at_boundary())
+ {
+ x_fe_face_values.reinit (cell, face);
+ const FEFaceValues<dim> &fe_face_values
+ = x_fe_face_values.get_present_fe_values ();
+
+ boundary_values.value_list (fe_face_values.get_quadrature_points(),
+ face_boundary_values);
+ advection_field.value_list (fe_face_values.get_quadrature_points(),
+ face_advection_directions);
+
+ for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+ if (fe_face_values.normal_vector(q_point) *
+ face_advection_directions[q_point]
+ < 0)
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ cell_matrix(i,j) -= (face_advection_directions[q_point] *
+ fe_face_values.normal_vector(q_point) *
+ fe_face_values.shape_value(i,q_point) *
+ fe_face_values.shape_value(j,q_point) *
+ fe_face_values.JxW(q_point));
+
+ cell_rhs(i) -= (face_advection_directions[q_point] *
+ fe_face_values.normal_vector(q_point) *
+ face_boundary_values[q_point] *
+ fe_face_values.shape_value(i,q_point) *
+ fe_face_values.JxW(q_point));
+ };
+ };
+
+
+ cell->get_dof_indices (local_dof_indices);
+
+ assembler_lock.acquire ();
+ for (unsigned int i=0; i<dofs_per_cell; ++i)
+ {
+ for (unsigned int j=0; j<dofs_per_cell; ++j)
+ system_matrix.add (local_dof_indices[i],
+ local_dof_indices[j],
+ cell_matrix(i,j));
+
+ system_rhs(local_dof_indices[i]) += cell_rhs(i);
+ };
+ assembler_lock.release ();
+ };
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::solve ()
+{
+ SolverControl solver_control (1000, 1e-12);
+ SolverBicgstab<> bicgstab (solver_control);
+
+ PreconditionJacobi<> preconditioner;
+ preconditioner.initialize(system_matrix, 1.0);
+
+ bicgstab.solve (system_matrix, solution, system_rhs,
+ preconditioner);
+
+ hanging_node_constraints.distribute (solution);
+}
+
+
+template <int dim>
+void AdvectionProblem<dim>::refine_grid ()
+{
+ Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
+
+ GradientEstimation::estimate (dof_handler,
+ solution,
+ estimated_error_per_cell);
+
+ GridRefinement::refine_and_coarsen_fixed_number (triangulation,
+ estimated_error_per_cell,
+ 0.5, 0.03);
+
+ triangulation.execute_coarsening_and_refinement ();
+}
+
+
+
+template <int dim>
+void AdvectionProblem<dim>::output_results (const unsigned int cycle) const
+{
+ std::string filename = "grid-";
+ filename += ('0' + cycle);
+ Assert (cycle < 10, ExcInternalError());
+
+ filename += ".eps";
+ std::ofstream output (filename.c_str());
+
+ GridOut grid_out;
+ grid_out.write_eps (triangulation, output);
+}
+
+
+template <int dim>
+void AdvectionProblem<dim>::run ()
+{
+ for (unsigned int cycle=0; cycle<6; ++cycle)
+ {
+ deallog << "Cycle " << cycle << ':' << std::endl;
+
+ if (cycle == 0)
+ {
+ GridGenerator::hyper_cube (triangulation, -1, 1);
+ triangulation.refine_global (4);
+ }
+ else
+ {
+ refine_grid ();
+ };
+
+
+ deallog << " Number of active cells: "
+ << triangulation.n_active_cells()
+ << std::endl;
+
+ setup_system ();
+
+ deallog << " Number of degrees of freedom: "
+ << dof_handler.n_dofs()
+ << std::endl;
+
+ assemble_system ();
+ solve ();
+ output_results (cycle);
+ };
+
+ DataOut<dim,hp::DoFHandler<dim> > data_out;
+ data_out.attach_dof_handler (dof_handler);
+ data_out.add_data_vector (solution, "solution");
+ data_out.build_patches ();
+
+ std::ofstream output ("final-solution.gmv");
+ data_out.write_gmv (output);
+}
+
+
+
+
+template <int dim>
+void
+GradientEstimation::estimate (const hp::DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution,
+ Vector<float> &error_per_cell)
+{
+ Assert (error_per_cell.size() == dof_handler.get_tria().n_active_cells(),
+ ExcInvalidVectorLength (error_per_cell.size(),
+ dof_handler.get_tria().n_active_cells()));
+
+ const unsigned int n_threads = multithread_info.n_default_threads;
+ std::vector<IndexInterval> index_intervals
+ = Threads::split_interval (0, dof_handler.get_tria().n_active_cells(),
+ n_threads);
+
+ Threads::ThreadGroup<> threads;
+ void (*estimate_interval_ptr) (const hp::DoFHandler<dim> &,
+ const Vector<double> &,
+ const IndexInterval &,
+ Vector<float> &)
+ = &GradientEstimation::template estimate_interval<dim>;
+ for (unsigned int i=0; i<n_threads; ++i)
+ threads += Threads::spawn (estimate_interval_ptr)(dof_handler, solution,
+ index_intervals[i],
+ error_per_cell);
+ threads.join_all ();
+}
+
+
+template <int dim>
+void
+GradientEstimation::estimate_interval (const hp::DoFHandler<dim> &dof_handler,
+ const Vector<double> &solution,
+ const IndexInterval &index_interval,
+ Vector<float> &error_per_cell)
+{
+ QMidpoint<dim> xmidpoint;
+ hp::QCollection<dim> midpoint_rule (xmidpoint);
+ hp::FEValues<dim> x_fe_midpoint_value (dof_handler.get_fe(),
+ midpoint_rule,
+ update_values | update_q_points);
+
+ Tensor<2,dim> Y;
+
+ typename hp::DoFHandler<dim>::active_cell_iterator cell, endc;
+
+ cell = dof_handler.begin_active();
+ advance (cell, static_cast<signed int>(index_interval.first));
+
+ endc = dof_handler.begin_active();
+ advance (endc, static_cast<signed int>(index_interval.second));
+
+ Vector<float>::iterator
+ error_on_this_cell = error_per_cell.begin() + index_interval.first;
+
+
+ std::vector<typename hp::DoFHandler<dim>::active_cell_iterator> active_neighbors;
+ active_neighbors.reserve (GeometryInfo<dim>::faces_per_cell *
+ GeometryInfo<dim>::subfaces_per_face);
+
+ for (; cell!=endc; ++cell, ++error_on_this_cell)
+ {
+ x_fe_midpoint_value.reinit (cell);
+ const FEValues<dim> &fe_midpoint_value = x_fe_midpoint_value.get_present_fe_values ();
+
+ Y.clear ();
+
+ Tensor<1,dim> projected_gradient;
+
+
+ active_neighbors.clear ();
+ for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no)
+ if (! cell->at_boundary(face_no))
+ {
+ const typename hp::DoFHandler<dim>::face_iterator
+ face = cell->face(face_no);
+ const typename hp::DoFHandler<dim>::cell_iterator
+ neighbor = cell->neighbor(face_no);
+
+ if (neighbor->active())
+ active_neighbors.push_back (neighbor);
+ else
+ {
+ if (dim == 1)
+ {
+ typename hp::DoFHandler<dim>::cell_iterator
+ neighbor_child = neighbor;
+ while (neighbor_child->has_children())
+ neighbor_child = neighbor_child->child (face_no==0 ? 1 : 0);
+
+ Assert (neighbor_child->neighbor(face_no==0 ? 1 : 0)==cell,
+ ExcInternalError());
+
+ active_neighbors.push_back (neighbor_child);
+ }
+ else
+ for (unsigned int subface_no=0; subface_no<face->n_children(); ++subface_no)
+ active_neighbors.push_back (
+ cell->neighbor_child_on_subface(face_no, subface_no));
+ };
+ };
+
+ const Point<dim> this_center = fe_midpoint_value.quadrature_point(0);
+
+ std::vector<double> this_midpoint_value(1);
+ fe_midpoint_value.get_function_values (solution, this_midpoint_value);
+
+
+ std::vector<double> neighbor_midpoint_value(1);
+ typename std::vector<typename hp::DoFHandler<dim>::active_cell_iterator>::const_iterator
+ neighbor_ptr = active_neighbors.begin();
+ for (; neighbor_ptr!=active_neighbors.end(); ++neighbor_ptr)
+ {
+ const typename hp::DoFHandler<dim>::active_cell_iterator
+ neighbor = *neighbor_ptr;
+
+ x_fe_midpoint_value.reinit (neighbor);
+ const FEValues<dim> &fe_midpoint_value = x_fe_midpoint_value.get_present_fe_values ();
+ const Point<dim> neighbor_center = fe_midpoint_value.quadrature_point(0);
+
+ fe_midpoint_value.get_function_values (solution,
+ neighbor_midpoint_value);
+
+ Point<dim> y = neighbor_center - this_center;
+ const double distance = std::sqrt(y.square());
+ y /= distance;
+
+ for (unsigned int i=0; i<dim; ++i)
+ for (unsigned int j=0; j<dim; ++j)
+ Y[i][j] += y[i] * y[j];
+
+ projected_gradient += (neighbor_midpoint_value[0] -
+ this_midpoint_value[0]) /
+ distance *
+ y;
+ };
+
+ AssertThrow (determinant(Y) != 0,
+ ExcInsufficientDirections());
+
+ const Tensor<2,dim> Y_inverse = invert(Y);
+
+ Point<dim> gradient;
+ contract (gradient, Y_inverse, projected_gradient);
+
+ *error_on_this_cell = (std::pow(cell->diameter(),
+ 1+1.0*dim/2) *
+ std::sqrt(gradient.square()));
+ };
+}
+
+
+
+int main ()
+{
+ logfile.precision(2);
+
+ deallog.attach(logfile);
+ deallog.depth_console(0);
+ deallog.threshold_double(1.e-10);
+
+ try
+ {
+ deallog.depth_console (0);
+
+ AdvectionProblem<2> advection_problem_2d;
+ advection_problem_2d.run ();
+ }
+ catch (std::exception &exc)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Exception on processing: " << std::endl
+ << exc.what() << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ }
+ catch (...)
+ {
+ std::cerr << std::endl << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ std::cerr << "Unknown exception!" << std::endl
+ << "Aborting!" << std::endl
+ << "----------------------------------------------------"
+ << std::endl;
+ return 1;
+ };
+
+ return 0;
+}
projects / dealii.git / commitdiff