From 67d39c1392e3a6e4c80097c0b3c344a30bd14528 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 4 Jul 2012 12:48:34 +0000 Subject: [PATCH] Several more tests. git-svn-id: https://svn.dealii.org/trunk@25679 0785d39b-7218-0410-832d-ea1e28bc413d --- tests/bits/fe_field_function_01_vector.cc | 110 +++++++++++++++ .../fe_field_function_01_vector/cmp/generic | 16 +++ tests/bits/fe_field_function_02_vector.cc | 109 +++++++++++++++ .../fe_field_function_02_vector/cmp/generic | 10 ++ tests/bits/fe_field_function_03_vector.cc | 88 ++++++++++++ .../fe_field_function_03_vector/cmp/generic | 19 +++ tests/bits/fe_field_function_04_vector.cc | 115 ++++++++++++++++ .../fe_field_function_04_vector/cmp/generic | 2 + tests/bits/fe_field_function_05_vector.cc | 127 ++++++++++++++++++ .../fe_field_function_05_vector/cmp/generic | 2 + 10 files changed, 598 insertions(+) create mode 100644 tests/bits/fe_field_function_01_vector.cc create mode 100644 tests/bits/fe_field_function_01_vector/cmp/generic create mode 100644 tests/bits/fe_field_function_02_vector.cc create mode 100644 tests/bits/fe_field_function_02_vector/cmp/generic create mode 100644 tests/bits/fe_field_function_03_vector.cc create mode 100644 tests/bits/fe_field_function_03_vector/cmp/generic create mode 100644 tests/bits/fe_field_function_04_vector.cc create mode 100644 tests/bits/fe_field_function_04_vector/cmp/generic create mode 100644 tests/bits/fe_field_function_05_vector.cc create mode 100644 tests/bits/fe_field_function_05_vector/cmp/generic diff --git a/tests/bits/fe_field_function_01_vector.cc b/tests/bits/fe_field_function_01_vector.cc new file mode 100644 index 0000000000..391a69923a --- /dev/null +++ b/tests/bits/fe_field_function_01_vector.cc @@ -0,0 +1,110 @@ +//---------------------------- template.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2011, 2012 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. +// +//---------------------------- template.cc --------------------------- + + +// Test the FEFieldFunction class + +#include "../tests.h" +#include + +// all include files you need here +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +class F : public Function +{ + public: + F() : Function(2) {} + virtual void vector_value (const Point &p, + Vector &v) const + { + v = 0; + v[0] = p.square(); + } +}; + +double abs_zero(double a) +{ + if( std::abs(a) < 1e-10) + return 0; + else + return a; +} + +template +void test() { + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(8/dim); + + FESystem fe(FE_Q (2),2); + DoFHandler dh(tria); + + dh.distribute_dofs(fe); + deallog << "Ndofs :" << dh.n_dofs() << std::endl; + + + + Vector v1(dh.n_dofs()), v2(dh.n_dofs()); + + VectorTools::interpolate(dh, F(), v1); + deallog << "V norm: " << v1.l2_norm() << std::endl; + + Functions::FEFieldFunction, Vector > + fef(dh, v1); + + VectorTools::interpolate(dh, fef, v2); + + v2.add(-1, v1); + deallog << "Interpolation error: " << abs_zero(v2.l2_norm()) + << std::endl; + + Vector error(tria.n_active_cells()); + QGauss quad(2); + VectorTools::integrate_difference(dh, v1, F(), + error, quad, + VectorTools::L2_norm); + deallog << "L2 Interpolation error: " + << abs_zero(error.l2_norm()) << std::endl; + error = 0; + VectorTools::integrate_difference(dh, v1, fef, + error, quad, + VectorTools::L2_norm); + deallog << "L2 Interpolation error with fef: " + << abs_zero(error.l2_norm()) << std::endl; + +} + +int main () +{ + std::ofstream logfile("fe_field_function_01_vector/output"); + deallog.attach(logfile); + deallog.depth_console(0); + + test<1>(); + test<2>(); + test<3>(); + + return 0; +} + diff --git a/tests/bits/fe_field_function_01_vector/cmp/generic b/tests/bits/fe_field_function_01_vector/cmp/generic new file mode 100644 index 0000000000..3115fec483 --- /dev/null +++ b/tests/bits/fe_field_function_01_vector/cmp/generic @@ -0,0 +1,16 @@ + +DEAL::Ndofs :1026 +DEAL::V norm: 10.1440 +DEAL::Interpolation error: 0.00000 +DEAL::L2 Interpolation error: 0.00000 +DEAL::L2 Interpolation error with fef: 0.00000 +DEAL::Ndofs :2178 +DEAL::V norm: 26.5652 +DEAL::Interpolation error: 0.00000 +DEAL::L2 Interpolation error: 0.00000 +DEAL::L2 Interpolation error with fef: 0.00000 +DEAL::Ndofs :1458 +DEAL::V norm: 32.6964 +DEAL::Interpolation error: 0.00000 +DEAL::L2 Interpolation error: 0.00000 +DEAL::L2 Interpolation error with fef: 0.00000 diff --git a/tests/bits/fe_field_function_02_vector.cc b/tests/bits/fe_field_function_02_vector.cc new file mode 100644 index 0000000000..c5db68e84f --- /dev/null +++ b/tests/bits/fe_field_function_02_vector.cc @@ -0,0 +1,109 @@ +//---------------------------- template.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2005, 2011, 2012 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. +// +//---------------------------- template.cc --------------------------- + + +// Test the FEFieldFunction class. This class was not thread-safe at one point +// because it keeps a cache on the side that was invalidated when different +// threads kept pouncing on it. Patrick Sodre wrote a fix for that that's +// based on a thread-local variable instead of the single cache. + +#include "../tests.h" +#include + +// all include files you need here +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +class F : public Function +{ + public: + F() : Function(2) {} + virtual void vector_value (const Point &p, + Vector &v) const + { + v = 0; + v[0] = p.square(); + } +}; + + +double abs_zero(double a) +{ + if( std::abs(a) < 1e-10) + return 0; + else + return a; +} + +template +void test() { + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(8/dim); + + FESystem fe(FE_Q (2),2); + DoFHandler dh(tria); + + dh.distribute_dofs(fe); + deallog << "Ndofs :" << dh.n_dofs() << std::endl; + + F ff; + + Vector v1(dh.n_dofs()), v2(dh.n_dofs()); + + VectorTools::interpolate(dh, ff, v1); + deallog << "V norm: " << v1.l2_norm() << std::endl; + + Functions::FEFieldFunction, Vector > + fef(dh, v1); + + // project the discrete function fef back + // onto the finite element space. this + // should be the identity operation and + // consequently subtracting the vector from + // itself should yield zero + { + ConstraintMatrix cm; + cm.close(); + VectorTools::project(dh, cm, QGauss(3), fef, v2); + } + v2.add(-1, v1); + + deallog << "Projection error: " << abs_zero(v2.l2_norm()) + << std::endl; + Assert (v2.l2_norm() < 1e-10, ExcInternalError()); +} + +int main () +{ + std::ofstream logfile("fe_field_function_02_vector/output"); + deallog.attach(logfile); + deallog.depth_console(0); + + test<1>(); + test<2>(); + test<3>(); + + return 0; +} + diff --git a/tests/bits/fe_field_function_02_vector/cmp/generic b/tests/bits/fe_field_function_02_vector/cmp/generic new file mode 100644 index 0000000000..021c0abe7b --- /dev/null +++ b/tests/bits/fe_field_function_02_vector/cmp/generic @@ -0,0 +1,10 @@ + +DEAL::Ndofs :1026 +DEAL::V norm: 10.1440 +DEAL::Projection error: 0.00000 +DEAL::Ndofs :2178 +DEAL::V norm: 26.5652 +DEAL::Projection error: 0.00000 +DEAL::Ndofs :1458 +DEAL::V norm: 32.6964 +DEAL::Projection error: 0.00000 diff --git a/tests/bits/fe_field_function_03_vector.cc b/tests/bits/fe_field_function_03_vector.cc new file mode 100644 index 0000000000..834205d8be --- /dev/null +++ b/tests/bits/fe_field_function_03_vector.cc @@ -0,0 +1,88 @@ +//------------------------------------------------------- +// $Id: fe_field_function_03.cc $ +// Version: $Name$ +// +// Copyright (C) 2005, 2011 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. +// +//------------------------------------------------------- + + +// Test the functionality of the laplacian in the FEFieldFunction class. + +#include "../tests.h" +#include + +// all include files you need here +#include +#include +#include +#include +#include +#include +#include +#include +#include + +template +void test() { + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(9/dim); + + FE_Q fe(2); + DoFHandler dh(tria); + + dh.distribute_dofs(fe); + deallog << "Ndofs :" << dh.n_dofs() << std::endl; + + Functions::CosineFunction ff; + + Vector v1(dh.n_dofs()), v2(dh.n_dofs()); + + VectorTools::interpolate(dh, ff, v1); + deallog << "V norm: " << v1.l2_norm() << std::endl; + + Functions::FEFieldFunction, Vector > + fef(dh, v1); + + //create the origin + Point p; + //compute the error of the laplacian in this point + deallog << "Value of the laplacian in 0:" << std::endl; + deallog << "correct value: " << ff.laplacian(p) <<", approximation: "<< fef.laplacian(p) << std::endl; + + //now we want to test the list version + Point p1 = Point::unit_vector(0); + p1 = p1 * 0.5; + Point p2 = p1 * 0.5; + std::vector > vec; + vec.push_back(p1); + vec.push_back(p2); + std::vector values_c(2); + std::vector values_a(2); + + //get the laplacians at these two points + ff.laplacian_list(vec, values_c); + fef.laplacian_list(vec, values_a); + deallog << "Value of the laplacian in 0.5*e1 and 0.25 * e1:" << std::endl; + deallog << " correct values: " <(); + test<2>(); + test<3>(); + + return 0; +} + diff --git a/tests/bits/fe_field_function_03_vector/cmp/generic b/tests/bits/fe_field_function_03_vector/cmp/generic new file mode 100644 index 0000000000..cd95fc4b27 --- /dev/null +++ b/tests/bits/fe_field_function_03_vector/cmp/generic @@ -0,0 +1,19 @@ + +DEAL::Ndofs :1025 +DEAL::V norm: 22.6385 +DEAL::Value of the laplacian in 0: +DEAL::correct value: -2.46740, approximation: -2.46743 +DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1: +DEAL:: correct values: -1.74472 -2.27958, approximations: -1.74739 -2.28103 +DEAL::Ndofs :1089 +DEAL::V norm: 16.5000 +DEAL::Value of the laplacian in 0: +DEAL::correct value: -4.93480, approximation: -4.92787 +DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1: +DEAL:: correct values: -3.48943 -4.55916, approximations: -3.57012 -4.59908 +DEAL::Ndofs :4913 +DEAL::V norm: 24.7815 +DEAL::Value of the laplacian in 0: +DEAL::correct value: -7.40220, approximation: -7.36064 +DEAL::Value of the laplacian in 0.5*e1 and 0.25 * e1: +DEAL:: correct values: -5.23415 -6.83874, approximations: -5.37564 -6.89283 diff --git a/tests/bits/fe_field_function_04_vector.cc b/tests/bits/fe_field_function_04_vector.cc new file mode 100644 index 0000000000..3586f30fc6 --- /dev/null +++ b/tests/bits/fe_field_function_04_vector.cc @@ -0,0 +1,115 @@ +//------------------------------------------------------- +// $Id: fe_field_function_04_vector.cc $ +// Version: $Name$ +// +// Copyright (C) 2005, 2011, 2012 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. +// +//------------------------------------------------------- + + +// FEFieldFunction ran into an assertion after +// Mapping::transform_real_to_unit_cell started throwing exceptions +// when it couldn't find the point on the reference cell that belongs +// to a given point, rather than just silently giving up + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +class F : public Function +{ + public: + F() : Function(2) {} + virtual void vector_value (const Point &p, + Vector &v) const + { + v = 0; + v[0] = p.square(); + } +}; + + +template +void test() +{ + const HyperBallBoundary boundary_description; + + Triangulation triangulation; + GridGenerator::hyper_ball (triangulation); + triangulation.set_boundary (0, boundary_description); + triangulation.refine_global (1); + + FESystem fe(FE_Q (2),2); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + Vector solution(dof_handler.n_dofs()); + VectorTools::interpolate (dof_handler, F(), solution); + + Functions::FEFieldFunction<2> fe_function (dof_handler, solution); + std::vector > points; + + // add a bunch of points. all + // points are inside the circle. + // the problem happens because we + // walk over a bunch of cells in + // the process of finding all of + // these points and then realize + // when we get to the one at the + // end that the coordinates for + // this point can't be found in the + // cell we have touched last (it's + // too far away from that cell, and + // the inverse mapping does not + // converge + for (unsigned int i=0; i<20; ++i) + for (unsigned int j=0; j<20; ++j) + points.push_back (Point(-0.7+i*0.07,-0.7+j*0.07)); + points.push_back (Point(-0.27999999999999992, -0.62999999999999989)); + + std::vector > m (points.size()); + fe_function.vector_value_list (points, m); + + for (unsigned int i=0; i(); + + return 0; +} + diff --git a/tests/bits/fe_field_function_04_vector/cmp/generic b/tests/bits/fe_field_function_04_vector/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/bits/fe_field_function_04_vector/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/bits/fe_field_function_05_vector.cc b/tests/bits/fe_field_function_05_vector.cc new file mode 100644 index 0000000000..5a1dde289c --- /dev/null +++ b/tests/bits/fe_field_function_05_vector.cc @@ -0,0 +1,127 @@ +//------------------------------------------------------- +// $Id: fe_field_function_05_vector.cc $ +// Version: $Name$ +// +// Copyright (C) 2005, 2011, 2012 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. +// +//------------------------------------------------------- + + +// FEFieldFunction ran into an assertion after +// Mapping::transform_real_to_unit_cell started throwing exceptions +// when it couldn't find the point on the reference cell that belongs +// to a given point, rather than just silently giving up +// +// this is a variant of _04 found by inspecting the code in +// FEFieldFunction but the (same) exception is triggered in a +// different place + +#include "../tests.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +template +class F : public Function +{ + public: + F() : Function(2) {} + virtual void vector_value (const Point &p, + Vector &v) const + { + v = 0; + v[0] = p.square(); + } +}; + + + +template +void test() +{ + const HyperBallBoundary boundary_description; + + Triangulation triangulation; + GridGenerator::hyper_ball (triangulation); + triangulation.set_boundary (0, boundary_description); + triangulation.refine_global (1); + + FESystem fe(FE_Q (2),2); + DoFHandler dof_handler(triangulation); + dof_handler.distribute_dofs(fe); + + // interpolate a quadratic function + // into the space; this function + // can be represented exactly, so + // that we can compare again later + Vector solution(dof_handler.n_dofs()); + VectorTools::interpolate (dof_handler, F(), solution); + + Functions::FEFieldFunction<2> fe_function (dof_handler, solution); + std::vector > points; + + // add only one points but also set + // the active cell to one that + // doesn't contain the current + // point. the problem happens + // because we walk over a bunch of + // cells in the process of finding + // all of these points and then + // realize when we get to the one + // at the end that the coordinates + // for this point can't be found in + // the cell we have touched last + // (it's too far away from that + // cell, and the inverse mapping + // does not converge + points.push_back (Point(-0.27999999999999992, -0.62999999999999989)); + fe_function.set_active_cell (typename DoFHandler::active_cell_iterator + (&triangulation, + 1, + 4, + &dof_handler)); + + std::vector > m (points.size()); + fe_function.vector_value_list (points, m); + + for (unsigned int i=0; i(); + + return 0; +} + diff --git a/tests/bits/fe_field_function_05_vector/cmp/generic b/tests/bits/fe_field_function_05_vector/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/bits/fe_field_function_05_vector/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK -- 2.39.5