From 85526607dbd43de76ce1c7d392ea8ced6e5184f7 Mon Sep 17 00:00:00 2001 From: Daniel Arndt Date: Sun, 26 Apr 2015 23:19:54 +0200 Subject: [PATCH] Implement point_gradient similarly to point_value --- include/deal.II/numerics/vector_tools.h | 124 +++++++++ .../deal.II/numerics/vector_tools.templates.h | 212 +++++++++++++++ source/numerics/CMakeLists.txt | 2 + .../numerics/vector_tools_point_gradient.cc | 24 ++ .../vector_tools_point_gradient.inst.in | 80 ++++++ tests/bits/point_gradient_01.cc | 212 +++++++++++++++ tests/bits/point_gradient_01.output | 22 ++ tests/bits/point_gradient_02.cc | 214 +++++++++++++++ tests/bits/point_gradient_02.output | 22 ++ tests/bits/point_gradient_hp_01.cc | 226 ++++++++++++++++ tests/bits/point_gradient_hp_01.output | 22 ++ tests/bits/point_gradient_hp_02.cc | 244 ++++++++++++++++++ tests/bits/point_gradient_hp_02.output | 40 +++ 13 files changed, 1444 insertions(+) create mode 100644 source/numerics/vector_tools_point_gradient.cc create mode 100644 source/numerics/vector_tools_point_gradient.inst.in create mode 100644 tests/bits/point_gradient_01.cc create mode 100644 tests/bits/point_gradient_01.output create mode 100644 tests/bits/point_gradient_02.cc create mode 100644 tests/bits/point_gradient_02.output create mode 100644 tests/bits/point_gradient_hp_01.cc create mode 100644 tests/bits/point_gradient_hp_01.output create mode 100644 tests/bits/point_gradient_hp_02.cc create mode 100644 tests/bits/point_gradient_hp_02.output diff --git a/include/deal.II/numerics/vector_tools.h b/include/deal.II/numerics/vector_tools.h index 78c3c1953a..7ce3c469e6 100644 --- a/include/deal.II/numerics/vector_tools.h +++ b/include/deal.II/numerics/vector_tools.h @@ -2106,6 +2106,130 @@ namespace VectorTools const InVector &fe_function, const Point &point); + /** + * Evaluate a possibly vector-valued finite element function defined by the + * given DoFHandler and nodal vector at the given point, and return the + * (vector) gradient of this function through the last argument. + * + * This is a wrapper function using a Q1-mapping for cell boundaries to call + * the other point_difference() function. + * + * @note If the cell in which the point is found is not locally owned, an + * exception of type VectorTools::ExcPointNotAvailableHere is thrown. + */ + template + void + point_gradient (const DoFHandler &dof, + const InVector &fe_function, + const Point &point, + std::vector > &value); + + /** + * Same as above for hp. + * + * @note If the cell in which the point is found is not locally owned, an + * exception of type VectorTools::ExcPointNotAvailableHere is thrown. + */ + template + void + point_gradient (const hp::DoFHandler &dof, + const InVector &fe_function, + const Point &point, + std::vector > &value); + + /** + * Evaluate a scalar finite element function defined by the given DoFHandler + * and nodal vector at the given point, and return the gradient of this + * function. + * + * Compared with the other function of the same name, this is a wrapper + * function using a Q1-mapping for cells. + * + * @note If the cell in which the point is found is not locally owned, an + * exception of type VectorTools::ExcPointNotAvailableHere is thrown. + */ + template + Tensor<1, spacedim> + point_gradient (const DoFHandler &dof, + const InVector &fe_function, + const Point &point); + + /** + * Same as above for hp. + * + * @note If the cell in which the point is found is not locally owned, an + * exception of type VectorTools::ExcPointNotAvailableHere is thrown. + */ + template + Tensor<1, spacedim> + point_gradient (const hp::DoFHandler &dof, + const InVector &fe_function, + const Point &point); + + /** + * Evaluate a possibly vector-valued finite element function defined by the + * given DoFHandler and nodal vector at the given point, and return the + * gradients of this function through the last argument. + * + * Compared with the other function of the same name, this function uses an + * arbitrary mapping for evaluation. + * + * @note If the cell in which the point is found is not locally owned, an + * exception of type VectorTools::ExcPointNotAvailableHere is thrown. + */ + template + void + point_gradient (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Point &point, + std::vector > &value); + + /** + * Same as above for hp. + * + * @note If the cell in which the point is found is not locally owned, an + * exception of type VectorTools::ExcPointNotAvailableHere is thrown. + */ + template + void + point_gradient (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const InVector &fe_function, + const Point &point, + std::vector > &value); + + /** + * Evaluate a scalar finite element function defined by the given DoFHandler + * and nodal vector at the given point, and return the gradient of this + * function. + * + * Compared with the other function of the same name, this function uses an + * arbitrary mapping for evaluation. + * + * @note If the cell in which the point is found is not locally owned, an + * exception of type VectorTools::ExcPointNotAvailableHere is thrown. + */ + template + Tensor<1, spacedim> + point_gradient (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Point &point); + + /** + * Same as above for hp. + * + * @note If the cell in which the point is found is not locally owned, an + * exception of type VectorTools::ExcPointNotAvailableHere is thrown. + */ + template + Tensor<1, spacedim> + point_gradient (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const InVector &fe_function, + const Point &point); + //@} /** * Mean value operations diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index e576dcd2ea..08ecab5eff 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -6793,6 +6793,218 @@ namespace VectorTools + template + void + point_gradient (const DoFHandler &dof, + const InVector &fe_function, + const Point &point, + std::vector > &gradients) + { + + point_gradient (StaticMappingQ1::mapping, + dof, + fe_function, + point, + gradients); + } + + + template + void + point_gradient (const hp::DoFHandler &dof, + const InVector &fe_function, + const Point &point, + std::vector > &gradients) + { + point_gradient(hp::StaticMappingQ1::mapping_collection, + dof, + fe_function, + point, + gradients); + } + + + template + Tensor<1, spacedim> + point_gradient (const DoFHandler &dof, + const InVector &fe_function, + const Point &point) + { + return point_gradient (StaticMappingQ1::mapping, + dof, + fe_function, + point); + } + + + template + Tensor<1, spacedim> + point_gradient (const hp::DoFHandler &dof, + const InVector &fe_function, + const Point &point) + { + return point_gradient(hp::StaticMappingQ1::mapping_collection, + dof, + fe_function, + point); + } + + + template + void + point_gradient (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Point &point, + std::vector > &gradient) + { + const FiniteElement &fe = dof.get_fe(); + + Assert(gradient.size() == fe.n_components(), + ExcDimensionMismatch(gradient.size(), fe.n_components())); + + // first find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const std::pair::active_cell_iterator, Point > + cell_point + = GridTools::find_active_cell_around_point (mapping, dof, point); + + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); + Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, + ExcInternalError()); + + const Quadrature + quadrature (GeometryInfo::project_to_unit_cell(cell_point.second)); + + FEValues fe_values(mapping, fe, quadrature, update_gradients); + fe_values.reinit(cell_point.first); + + // then use this to get the gradients of + // the given fe_function at this point + std::vector > > u_gradient(1, std::vector > (fe.n_components())); + fe_values.get_function_gradients(fe_function, u_gradient); + + gradient = u_gradient[0]; + } + + + template + void + point_gradient (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const InVector &fe_function, + const Point &point, + std::vector > &gradient) + { + const hp::FECollection &fe = dof.get_fe(); + + Assert(gradient.size() == fe.n_components(), + ExcDimensionMismatch(gradient.size(), fe.n_components())); + + // first find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const std::pair::active_cell_iterator, Point > + cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); + Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, + ExcInternalError()); + + const Quadrature + quadrature (GeometryInfo::project_to_unit_cell(cell_point.second)); + hp::FEValues hp_fe_values(mapping, fe, hp::QCollection(quadrature), update_gradients); + hp_fe_values.reinit(cell_point.first); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + // then use this to get the gradients of + // the given fe_function at this point + std::vector > > u_gradient(1, std::vector > (fe.n_components())); + fe_values.get_function_gradients(fe_function, u_gradient); + + gradient = u_gradient[0]; + } + + + template + Tensor<1, spacedim> + point_gradient (const Mapping &mapping, + const DoFHandler &dof, + const InVector &fe_function, + const Point &point) + { + const FiniteElement &fe = dof.get_fe(); + + Assert(fe.n_components() == 1, + ExcMessage ("Finite element is not scalar as is necessary for this function")); + + // first find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const std::pair::active_cell_iterator, Point > + cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); + Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, + ExcInternalError()); + + const Quadrature + quadrature (GeometryInfo::project_to_unit_cell(cell_point.second)); + FEValues fe_values(mapping, fe, quadrature, update_gradients); + fe_values.reinit(cell_point.first); + + // then use this to get the gradients of + // the given fe_function at this point + std::vector > u_gradient(1); + fe_values.get_function_gradients(fe_function, u_gradient); + + return u_gradient[0]; + } + + + template + Tensor<1, spacedim> + point_gradient (const hp::MappingCollection &mapping, + const hp::DoFHandler &dof, + const InVector &fe_function, + const Point &point) + { + const hp::FECollection &fe = dof.get_fe(); + + Assert(fe.n_components() == 1, + ExcMessage ("Finite element is not scalar as is necessary for this function")); + + // first find the cell in which this point + // is, initialize a quadrature rule with + // it, and then a FEValues object + const std::pair::active_cell_iterator, Point > + cell_point = GridTools::find_active_cell_around_point (mapping, dof, point); + + AssertThrow(cell_point.first->is_locally_owned(), + ExcPointNotAvailableHere()); + Assert(GeometryInfo::distance_to_unit_cell(cell_point.second) < 1e-10, + ExcInternalError()); + + const Quadrature + quadrature (GeometryInfo::project_to_unit_cell(cell_point.second)); + hp::FEValues hp_fe_values(mapping, fe, hp::QCollection(quadrature), update_gradients); + hp_fe_values.reinit(cell_point.first); + const FEValues &fe_values = hp_fe_values.get_present_fe_values(); + + // then use this to get the gradients of + // the given fe_function at this point + std::vector > u_gradient(1); + fe_values.get_function_gradients(fe_function, u_gradient); + + return u_gradient[0]; + } + + + template void subtract_mean_value(VECTOR &v, diff --git a/source/numerics/CMakeLists.txt b/source/numerics/CMakeLists.txt index c1f0cc2868..04e28cdd9f 100644 --- a/source/numerics/CMakeLists.txt +++ b/source/numerics/CMakeLists.txt @@ -41,6 +41,7 @@ SET(_src vector_tools_interpolate.cc vector_tools_mean_value.cc vector_tools_point_value.cc + vector_tools_point_gradient.cc vector_tools_project.cc vector_tools_rhs.cc ) @@ -67,6 +68,7 @@ SET(_inst vector_tools_interpolate.inst.in vector_tools_mean_value.inst.in vector_tools_point_value.inst.in + vector_tools_point_gradient.inst.in vector_tools_project.inst.in vector_tools_rhs.inst.in ) diff --git a/source/numerics/vector_tools_point_gradient.cc b/source/numerics/vector_tools_point_gradient.cc new file mode 100644 index 0000000000..dc6867f701 --- /dev/null +++ b/source/numerics/vector_tools_point_gradient.cc @@ -0,0 +1,24 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +#include + +DEAL_II_NAMESPACE_OPEN + +// ---------------------------- explicit instantiations -------------------- +#include "vector_tools_point_gradient.inst" + +DEAL_II_NAMESPACE_CLOSE diff --git a/source/numerics/vector_tools_point_gradient.inst.in b/source/numerics/vector_tools_point_gradient.inst.in new file mode 100644 index 0000000000..0844b7ce62 --- /dev/null +++ b/source/numerics/vector_tools_point_gradient.inst.in @@ -0,0 +1,80 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 1998 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +for (VEC : SERIAL_VECTORS ; deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) + { +#if deal_II_dimension == deal_II_space_dimension + + namespace VectorTools \{ + + template + void point_gradient ( + const hp::DoFHandler&, + const VEC&, + const Point&, + std::vector >&); + + template + Tensor<1,deal_II_space_dimension> point_gradient ( + const hp::DoFHandler&, + const VEC&, + const Point&); + + template + void point_gradient ( + const hp::MappingCollection&, + const hp::DoFHandler&, + const VEC&, + const Point&, + std::vector >&); + + template + Tensor<1,deal_II_space_dimension> point_gradient ( + const hp::MappingCollection&, + const hp::DoFHandler&, + const VEC&, + const Point&); + + template + void point_gradient ( + const DoFHandler&, + const VEC&, + const Point&, + std::vector >&); + + template + Tensor<1,deal_II_space_dimension> point_gradient ( + const DoFHandler&, + const VEC&, + const Point&); + + template + void point_gradient ( + const Mapping&, + const DoFHandler&, + const VEC&, + const Point&, + std::vector >&); + + template + Tensor<1,deal_II_space_dimension> point_gradient ( + const Mapping&, + const DoFHandler&, + const VEC&, + const Point&); + + \} +#endif + } diff --git a/tests/bits/point_gradient_01.cc b/tests/bits/point_gradient_01.cc new file mode 100644 index 0000000000..70238446ae --- /dev/null +++ b/tests/bits/point_gradient_01.cc @@ -0,0 +1,212 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check VectorTools::point_gradient + + + +#include "../tests.h" +#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); + } + + + virtual Tensor<1, dim> gradient (const Point &p, + const unsigned int component) const + { + Tensor<1, dim> return_value; + for (unsigned int i=0; i &p, + std::vector > &gradients) const + { + gradients[0] = gradient(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); + } + + virtual double gradient (const Point &p, + const unsigned int component) const + { + return std::exp (p(0)); + } + + virtual void vector_gradient (const Point &p, + std::vector > &gradients) const + { + gradients(0) = gradient(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); + + FE_Q element(3); + DoFHandler dof(tria); + dof.distribute_dofs(element); + + // 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.n_dofs()); + VectorTools::interpolate (dof, function, v); + + // for the following points, check the + // function gradient, output it, and + // verify that the gradient 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 > gradient(1); + for (unsigned int i=0; i<3; ++i) + { + VectorTools::point_gradient (dof, v, p[i], gradient); + deallog << -gradient[0] << std::endl; + + Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4, + ExcInternalError()); + + const Tensor<1, dim> scalar_gradient = VectorTools::point_gradient (dof, v, p[i]); + Assert (std::abs((gradient[0] - scalar_gradient).norm_square()) < 1e-4, + ExcInternalError()); + } + } + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile ("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_gradient_01.output b/tests/bits/point_gradient_01.output new file mode 100644 index 0000000000..da8670515e --- /dev/null +++ b/tests/bits/point_gradient_01.output @@ -0,0 +1,22 @@ + +DEAL:1d::0.000 +DEAL:1d::-1.000 +DEAL:1d::-0.6667 +DEAL:1d::2.150e-07 +DEAL:1d::1.111 +DEAL:1d::0.7854 +DEAL:1d::OK +DEAL:2d::0.000 0.000 +DEAL:2d::-1.000 -1.000 +DEAL:2d::-0.6667 -0.6667 +DEAL:2d::-0.0008683 0.0008683 +DEAL:2d::0.7855 0.7855 +DEAL:2d::0.6802 0.6802 +DEAL:2d::OK +DEAL:3d::2.463e-16 -2.463e-16 0.000 +DEAL:3d::-1.000 -1.000 -1.000 +DEAL:3d::-0.6667 -0.6667 -0.6667 +DEAL:3d::-0.0008683 0.0008683 0.0008683 +DEAL:3d::0.5585 0.5585 0.5585 +DEAL:3d::0.5874 0.5874 0.5874 +DEAL:3d::OK diff --git a/tests/bits/point_gradient_02.cc b/tests/bits/point_gradient_02.cc new file mode 100644 index 0000000000..0a0e6540b9 --- /dev/null +++ b/tests/bits/point_gradient_02.cc @@ -0,0 +1,214 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check VectorTools::point_gradient, alternative algorithm with mapping + + + +#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); + } + + + virtual Tensor<1, dim> gradient (const Point &p, + const unsigned int component) const + { + Tensor<1, dim> return_value; + for (unsigned int i=0; i &p, + std::vector > &gradients) const + { + gradients[0] = gradient(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); + } + + virtual double gradient (const Point &p, + const unsigned int component) const + { + return std::exp (p(0)); + } + + virtual void vector_gradient (const Point &p, + std::vector > &gradients) const + { + gradients(0) = gradient(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); + + FE_Q element(3); + DoFHandler dof(tria); + MappingQ1 mapping; + dof.distribute_dofs(element); + + // 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.n_dofs()); + VectorTools::interpolate (dof, function, v); + + // for the following points, check the + // function gradient, output it, and + // verify that the gradient 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 > gradient(1); + for (unsigned int i=0; i<3; ++i) + { + VectorTools::point_gradient (mapping, dof, v, p[i], gradient); + deallog << -gradient[0] << std::endl; + + Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4, + ExcInternalError()); + + const Tensor<1, dim> scalar_gradient = VectorTools::point_gradient (dof, v, p[i]); + Assert (std::abs((gradient[0] - scalar_gradient).norm_square()) < 1e-4, + ExcInternalError()); + } + } + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile ("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_gradient_02.output b/tests/bits/point_gradient_02.output new file mode 100644 index 0000000000..da8670515e --- /dev/null +++ b/tests/bits/point_gradient_02.output @@ -0,0 +1,22 @@ + +DEAL:1d::0.000 +DEAL:1d::-1.000 +DEAL:1d::-0.6667 +DEAL:1d::2.150e-07 +DEAL:1d::1.111 +DEAL:1d::0.7854 +DEAL:1d::OK +DEAL:2d::0.000 0.000 +DEAL:2d::-1.000 -1.000 +DEAL:2d::-0.6667 -0.6667 +DEAL:2d::-0.0008683 0.0008683 +DEAL:2d::0.7855 0.7855 +DEAL:2d::0.6802 0.6802 +DEAL:2d::OK +DEAL:3d::2.463e-16 -2.463e-16 0.000 +DEAL:3d::-1.000 -1.000 -1.000 +DEAL:3d::-0.6667 -0.6667 -0.6667 +DEAL:3d::-0.0008683 0.0008683 0.0008683 +DEAL:3d::0.5585 0.5585 0.5585 +DEAL:3d::0.5874 0.5874 0.5874 +DEAL:3d::OK diff --git a/tests/bits/point_gradient_hp_01.cc b/tests/bits/point_gradient_hp_01.cc new file mode 100644 index 0000000000..4319a3e52f --- /dev/null +++ b/tests/bits/point_gradient_hp_01.cc @@ -0,0 +1,226 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check VectorTools::point_gradient 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); + } + + + virtual Tensor<1, dim> gradient (const Point &p, + const unsigned int component) const + { + Tensor<1, dim> return_value; + for (unsigned int i=0; i &p, + std::vector > &gradients) const + { + gradients[0] = gradient(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); + } + + virtual double gradient (const Point &p, + const unsigned int component) const + { + return std::exp (p(0)); + } + + virtual void vector_gradient (const Point &p, + std::vector > &gradients) const + { + gradients(0) = gradient(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 (Testing::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 gradient, output it, and + // verify that the gradient 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 > gradient(1); + for (unsigned int i=0; i<3; ++i) + { + VectorTools::point_gradient (dof_handler, v, p[i], gradient); + deallog << -gradient[0] << std::endl; + + Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4, + ExcInternalError()); + + const Tensor<1, dim> scalar_gradient = VectorTools::point_gradient (dof_handler, v, p[i]); + Assert (std::abs((gradient[0] - scalar_gradient).norm_square()) < 1e-4, + ExcInternalError()); + } + } + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile ("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_gradient_hp_01.output b/tests/bits/point_gradient_hp_01.output new file mode 100644 index 0000000000..d6adfac86a --- /dev/null +++ b/tests/bits/point_gradient_hp_01.output @@ -0,0 +1,22 @@ + +DEAL:1d::0.000 +DEAL:1d::-1.000 +DEAL:1d::-0.6667 +DEAL:1d::2.150e-07 +DEAL:1d::1.111 +DEAL:1d::0.7854 +DEAL:1d::OK +DEAL:2d::-3.553e-15 0.000 +DEAL:2d::-1.000 -1.000 +DEAL:2d::-0.6667 -0.6667 +DEAL:2d::-4.747e-06 4.747e-06 +DEAL:2d::0.7854 0.7854 +DEAL:2d::0.6802 0.6802 +DEAL:2d::OK +DEAL:3d::-3.031e-15 -5.216e-16 0.000 +DEAL:3d::-1.000 -1.000 -1.000 +DEAL:3d::-0.6667 -0.6667 -0.6667 +DEAL:3d::-4.747e-06 4.747e-06 4.747e-06 +DEAL:3d::0.5553 0.5553 0.5553 +DEAL:3d::0.5890 0.5890 0.5890 +DEAL:3d::OK diff --git a/tests/bits/point_gradient_hp_02.cc b/tests/bits/point_gradient_hp_02.cc new file mode 100644 index 0000000000..4bcbe1c3e5 --- /dev/null +++ b/tests/bits/point_gradient_hp_02.cc @@ -0,0 +1,244 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2004 - 2014 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// check VectorTools::point_gradient 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); + } + + + virtual Tensor<1, dim> gradient (const Point &p, + const unsigned int component) const + { + Tensor<1, dim> return_value; + for (unsigned int i=0; i &p, + std::vector > &gradients) const + { + gradients[0] = gradient(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); + } + + virtual double gradient (const Point &p, + const unsigned int component) const + { + return std::exp (p(0)); + } + + virtual void vector_gradient (const Point &p, + std::vector > &gradients) const + { + gradients(0) = gradient(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 (Testing::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 gradient, output it, and + // verify that the gradient 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 > gradient(1); + for (unsigned int i=0; i<3; ++i) + { + VectorTools::point_gradient (mapping_1, dof_handler, v, p[i], gradient); + deallog << -gradient[0] << std::endl; + + Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4, + ExcInternalError()); + VectorTools::point_gradient (mapping_2, dof_handler, v, p[i], gradient); + deallog << -gradient[0] << std::endl; + + Assert (std::abs((gradient[0] - function.gradient(p[i])).norm_square()) < 1e-4, + ExcInternalError()); + + const Tensor<1, dim> scalar_gradient_1 = VectorTools::point_gradient (mapping_1, dof_handler, v, p[i]); + Assert (std::abs((gradient[0] - scalar_gradient_1).norm_square()) < 1e-4, + ExcInternalError()); + + const Tensor<1, dim> scalar_gradient_2 = VectorTools::point_gradient (mapping_2, dof_handler, v, p[i]); + Assert (std::abs((gradient[0] - scalar_gradient_2).norm_square()) < 1e-4, + ExcInternalError()); + } + } + + deallog << "OK" << std::endl; +} + + +int main () +{ + std::ofstream logfile ("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_gradient_hp_02.output b/tests/bits/point_gradient_hp_02.output new file mode 100644 index 0000000000..600d18b6de --- /dev/null +++ b/tests/bits/point_gradient_hp_02.output @@ -0,0 +1,40 @@ + +DEAL:1d::0.000 +DEAL:1d::0.000 +DEAL:1d::-1.000 +DEAL:1d::-1.000 +DEAL:1d::-0.6667 +DEAL:1d::-0.6667 +DEAL:1d::2.150e-07 +DEAL:1d::2.150e-07 +DEAL:1d::1.111 +DEAL:1d::1.111 +DEAL:1d::0.7854 +DEAL:1d::0.7854 +DEAL:1d::OK +DEAL:2d::-3.553e-15 0.000 +DEAL:2d::-3.553e-15 0.000 +DEAL:2d::-1.000 -1.000 +DEAL:2d::-1.000 -1.000 +DEAL:2d::-0.6667 -0.6667 +DEAL:2d::-0.6667 -0.6667 +DEAL:2d::-4.747e-06 4.747e-06 +DEAL:2d::-4.747e-06 4.747e-06 +DEAL:2d::0.7854 0.7854 +DEAL:2d::0.7854 0.7854 +DEAL:2d::0.6802 0.6802 +DEAL:2d::0.6802 0.6802 +DEAL:2d::OK +DEAL:3d::-3.031e-15 -5.216e-16 0.000 +DEAL:3d::-3.031e-15 -5.216e-16 0.000 +DEAL:3d::-1.000 -1.000 -1.000 +DEAL:3d::-1.000 -1.000 -1.000 +DEAL:3d::-0.6667 -0.6667 -0.6667 +DEAL:3d::-0.6667 -0.6667 -0.6667 +DEAL:3d::-4.747e-06 4.747e-06 4.747e-06 +DEAL:3d::-4.747e-06 4.747e-06 4.747e-06 +DEAL:3d::0.5553 0.5553 0.5553 +DEAL:3d::0.5553 0.5553 0.5553 +DEAL:3d::0.5890 0.5890 0.5890 +DEAL:3d::0.5890 0.5890 0.5890 +DEAL:3d::OK -- 2.39.5