From d734f025b4dc24c5125d9074fa4b8d6c30595309 Mon Sep 17 00:00:00 2001 From: goll Date: Thu, 26 Jan 2012 09:46:44 +0000 Subject: [PATCH] Tests to test point_value in the hp setting. git-svn-id: https://svn.dealii.org/trunk@24930 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/bits/point_value_hp_01.cc | 183 +++++++++++++++++++++ tests/bits/point_value_hp_01/cmp/generic | 22 +++ tests/bits/point_value_hp_02.cc | 201 +++++++++++++++++++++++ tests/bits/point_value_hp_02/cmp/generic | 40 +++++ 4 files changed, 446 insertions(+) create mode 100644 tests/bits/point_value_hp_01.cc create mode 100644 tests/bits/point_value_hp_01/cmp/generic create mode 100644 tests/bits/point_value_hp_02.cc create mode 100644 tests/bits/point_value_hp_02/cmp/generic diff --git a/tests/bits/point_value_hp_01.cc b/tests/bits/point_value_hp_01.cc new file mode 100644 index 0000000000..aebfe10c07 --- /dev/null +++ b/tests/bits/point_value_hp_01.cc @@ -0,0 +1,183 @@ +//---------------------------- point_value_hp_01.cc --------------------------- +// Copyright (C) 2004, 2006, 2007 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. +// +//---------------------------- point_value_hp_01.cc --------------------------- + +// check VectorTools::point_value for a hp dofhandler. + + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +template +class MySquareFunction : public Function +{ + public: + MySquareFunction () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const + { return (component+1)*p.square()+1; } + + virtual void vector_value (const Point &p, + Vector &values) const + { values(0) = value(p,0); } +}; + + +template +class MyExpFunction : public Function +{ + public: + MyExpFunction () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const + { return std::exp (p(0)); } + + virtual void vector_value (const Point &p, + Vector &values) const + { values(0) = value(p,0); } +}; + + + +template +void make_mesh (Triangulation &tria) +{ + + GridGenerator::hyper_cube(tria, -1, 1); + + // refine the mesh in a random way so as to + // generate as many cells with + // hanging nodes as possible + tria.refine_global (4-dim); + const double steps[4] = { /*d=0*/ 0, 7, 3, 3 }; + for (unsigned int i=0; i::active_cell_iterator + cell = tria.begin_active(); + for (unsigned int index=0; cell != tria.end(); ++cell, ++index) + if (index % (3*dim) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement (); + } +} + + + + +template +void +check () +{ + Triangulation tria; + make_mesh (tria); + + hp::FECollection fe; + fe.push_back (FE_Q(3)); + fe.push_back (FE_Q(4)); + fe.push_back (FE_Q(5)); + + hp::DoFHandler dof_handler (tria); + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell->set_active_fe_index (std::rand() % fe.size()); + } + + dof_handler.distribute_dofs (fe); + + // test with two different functions: one + // that is exactly representable on the + // chosen finite element space, and one + // that isn't + for (unsigned int i=0; i<2; ++i) + { + static const MySquareFunction function_1; + static const Functions::CosineFunction function_2; + + const Function & + function = (i==0 ? + static_cast&>(function_1) : + static_cast&>(function_2)); + + Vector v (dof_handler.n_dofs()); + VectorTools::interpolate (dof_handler, function, v); + + // for the following points, check the + // function value, output it, and + // verify that the value retrieved from + // the interpolated function is close + // enough to that of the real function + // + // also verify that the actual value is + // roughly correct + Point p[3]; + for (unsigned int d=0; d value(1); + for (unsigned int i=0; i<3; ++i) + { + VectorTools::point_value (dof_handler, v, p[i], value); + deallog << -value(0) << std::endl; + + Assert (std::abs(value(0) - function.value(p[i])) < 1e-4, + ExcInternalError()); + + const double scalar_value = VectorTools::point_value (dof_handler, v, p[i]); + Assert (std::abs(value(0) - scalar_value) < 1e-4, + ExcInternalError()); + } + } + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile ("point_value_hp_01/output"); + deallog << std::setprecision (4); + deallog.attach(logfile); + deallog.depth_console (0); + + deallog.push ("1d"); + check<1> (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/bits/point_value_hp_01/cmp/generic b/tests/bits/point_value_hp_01/cmp/generic new file mode 100644 index 0000000000..763c6337fe --- /dev/null +++ b/tests/bits/point_value_hp_01/cmp/generic @@ -0,0 +1,22 @@ + +DEAL:1d::-1.000 +DEAL:1d::-1.250 +DEAL:1d::-1.111 +DEAL:1d::-1.000 +DEAL:1d::-0.7071 +DEAL:1d::-0.8660 +DEAL:1d::OK +DEAL:2d::-1.000 +DEAL:2d::-1.500 +DEAL:2d::-1.222 +DEAL:2d::-1.000 +DEAL:2d::-0.5000 +DEAL:2d::-0.7500 +DEAL:2d::OK +DEAL:3d::-1.000 +DEAL:3d::-1.750 +DEAL:3d::-1.333 +DEAL:3d::-1.000 +DEAL:3d::-0.3536 +DEAL:3d::-0.6495 +DEAL:3d::OK diff --git a/tests/bits/point_value_hp_02.cc b/tests/bits/point_value_hp_02.cc new file mode 100644 index 0000000000..c14af9a4b1 --- /dev/null +++ b/tests/bits/point_value_hp_02.cc @@ -0,0 +1,201 @@ +//---------------------------- point_value_hp_02.cc --------------------------- +// Copyright (C) 2004, 2006, 2007 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. +// +//---------------------------- point_value_hp_02.cc --------------------------- + +// check VectorTools::point_value for hp, alternative algorithm with mapping + + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + + +template +class MySquareFunction : public Function +{ + public: + MySquareFunction () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const + { return (component+1)*p.square()+1; } + + virtual void vector_value (const Point &p, + Vector &values) const + { values(0) = value(p,0); } +}; + + +template +class MyExpFunction : public Function +{ + public: + MyExpFunction () : Function () {} + + virtual double value (const Point &p, + const unsigned int component) const + { return std::exp (p(0)); } + + virtual void vector_value (const Point &p, + Vector &values) const + { values(0) = value(p,0); } +}; + + + +template +void make_mesh (Triangulation &tria) +{ + + GridGenerator::hyper_cube(tria, -1, 1); + + // refine the mesh in a random way so as to + // generate as many cells with + // hanging nodes as possible + tria.refine_global (4-dim); + const double steps[4] = { /*d=0*/ 0, 7, 3, 3 }; + for (unsigned int i=0; i::active_cell_iterator + cell = tria.begin_active(); + for (unsigned int index=0; cell != tria.end(); ++cell, ++index) + if (index % (3*dim) == 0) + cell->set_refine_flag(); + tria.execute_coarsening_and_refinement (); + } +} + + + + +template +void +check () +{ + Triangulation tria; + make_mesh (tria); + + hp::FECollection fe; + fe.push_back (FE_Q(3)); + fe.push_back (FE_Q(4)); + fe.push_back (FE_Q(5)); + + hp::MappingCollection mapping_1; + mapping_1.push_back(MappingQ1()); + mapping_1.push_back(MappingQ1()); + mapping_1.push_back(MappingQ1()); + + hp::MappingCollection mapping_2; + mapping_2.push_back(MappingQ1()); + + hp::DoFHandler dof_handler (tria); + + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell) + { + cell->set_active_fe_index (std::rand() % fe.size()); + } + + dof_handler.distribute_dofs (fe); + + // test with two different functions: one + // that is exactly representable on the + // chosen finite element space, and one + // that isn't + for (unsigned int i=0; i<2; ++i) + { + static const MySquareFunction function_1; + static const Functions::CosineFunction function_2; + + const Function & + function = (i==0 ? + static_cast&>(function_1) : + static_cast&>(function_2)); + + Vector v (dof_handler.n_dofs()); + VectorTools::interpolate (dof_handler, function, v); + + // for the following points, check the + // function value, output it, and + // verify that the value retrieved from + // the interpolated function is close + // enough to that of the real function + // + // also verify that the actual value is + // roughly correct + Point p[3]; + for (unsigned int d=0; d value(1); + for (unsigned int i=0; i<3; ++i) + { + VectorTools::point_value (mapping_1, dof_handler, v, p[i], value); + deallog << -value(0) << std::endl; + + Assert (std::abs(value(0) - function.value(p[i])) < 1e-4, + ExcInternalError()); + VectorTools::point_value (mapping_2, dof_handler, v, p[i], value); + deallog << -value(0) << std::endl; + + Assert (std::abs(value(0) - function.value(p[i])) < 1e-4, + ExcInternalError()); + + const double scalar_value_1 = VectorTools::point_value (mapping_1, dof_handler, v, p[i]); + Assert (std::abs(value(0) - scalar_value_1) < 1e-4, + ExcInternalError()); + + const double scalar_value_2 = VectorTools::point_value (mapping_2, dof_handler, v, p[i]); + Assert (std::abs(value(0) - scalar_value_2) < 1e-4, + ExcInternalError()); + } + } + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile ("point_value_hp_02/output"); + deallog << std::setprecision (4); + deallog.attach(logfile); + deallog.depth_console (0); + + deallog.push ("1d"); + check<1> (); + deallog.pop (); + deallog.push ("2d"); + check<2> (); + deallog.pop (); + deallog.push ("3d"); + check<3> (); + deallog.pop (); +} diff --git a/tests/bits/point_value_hp_02/cmp/generic b/tests/bits/point_value_hp_02/cmp/generic new file mode 100644 index 0000000000..2d85d565af --- /dev/null +++ b/tests/bits/point_value_hp_02/cmp/generic @@ -0,0 +1,40 @@ + +DEAL:1d::-1.000 +DEAL:1d::-1.000 +DEAL:1d::-1.250 +DEAL:1d::-1.250 +DEAL:1d::-1.111 +DEAL:1d::-1.111 +DEAL:1d::-1.000 +DEAL:1d::-1.000 +DEAL:1d::-0.7071 +DEAL:1d::-0.7071 +DEAL:1d::-0.8660 +DEAL:1d::-0.8660 +DEAL:1d::OK +DEAL:2d::-1.000 +DEAL:2d::-1.000 +DEAL:2d::-1.500 +DEAL:2d::-1.500 +DEAL:2d::-1.222 +DEAL:2d::-1.222 +DEAL:2d::-1.000 +DEAL:2d::-1.000 +DEAL:2d::-0.5000 +DEAL:2d::-0.5000 +DEAL:2d::-0.7500 +DEAL:2d::-0.7500 +DEAL:2d::OK +DEAL:3d::-1.000 +DEAL:3d::-1.000 +DEAL:3d::-1.750 +DEAL:3d::-1.750 +DEAL:3d::-1.333 +DEAL:3d::-1.333 +DEAL:3d::-1.000 +DEAL:3d::-1.000 +DEAL:3d::-0.3536 +DEAL:3d::-0.3536 +DEAL:3d::-0.6495 +DEAL:3d::-0.6495 +DEAL:3d::OK -- 2.39.5