From: Wolfgang Bangerth Date: Mon, 25 Sep 2006 21:25:21 +0000 (+0000) Subject: Modify this testcase to only check the hanging nodes the test's name suggest. Correct... X-Git-Tag: v8.0.0~11013 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=aacd81d9434e5e5ee4b21ea61fb206d40e78e790;p=dealii.git Modify this testcase to only check the hanging nodes the test's name suggest. Correctness changes will be tested in a separate set of tests based on the deal.II/interpolate_* and deal.II/project_* tests at one point. git-svn-id: https://svn.dealii.org/trunk@13960 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/hp/hp_hanging_nodes_02.cc b/tests/hp/hp_hanging_nodes_02.cc index 28ac7dc61a..f1c365322b 100644 --- a/tests/hp/hp_hanging_nodes_02.cc +++ b/tests/hp/hp_hanging_nodes_02.cc @@ -12,11 +12,9 @@ //---------------------------- 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 #include @@ -53,121 +51,24 @@ std::ofstream logfile("hp_hanging_nodes_02/output"); template -class Linear : public Function +void run (bool random_p, + unsigned int *indx) { - public: - Linear () : Function() {}; - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - - - -template -double Linear::value (const Point &p, - const unsigned int) const -{ - double res = p(0); - for (unsigned int d = 1; d < dim; ++d) - res *= p(d); - - return res; -} - - - -template -class Quadratic : public Function -{ - public: - Quadratic () : Function() {}; - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - - - -template -double Quadratic::value (const Point &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 -class Cubic : public Function -{ - public: - Cubic () : Function() {}; - - virtual double value (const Point &p, - const unsigned int component = 0) const; -}; - + Triangulation triangulation; + hp::FECollection fe; + hp::DoFHandler dof_handler(triangulation); + ConstraintMatrix hanging_node_constraints; + FE_Q fe_1 (indx[0]), + fe_2 (indx[1]), + fe_3 (indx[2]), + fe_4 (indx[3]); -template -double Cubic::value (const Point &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 -class LaplaceProblem -{ - public: - LaplaceProblem (); - - void run (const Function &f_test, - bool random, - unsigned int *indx); - - private: - void make_grid_and_dofs (const bool random_p); - void assemble_system (const Function &f_test); - void solve (); - void eval_error (const Function &f_test); - void output_results () const; - - Triangulation triangulation; - hp::FECollection fe; - hp::DoFHandler dof_handler; - - SparsityPattern sparsity_pattern; - SparseMatrix system_matrix; - - // Although we do not have h-refinement, - // hanging nodes will inevitably appear - // due to different polynomial degrees. - ConstraintMatrix hanging_node_constraints; - - Vector solution; - Vector system_rhs; -}; - - -template -LaplaceProblem::LaplaceProblem () : - dof_handler (triangulation) -{} - - -template -void LaplaceProblem::make_grid_and_dofs (const bool random_p) -{ GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (5-dim); deallog << "Number of active cells: " @@ -207,15 +108,6 @@ void LaplaceProblem::make_grid_and_dofs (const bool random_p) << 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. @@ -226,155 +118,15 @@ void LaplaceProblem::make_grid_and_dofs (const bool random_p) 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 -void LaplaceProblem::assemble_system (const Function &f_test) -{ - hp::QCollection quadrature_formula(QGauss(6)); - hp::FEValues 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 cell_matrix (max_dofs_per_cell, max_dofs_per_cell); - Vector cell_rhs (max_dofs_per_cell); - - std::vector local_dof_indices (max_dofs_per_cell); - - typename hp::DoFHandler::active_cell_iterator - cell = dof_handler.begin_active(), - endc = dof_handler.end(); - for (; cell!=endc; ++cell) - { - x_fe_values.reinit (cell); - - const FEValues &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 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 -void LaplaceProblem::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 -void LaplaceProblem::eval_error (const Function &f_test) -{ - hp::QCollection quadrature(QGauss(6)); - Vector 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 -void LaplaceProblem::output_results () const -{ - DataOut > 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 -void LaplaceProblem::run (const Function &f_test, - bool random, - unsigned int *indx) -{ - FE_Q 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 -void run_test (const Function &f_test, - unsigned int *indx) +void run_test (unsigned int *indx) { - LaplaceProblem laplace_problem_1; - laplace_problem_1.run (f_test, true, indx); - - LaplaceProblem laplace_problem_2; - laplace_problem_2.run (f_test, false, indx); + run (true, indx); + run (false, indx); } @@ -393,25 +145,17 @@ int main () }; - 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; }