From: bangerth Date: Fri, 11 Aug 2006 17:11:55 +0000 (+0000) Subject: Add tests to check that the projection of an element yields the correct approximation. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=d30dca4fc2f7a968c36f8826ef954d19dfc0af0e;p=dealii-svn.git Add tests to check that the projection of an element yields the correct approximation. git-svn-id: https://svn.dealii.org/trunk@13672 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/tests/deal.II/Makefile b/tests/deal.II/Makefile index 3fdea2819a..8a469d27bf 100644 --- a/tests/deal.II/Makefile +++ b/tests/deal.II/Makefile @@ -54,7 +54,7 @@ tests_x = block_matrices \ interpolate_* \ have_same_coarse_mesh_* \ get_finest_common_cells_* \ - project_* + project_*[0-9] # from above list of regular expressions, generate the real set of # tests diff --git a/tests/deal.II/project_common.cc b/tests/deal.II/project_common.cc new file mode 100644 index 0000000000..faebd59cf6 --- /dev/null +++ b/tests/deal.II/project_common.cc @@ -0,0 +1,156 @@ +//---------------------------- project_common.cc --------------------------- +// $Id: project_common.cc 12732 2006-03-28 23:15:45Z wolf $ +// Version: $Name$ +// +// Copyright (C) 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. +// +//---------------------------- project_common.cc --------------------------- + + +// common framework to check whether an element of polynomial order p can +// represent functions of order q + +#include "../tests.h" +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +template +void test (); + + +template +class F : public Function +{ + public: + F (const unsigned int q, + const unsigned int n_components) + : + Function(n_components), + q(q) + {} + + virtual double value (const Point &p, + const unsigned int component) const + { + Assert ((component == 0) && (this->n_components == 1), + ExcInternalError()); + double val = 0; + for (unsigned int d=0; d &p, + Vector &v) const + { + for (unsigned int c=0; c +void test_no_hanging_nodes (const FiniteElement &fe, + const unsigned int p) +{ + Triangulation triangulation; + GridGenerator::hyper_cube (triangulation); + triangulation.refine_global (3); + + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs (fe); + + deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl; + + // there are no constraints here + ConstraintMatrix constraints; + constraints.close (); + + Vector interpolant (dof_handler.n_dofs()); + Vector error (triangulation.n_active_cells()); + for (unsigned int q=0; q<=p+2; ++q) + { + // project the function + VectorTools::project (dof_handler, + constraints, + QGauss(p+2), + F (q, fe.n_components()), + interpolant); + + // then compute the interpolation error + VectorTools::integrate_difference (dof_handler, + interpolant, + F (q, fe.n_components()), + error, + QGauss(std::max(p,q)+1), + VectorTools::L2_norm); + deallog << fe.get_name() << ", P_" << q + << ", rel. error=" << error.l2_norm() / interpolant.l2_norm() + << std::endl; + + if (q<=p) + Assert (error.l2_norm() < 1e-12*interpolant.l2_norm(), + ExcInternalError()); + } +} + + + +int main () +{ + std::ofstream logfile(logname); + logfile.precision (3); + + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + test<1>(); + test<2>(); + test<3>(); +} + diff --git a/tests/deal.II/project_q_01.cc b/tests/deal.II/project_q_01.cc new file mode 100644 index 0000000000..80042c6d2f --- /dev/null +++ b/tests/deal.II/project_q_01.cc @@ -0,0 +1,29 @@ +//---------------------------- project_q_01.cc --------------------------- +// $Id: project_q_01.cc 12732 2006-03-28 23:15:45Z wolf $ +// Version: $Name$ +// +// Copyright (C) 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. +// +//---------------------------- project_q_01.cc --------------------------- + + +// check that VectorTools::project works for Q elements correctly on +// a uniformly refined mesh for functions of degree q + +char logname[] = "project_q_01/output"; + + +#include "project_common.cc" + + +template +void test () +{ + for (unsigned int p=1; p<6-dim; ++p) + test_no_hanging_nodes (FE_Q(p), p); +} diff --git a/tests/deal.II/project_q_01/cmp/generic b/tests/deal.II/project_q_01/cmp/generic new file mode 100644 index 0000000000..35cf8fe9c1 --- /dev/null +++ b/tests/deal.II/project_q_01/cmp/generic @@ -0,0 +1,47 @@ + +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], P_2, rel. error=0.000278 +DEAL::FESystem<1>[FE_Q<1>(1)^2-FE_Q<1>(2)], P_3, rel. error=0.000654 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_3, rel. error=5.15e-06 +DEAL::FESystem<1>[FE_Q<1>(2)^2-FE_Q<1>(3)], P_4, rel. error=1.52e-05 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_4, rel. error=1.20e-07 +DEAL::FESystem<1>[FE_Q<1>(3)^2-FE_Q<1>(4)], P_5, rel. error=4.27e-07 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_0, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_1, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_2, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_3, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_4, rel. error=0 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_5, rel. error=3.21e-09 +DEAL::FESystem<1>[FE_Q<1>(4)^2-FE_Q<1>(5)], P_6, rel. error=1.33e-08 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], P_2, rel. error=8.74e-05 +DEAL::FESystem<2>[FE_Q<2>(1)^2-FE_Q<2>(2)], P_3, rel. error=0.000201 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_3, rel. error=1.02e-06 +DEAL::FESystem<2>[FE_Q<2>(2)^2-FE_Q<2>(3)], P_4, rel. error=3.01e-06 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_0, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_1, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_2, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_3, rel. error=0 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_4, rel. error=2.22e-08 +DEAL::FESystem<2>[FE_Q<2>(3)^2-FE_Q<2>(4)], P_5, rel. error=7.82e-08 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], P_2, rel. error=2.52e-05 +DEAL::FESystem<3>[FE_Q<3>(1)^2-FE_Q<3>(2)], P_3, rel. error=5.74e-05 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_0, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_1, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_2, rel. error=0 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_3, rel. error=2.04e-07 +DEAL::FESystem<3>[FE_Q<3>(2)^2-FE_Q<3>(3)], P_4, rel. error=6.07e-07 diff --git a/tests/deal.II/project_rt_01.cc b/tests/deal.II/project_rt_01.cc index b2b2040226..abd3beb019 100644 --- a/tests/deal.II/project_rt_01.cc +++ b/tests/deal.II/project_rt_01.cc @@ -12,121 +12,19 @@ //---------------------------- project_rt_01.cc --------------------------- -// check that VectorTools::project works for RT elements correctly on +// check that VectorTools::project works for Q elements correctly on // a uniformly refined mesh for functions of degree q -#include "../tests.h" -#include -#include -#include -#include +char logname[] = "project_rt_01/output"; -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include - - -template -class F : public Function -{ - public: - F (const unsigned int q, - const unsigned int n_components) - : - Function(n_components), - q(q) - {} - - virtual void vector_value (const Point &p, - Vector &v) const - { - for (unsigned int c=0; c void test () { - Triangulation triangulation; - GridGenerator::hyper_cube (triangulation); - triangulation.refine_global (3); - - for (unsigned int p=0; p<6-dim; ++p) - { - FE_RaviartThomas fe(p); - DoFHandler dof_handler(triangulation); - dof_handler.distribute_dofs (fe); - - deallog << "n_dofs=" << dof_handler.n_dofs() << std::endl; - - ConstraintMatrix constraints; - constraints.close (); - - Vector interpolant (dof_handler.n_dofs()); - Vector error (triangulation.n_active_cells()); - for (unsigned int q=0; q<=p+2; ++q) - { - // project the function - VectorTools::project (dof_handler, - constraints, - QGauss(p+2), - F (q, fe.n_components()), - interpolant); - - // then compute the interpolation error - VectorTools::integrate_difference (dof_handler, - interpolant, - F (q, fe.n_components()), - error, - QGauss(std::max(p,q)+1), - VectorTools::L2_norm); - deallog << fe.get_name() << ", P_" << q - << ", rel. error=" << error.l2_norm() / interpolant.l2_norm() - << std::endl; - - if (q<=p) - Assert (error.l2_norm() < 1e-12*interpolant.l2_norm(), - ExcInternalError()); - } - } -} - - - -int main () -{ - std::ofstream logfile("project_rt_01/output"); - logfile.precision (3); - - deallog.attach(logfile); - deallog.depth_console(0); - deallog.threshold_double(1.e-10); - - test<2>(); - test<3>(); + if (dim != 1) + for (unsigned int p=0; p<6-dim; ++p) + test_no_hanging_nodes (FE_RaviartThomas(p), p); } -