From: kayser-herold Date: Sat, 29 Jul 2006 23:43:28 +0000 (+0000) Subject: New test, which L_2 projects a polynomial function X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8a5372e139601ca5e36d8caf7354613c38a4decb;p=dealii-svn.git New test, which L_2 projects a polynomial function to a hp FE discretisation and afterwards evaluates the L_2 norm of the error. The L_2 norm of the error should be in the range of the FP accuracy. git-svn-id: https://svn.dealii.org/trunk@13513 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/hp/hp_hanging_nodes_02.cc b/tests/hp/hp_hanging_nodes_02.cc new file mode 100644 index 0000000000..e41eff51db --- /dev/null +++ b/tests/hp/hp_hanging_nodes_02.cc @@ -0,0 +1,417 @@ +//---------------------------- step-3.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-3.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. + +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +#include + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include + +#include +#include +#include + +std::ofstream logfile("hp_hanging_nodes_02/output"); + + +template +class Linear : public Function +{ + 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; +}; + + + +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); + + 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 (3); + deallog << "Number of active cells: " + << triangulation.n_active_cells() + << std::endl; + deallog << "Total number of cells: " + << triangulation.n_cells() + << std::endl; + + // Now to the p-Method. Assign + // random active_fe_indices to the + // different cells. + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active (), + endc = dof_handler.end (); + if (random_p) + { + for (; cell != endc; ++cell) + { + cell->set_active_fe_index ((int)(4.0 * (double) random () / (double) RAND_MAX)); + } + } + else + { + unsigned int cell_no = 0; + for (; cell != endc; ++cell) + { + if (cell_no >= triangulation.n_active_cells () / 2) + cell->set_active_fe_index (1); + else + cell->set_active_fe_index (0); + } + } + + + dof_handler.distribute_dofs (fe); + deallog << "Number of degrees of freedom: " + << 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.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + + 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 (1000, 1e-12); + 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 (); +} + + +template +void run_test (const Function &f_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); +} + + + +int main () +{ + logfile.precision(2); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + unsigned int index[] = + { + 1,2,3,4,5,6,7 + }; + + + 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])); + + deallog << "Testing Order 2" << std::endl; + run_test<2> (test2d_2, &(index[1])); + run_test<3> (test3d_2, &(index[1])); + + deallog << "Testing Order 3" << std::endl; + run_test<2> (test2d_3, &(index[2])); + run_test<3> (test3d_3, &(index[2])); + + return 0; +}