From: oliver Date: Thu, 21 Oct 2004 16:25:59 +0000 (+0000) Subject: Added a testcase for the FE_Q class in 3D, which tests the automatically X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=f31ea314d2eb189587333ae14dfa18dfa9ec3f0f;p=dealii-svn.git Added a testcase for the FE_Q class in 3D, which tests the automatically generated constraint matrices. git-svn-id: https://svn.dealii.org/trunk@9718 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/bits/fe_q_constraints.cc b/tests/bits/fe_q_constraints.cc new file mode 100644 index 0000000000..88e5f300fc --- /dev/null +++ b/tests/bits/fe_q_constraints.cc @@ -0,0 +1,272 @@ +//--------------------- fe_q_constraints.cc -------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2003, 2004 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. +// +//--------------------- fe_q_constraints.cc -------------------------- + + +// This file implements a simple test procedure for the algorithm, +// which constructs the constraint matrices for FE_Q<3> elements +// of arbitrary order. After interpolating a polynomial with the +// degree of the ansatzspace, the hanging node constraints are +// distributed to the solution vector. This procedure should not +// change the solution vector. Hence the difference between the +// exact function and the FE function is measured before and +// after the distribution of the hanging node constraints. + +#include "../tests.h" +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +template +class TestFunction : public Function +{ +private: + std::vector > base; + const unsigned int p_order; + +public: + TestFunction (const unsigned int p_order); + + virtual double value (const Point &p, + const unsigned int component = 0) const; +}; + + +template +TestFunction::TestFunction (const unsigned int p_order) : + p_order (p_order) +{ + std::vector coeff(p_order); + + for (unsigned int d = 0; d < dim; ++d) + { + for (unsigned int c = 0; c < p_order; ++c) + coeff[c] = (double) rand () / (double) RAND_MAX; + base.push_back (Polynomials::Polynomial (coeff)); + } +} + + +template +double TestFunction::value (const Point &p, + const unsigned int /*component*/) const +{ + double val = base[0].value(p(0)); + for (unsigned int i = 1; i < dim; ++i) + val *= base[i].value(p(i)); + return val; +} + + +template +class TestFEQConstraints +{ + public: + TestFEQConstraints (unsigned int p_order, unsigned int refinements); + void run (); + + private: + const unsigned int refinements; + const unsigned int p_order; + + void refine_grid_random (); + void make_grid_and_dofs (); + void output_grid () const; + + void test (); + + Triangulation triangulation; + FE_Q fe; + DoFHandler dof_handler; + ConstraintMatrix hanging_node_constraints; +}; + +template +TestFEQConstraints::TestFEQConstraints (unsigned int p_order, + unsigned int refinements) : + refinements (refinements), + p_order (p_order), + fe (p_order), + dof_handler (triangulation) +{ +} + + +// Actually this function creates with pseudo-random numbers errors, +// which are then used to refine the grid. This should be a tough test for +// the hanging nodes. +template +void TestFEQConstraints::refine_grid_random () +{ + const unsigned int n_cells = triangulation.n_active_cells(); + Vector estimated_error_per_cell (n_cells); + + for (unsigned int i = 0; i < n_cells; ++i) + estimated_error_per_cell(i) = (double) rand () / (double) RAND_MAX; + + GridRefinement::refine_and_coarsen_fixed_number (triangulation, + estimated_error_per_cell, + 0.3, 0.03); + triangulation.execute_coarsening_and_refinement (); +} + + +template +void TestFEQConstraints::make_grid_and_dofs () +{ + GridGenerator::hyper_cube (triangulation, -1, 1); + triangulation.refine_global (2); + + for (unsigned int i = 0; i < refinements; ++i) + refine_grid_random (); + + dof_handler.distribute_dofs (fe); + + deallog << "---------------------------------------------------------" + << std::endl; + deallog << "P-Order: " << p_order + << " Number of degrees of freedom: " + << dof_handler.n_dofs() + << std::endl; + + hanging_node_constraints.clear (); + DoFTools::make_hanging_node_constraints (dof_handler, + hanging_node_constraints); + hanging_node_constraints.close (); +} + + +template +void TestFEQConstraints::output_grid () const +{ + std::ofstream output ("test_feq3.eps"); + GridOut grid_out; + grid_out.write_eps (triangulation, output); +} + + +template +void TestFEQConstraints::test () +{ + TestFunction test_function (p_order); + double l2test, + l2norm1, + l2norm2, + l2error; + + Vector solution; + solution.reinit (dof_handler.n_dofs ()); + hanging_node_constraints.distribute (solution); + + QGauss quadrature (p_order + 1); + + Vector norm_per_cell (triangulation.n_active_cells ()); + VectorTools::interpolate (dof_handler, test_function, solution); + + // First error check. Simply the interpolation error of the used FE-Space + // on the given triangulation. + VectorTools::integrate_difference (dof_handler, solution, + ZeroFunction(1), + norm_per_cell, + quadrature, + VectorTools::L2_norm); + l2test = norm_per_cell.l2_norm (); + deallog << "L2-Norm of test function " + << l2test << std::endl; + + // First error check. Simply the interpolation error of the used FE-Space + // on the given triangulation. + VectorTools::integrate_difference (dof_handler, solution, + test_function, + norm_per_cell, + quadrature, + VectorTools::L2_norm); + l2norm1 = norm_per_cell.l2_norm () / l2test; + + // Second error check. Interpolation error, after redistribution of the + // values onto the DoFs on the hanging nodes. + hanging_node_constraints.distribute (solution); + VectorTools::integrate_difference (dof_handler, solution, + test_function, + norm_per_cell, + quadrature, + VectorTools::L2_norm); + l2norm2 = norm_per_cell.l2_norm () / l2test; + l2error = fabs (l2norm1 - l2norm2); + deallog << "Normed L2-Error 1: " << l2norm1 + << " Normed L2-Error 2: " << l2norm2 << std::endl + << "Normed L2-Diff: " << l2error << " "; + if (l2error < 1.0e-16) + deallog << "OK" << std::endl; + else + deallog << "Error !" << std::endl; +} + + +template +void TestFEQConstraints::run () +{ + make_grid_and_dofs (); +// output_grid (); + test (); +} + + +int main () +{ + std::ofstream logfile("fe_q_constraints.output"); + deallog.attach(logfile); + deallog.depth_console(0); + + unsigned int ref_level[] = {5, 4, 3, 3, 2, 2}; + + // Initialise the random generator with an arbitrary number. This + // should ensure reproducible results. + srand (4375384); + + // With these parameters, the elements are tested up to a polynomial + // degree of 6. The triangulation is refined 2 times uniformly and + // several times (pseudo-) randomly. + TestFEQConstraints<3> *test; + for (unsigned int p = 1; p < 7; ++p) + { + test = new TestFEQConstraints<3> (p, ref_level[p-1]); + test->run (); + delete test; + } + + return 0; +}