//---------------------------- hp_hanging_node_constraints_02.cc ---------------------------
-// modified step-3 to be a L_2 projection. This function
-// L_2-projects a bilinear, biquadratic or bicubic function
-// onto a square shaped domain which is filled with
-// hp elements. Afterwards the L_2 error is computed, which
-// should be in the range of the FP approximation order.
+// compute hanging node constraints (before and after the processing
+// we do in ConstraintMatrix::close()) for some grids with and without
+// random distribution of FEs
#include <base/logstream.h>
#include <grid/tria.h>
template <int dim>
-class Linear : public Function<dim>
+void run (bool random_p,
+ unsigned int *indx)
{
- public:
- Linear () : Function<dim>() {};
-
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-};
-
-
-
-template <int dim>
-double Linear<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- double res = p(0);
- for (unsigned int d = 1; d < dim; ++d)
- res *= p(d);
-
- return res;
-}
-
-
-
-template <int dim>
-class Quadratic : public Function<dim>
-{
- public:
- Quadratic () : Function<dim>() {};
-
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-};
-
-
-
-template <int dim>
-double Quadratic<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- double res = p(0) * p(0);
- for (unsigned int d = 1; d < dim; ++d)
- res *= p(d) * p(d);
-
- return res;
-}
-
-
-template <int dim>
-class Cubic : public Function<dim>
-{
- public:
- Cubic () : Function<dim>() {};
-
- virtual double value (const Point<dim> &p,
- const unsigned int component = 0) const;
-};
-
+ Triangulation<dim> triangulation;
+ hp::FECollection<dim> fe;
+ hp::DoFHandler<dim> dof_handler(triangulation);
+ ConstraintMatrix hanging_node_constraints;
+ FE_Q<dim> fe_1 (indx[0]),
+ fe_2 (indx[1]),
+ fe_3 (indx[2]),
+ fe_4 (indx[3]);
-template <int dim>
-double Cubic<dim>::value (const Point<dim> &p,
- const unsigned int) const
-{
- double res = p(0) * p(0) * p(0);
- for (unsigned int d = 1; d < dim; ++d)
- res *= p(d) * p(d) * p(d);
+ fe.push_back (fe_1);
+ fe.push_back (fe_2);
+ fe.push_back (fe_3);
+ fe.push_back (fe_4);
- return res;
-}
-
-
-template <int dim>
-class LaplaceProblem
-{
- public:
- LaplaceProblem ();
-
- void run (const Function<dim> &f_test,
- bool random,
- unsigned int *indx);
-
- private:
- void make_grid_and_dofs (const bool random_p);
- void assemble_system (const Function<dim> &f_test);
- void solve ();
- void eval_error (const Function<dim> &f_test);
- void output_results () const;
-
- Triangulation<dim> triangulation;
- hp::FECollection<dim> fe;
- hp::DoFHandler<dim> dof_handler;
-
- SparsityPattern sparsity_pattern;
- SparseMatrix<double> system_matrix;
-
- // Although we do not have h-refinement,
- // hanging nodes will inevitably appear
- // due to different polynomial degrees.
- ConstraintMatrix hanging_node_constraints;
-
- Vector<double> solution;
- Vector<double> system_rhs;
-};
-
-
-template <int dim>
-LaplaceProblem<dim>::LaplaceProblem () :
- dof_handler (triangulation)
-{}
-
-
-template <int dim>
-void LaplaceProblem<dim>::make_grid_and_dofs (const bool random_p)
-{
GridGenerator::hyper_cube (triangulation, -1, 1);
triangulation.refine_global (5-dim);
deallog << "Number of active cells: "
<< dof_handler.n_dofs()
<< std::endl;
- solution.reinit (dof_handler.n_dofs());
- system_rhs.reinit (dof_handler.n_dofs());
-
- // Create sparsity pattern.
- 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);
-
// Create constraints which stem from
// the different polynomial degrees on
// the different elements.
hanging_node_constraints.print (deallog.get_file_stream ());
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 Function<dim> &f_test)
-{
- hp::QCollection<dim> quadrature_formula(QGauss<dim>(6));
- hp::FEValues<dim> x_fe_values (fe, quadrature_formula,
- update_values | update_gradients |
- update_JxW_values | update_q_points);
-
- const unsigned int max_dofs_per_cell = fe.max_dofs_per_cell ();
- const unsigned int n_q_points = quadrature_formula[0].n_quadrature_points;
-
- FullMatrix<double> cell_matrix (max_dofs_per_cell, max_dofs_per_cell);
- Vector<double> cell_rhs (max_dofs_per_cell);
-
- std::vector<unsigned int> local_dof_indices (max_dofs_per_cell);
-
- 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();
-
- cell_matrix = 0;
- cell_rhs = 0;
-
- const unsigned int dofs_per_cell = cell->get_fe ().dofs_per_cell;
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int j=0; j<dofs_per_cell; ++j)
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- cell_matrix(i,j) += (fe_values.shape_value (i, q_point) *
- fe_values.shape_value (j, q_point) *
- fe_values.JxW (q_point));
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
- {
- Point<dim> p_q = fe_values.quadrature_point (q_point);
- double f = f_test.value (p_q);
- cell_rhs(i) += (fe_values.shape_value (i, q_point) *
- f *
- fe_values.JxW (q_point));
- }
-
- local_dof_indices.resize (dofs_per_cell);
- 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));
-
- for (unsigned int i=0; i<dofs_per_cell; ++i)
- system_rhs(local_dof_indices[i]) += cell_rhs(i);
- }
-
- // Include hanging nodes.
- hanging_node_constraints.condense (system_matrix);
- hanging_node_constraints.condense (system_rhs);
-}
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::solve ()
-{
- SolverControl solver_control (system_rhs.size(), 1e-6);
- SolverCG<> cg (solver_control);
-
- cg.solve (system_matrix, solution, system_rhs,
- PreconditionIdentity());
-
- hanging_node_constraints.distribute (solution);
-}
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::eval_error (const Function<dim> &f_test)
-{
- hp::QCollection<dim> quadrature(QGauss<dim>(6));
- Vector<double> cellwise_errors (triangulation.n_active_cells());
- VectorTools::integrate_difference (dof_handler, solution, f_test,
- cellwise_errors, quadrature,
- VectorTools::L2_norm);
- const double p_l2_error = cellwise_errors.l2_norm();
-
- deallog << "L2_Error : " << p_l2_error << std::endl;
-}
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::output_results () const
-{
- 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 ();
-
- data_out.write_gnuplot (deallog.get_file_stream());
-}
-
-
-
-template <int dim>
-void LaplaceProblem<dim>::run (const Function<dim> &f_test,
- bool random,
- unsigned int *indx)
-{
- FE_Q<dim> fe_1 (indx[0]),
- fe_2 (indx[1]),
- fe_3 (indx[2]),
- fe_4 (indx[3]);
-
- fe.push_back (fe_1);
- fe.push_back (fe_2);
- fe.push_back (fe_3);
- fe.push_back (fe_4);
-
- make_grid_and_dofs (random);
- assemble_system (f_test);
- solve ();
- eval_error (f_test);
- output_results ();
+ hanging_node_constraints.print (deallog.get_file_stream ());
}
template <int dim>
-void run_test (const Function<dim> &f_test,
- unsigned int *indx)
+void run_test (unsigned int *indx)
{
- LaplaceProblem<dim> laplace_problem_1;
- laplace_problem_1.run (f_test, true, indx);
-
- LaplaceProblem<dim> laplace_problem_2;
- laplace_problem_2.run (f_test, false, indx);
+ run<dim> (true, indx);
+ run<dim> (false, indx);
}
};
- Linear<2> test2d_1;
- Quadratic<2> test2d_2;
- Cubic<2> test2d_3;
-
- Linear<3> test3d_1;
- Quadratic<3> test3d_2;
- Cubic<3> test3d_3;
-
deallog << "Testing Order 1" << std::endl;
- run_test<2> (test2d_1, &(index[0]));
- run_test<3> (test3d_1, &(index[0]));
+ run_test<2> (&(index[0]));
+ run_test<3> (&(index[0]));
deallog << "Testing Order 2" << std::endl;
- run_test<2> (test2d_2, &(index[1]));
- run_test<3> (test3d_2, &(index[1]));
+ run_test<2> (&(index[1]));
+ run_test<3> (&(index[1]));
deallog << "Testing Order 3" << std::endl;
- run_test<2> (test2d_3, &(index[2]));
- run_test<3> (test3d_3, &(index[2]));
+ run_test<2> (&(index[2]));
+ run_test<3> (&(index[2]));
return 0;
}