From bb0735cd65243d2f41fffcd556a3a311a3c6ed8c Mon Sep 17 00:00:00 2001 From: Timo Heister Date: Wed, 2 Oct 2019 13:20:15 -0400 Subject: [PATCH] FEInterfaceValues: gradients and hessians add average_gradient(), average_hessian(), jump_gradient() --- include/deal.II/fe/fe_interface_values.h | 130 ++++++++++++++ tests/feinterface/fe_interface_values_08.cc | 170 ++++++++++++++++++ .../feinterface/fe_interface_values_08.output | 50 ++++++ 3 files changed, 350 insertions(+) create mode 100644 tests/feinterface/fe_interface_values_08.cc create mode 100644 tests/feinterface/fe_interface_values_08.output diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h index 9e611039cc..50100a27c8 100644 --- a/include/deal.II/fe/fe_interface_values.h +++ b/include/deal.II/fe/fe_interface_values.h @@ -307,6 +307,47 @@ public: const unsigned int q_point, const unsigned int component = 0) const; + /** + * Return the average of the gradient $\{\nabla u \} = \frac{1}{2}\nabla + * u_{\text{cell0}} + \frac{1}{2} \nabla u_{\text{cell1}}$ on the interface + * for the shape function @p interface_dof_index at the quadrature point @p + * q_point of component @p component. + * + * If this is a boundary face (at_boundary() returns true), then + * $\{\nabla u\}=\nabla u_{\text{cell0}}$. + */ + Tensor<1, dim> + average_gradient(const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component = 0) const; + + /** + * Return the average of the hessian $\{\nabla^2 u \} = \frac{1}{2}\nabla^2 + * u_{\text{cell0}} + \frac{1}{2} \nabla^2 u_{\text{cell1}}$ on the interface + * for the shape function @p interface_dof_index at the quadrature point @p + * q_point of component @p component. + * + * If this is a boundary face (at_boundary() returns true), then + * $\{\nabla^2 u\}=\nabla^2 u_{\text{cell0}}$. + */ + Tensor<2, dim> + average_hessian(const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component = 0) const; + + /** + * Return the jump if the gradient $[\nabla u]=\nabla u_{\text{cell0}} - + * \nabla u_{\text{cell1}}$ on the interface for the shape function @p + * interface_dof_index at the quadrature point @p q_point of component @p + * component. + * + * If this is a boundary face (at_boundary() returns true), then + * $[\nabla u]=\nabla u_{\text{cell0}}$. + */ + Tensor<1, dim> + jump_gradient(const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component = 0) const; /** * @} @@ -705,6 +746,95 @@ FEInterfaceValues::average( return value; } + + +template +Tensor<1, dim> +FEInterfaceValues::average_gradient( + const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component) const +{ + const auto dofpair = dofmap[interface_dof_index]; + + if (at_boundary()) + return get_fe_face_values(0).shape_grad_component(dofpair[0], + q_point, + component); + + Tensor<1, dim> value; + + if (dofpair[0] != numbers::invalid_unsigned_int) + value += 0.5 * get_fe_face_values(0).shape_grad_component(dofpair[0], + q_point, + component); + if (dofpair[1] != numbers::invalid_unsigned_int) + value += 0.5 * get_fe_face_values(1).shape_grad_component(dofpair[1], + q_point, + component); + + return value; +} + + + +template +Tensor<2, dim> +FEInterfaceValues::average_hessian( + const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component) const +{ + const auto dofpair = dofmap[interface_dof_index]; + + if (at_boundary()) + return get_fe_face_values(0).shape_hessian_component(dofpair[0], + q_point, + component); + + Tensor<2, dim> value; + + if (dofpair[0] != numbers::invalid_unsigned_int) + value += 0.5 * get_fe_face_values(0).shape_hessian_component(dofpair[0], + q_point, + component); + if (dofpair[1] != numbers::invalid_unsigned_int) + value += 0.5 * get_fe_face_values(1).shape_hessian_component(dofpair[1], + q_point, + component); + + return value; +} + +template +Tensor<1, dim> +FEInterfaceValues::jump_gradient( + const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component) const +{ + const auto dofpair = dofmap[interface_dof_index]; + + if (at_boundary()) + return get_fe_face_values(0).shape_grad_component(dofpair[0], + q_point, + component); + + Tensor<1, dim> value; + + if (dofpair[0] != numbers::invalid_unsigned_int) + value += 1.0 * get_fe_face_values(0).shape_grad_component(dofpair[0], + q_point, + component); + if (dofpair[1] != numbers::invalid_unsigned_int) + value += -1.0 * get_fe_face_values(1).shape_grad_component(dofpair[1], + q_point, + component); + + return value; +} + + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/tests/feinterface/fe_interface_values_08.cc b/tests/feinterface/fe_interface_values_08.cc new file mode 100644 index 0000000000..8568b49c23 --- /dev/null +++ b/tests/feinterface/fe_interface_values_08.cc @@ -0,0 +1,170 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2018 - 2019 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.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// evaluate average_gradient average_hessian jump_gradient of FEInterfaceValues + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +template +void +make_2_cells(Triangulation &tria); + +template <> +void make_2_cells<2>(Triangulation<2> &tria) +{ + const unsigned int dim = 2; + std::vector repetitions = {2, 1}; + Point p1; + Point p2(2.0, 1.0); + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); +} + +template <> +void make_2_cells<3>(Triangulation<3> &tria) +{ + const unsigned int dim = 3; + std::vector repetitions = {2, 1, 1}; + Point p1; + Point p2(2.0, 1.0, 1.0); + + GridGenerator::subdivided_hyper_rectangle(tria, repetitions, p1, p2); +} + + +template +void +test(const FiniteElement &fe) +{ + Triangulation tria; + make_2_cells(tria); + + DoFHandler dofh(tria); + deallog << fe.get_name() << std::endl; + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_values | update_gradients | + update_hessians | update_quadrature_points | + update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + + auto cell = dofh.begin(); + + { + unsigned int f = 0; + Assert(cell->at_boundary(f), ExcInternalError()); + + fiv.reinit(cell, f); + + const unsigned int n_dofs = fiv.n_current_interface_dofs(); + Vector cell_vector(n_dofs); + + const auto &q_points = fiv.get_quadrature_points(); + + cell_vector = 0.0; + for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) + for (unsigned int i = 0; i < n_dofs; ++i) + cell_vector(i) += + fiv.average_gradient(i, qpoint).norm() * fiv.get_JxW_values()[qpoint]; + deallog << "average_gradient.norm(): " << cell_vector; + + cell_vector = 0.0; + for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) + for (unsigned int i = 0; i < n_dofs; ++i) + cell_vector(i) += + fiv.average_hessian(i, qpoint).norm() * fiv.get_JxW_values()[qpoint]; + deallog << "average_hessian.norm(): " << cell_vector; + + cell_vector = 0.0; + for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) + for (unsigned int i = 0; i < n_dofs; ++i) + cell_vector(i) += + fiv.jump_gradient(i, qpoint).norm() * fiv.get_JxW_values()[qpoint]; + deallog << "jump_gradient.norm(): " << cell_vector; + } + + + for (unsigned int f = 0; f < GeometryInfo::faces_per_cell; ++f) + if (!cell->at_boundary(f)) + { + fiv.reinit(cell, + f, + numbers::invalid_unsigned_int, + cell->neighbor(f), + cell->neighbor_of_neighbor(f), + numbers::invalid_unsigned_int); + + const unsigned int n_dofs = fiv.n_current_interface_dofs(); + Vector cell_vector(n_dofs); + + const auto &q_points = fiv.get_quadrature_points(); + + cell_vector = 0.0; + for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) + for (unsigned int i = 0; i < n_dofs; ++i) + cell_vector(i) += fiv.average_gradient(i, qpoint).norm() * + fiv.get_JxW_values()[qpoint]; + deallog << "average_gradient.norm(): " << cell_vector; + + cell_vector = 0.0; + for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) + for (unsigned int i = 0; i < n_dofs; ++i) + cell_vector(i) += fiv.average_hessian(i, qpoint).norm() * + fiv.get_JxW_values()[qpoint]; + deallog << "average_hessian.norm(): " << cell_vector; + + cell_vector = 0.0; + for (unsigned int qpoint = 0; qpoint < q_points.size(); ++qpoint) + for (unsigned int i = 0; i < n_dofs; ++i) + cell_vector(i) += fiv.jump_gradient(i, qpoint).norm() * + fiv.get_JxW_values()[qpoint]; + deallog << "jump_gradient.norm(): " << cell_vector; + } +} + + + +int +main() +{ + initlog(); + test<2>(FE_Q<2>(1)); + test<2>(FE_DGQ<2>(0)); + test<2>(FE_DGQ<2>(1)); + test<2>(FE_DGQ<2>(2)); + test<3>(FE_Q<3>(1)); + test<3>(FE_DGQ<3>(0)); + test<3>(FE_DGQ<3>(1)); +} diff --git a/tests/feinterface/fe_interface_values_08.output b/tests/feinterface/fe_interface_values_08.output new file mode 100644 index 0000000000..6487816da7 --- /dev/null +++ b/tests/feinterface/fe_interface_values_08.output @@ -0,0 +1,50 @@ + +DEAL::FE_Q<2>(1) +DEAL::average_gradient.norm(): 1.14783 0.500000 1.14783 0.500000 +DEAL::average_hessian.norm(): 1.41421 1.41421 1.41421 1.41421 +DEAL::jump_gradient.norm(): 1.14783 0.500000 1.14783 0.500000 +DEAL::average_gradient.norm(): 0.250000 1.00000 0.250000 1.00000 0.250000 0.250000 +DEAL::average_hessian.norm(): 0.707107 0.00000 0.707107 0.00000 0.707107 0.707107 +DEAL::jump_gradient.norm(): 0.500000 1.00000 0.500000 1.00000 0.500000 0.500000 +DEAL::FE_DGQ<2>(0) +DEAL::average_gradient.norm(): 0.00000 +DEAL::average_hessian.norm(): 0.00000 +DEAL::jump_gradient.norm(): 0.00000 +DEAL::average_gradient.norm(): 0.00000 0.00000 +DEAL::average_hessian.norm(): 0.00000 0.00000 +DEAL::jump_gradient.norm(): 0.00000 0.00000 +DEAL::FE_DGQ<2>(1) +DEAL::average_gradient.norm(): 1.14783 0.500000 1.14783 0.500000 +DEAL::average_hessian.norm(): 1.41421 1.41421 1.41421 1.41421 +DEAL::jump_gradient.norm(): 1.14783 0.500000 1.14783 0.500000 +DEAL::average_gradient.norm(): 0.250000 0.573917 0.250000 0.573917 0.573917 0.250000 0.573917 0.250000 +DEAL::average_hessian.norm(): 0.707107 0.707107 0.707107 0.707107 0.707107 0.707107 0.707107 0.707107 +DEAL::jump_gradient.norm(): 0.500000 1.14783 0.500000 1.14783 1.14783 0.500000 1.14783 0.500000 +DEAL::FE_DGQ<2>(2) +DEAL::average_gradient.norm(): 1.52420 0.860663 0.215166 3.17925 2.66667 0.666667 1.52420 0.860663 0.215166 +DEAL::average_hessian.norm(): 7.17398 7.68564 2.12446 12.5704 13.4538 4.36931 7.17398 7.68564 2.12446 +DEAL::jump_gradient.norm(): 1.52420 0.860663 0.215166 3.17925 2.66667 0.666667 1.52420 0.860663 0.215166 +DEAL::average_gradient.norm(): 0.107583 0.430331 0.762102 0.333333 1.33333 1.58962 0.107583 0.430331 0.762102 0.762102 0.430331 0.107583 1.58962 1.33333 0.333333 0.762102 0.430331 0.107583 +DEAL::average_hessian.norm(): 1.06223 3.84282 3.58699 2.18466 6.72690 6.28519 1.06223 3.84282 3.58699 3.58699 3.84282 1.06223 6.28519 6.72690 2.18466 3.58699 3.84282 1.06223 +DEAL::jump_gradient.norm(): 0.215166 0.860663 1.52420 0.666667 2.66667 3.17925 0.215166 0.860663 1.52420 1.52420 0.860663 0.215166 3.17925 2.66667 0.666667 1.52420 0.860663 0.215166 +DEAL::FE_Q<3>(1) +DEAL::average_gradient.norm(): 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 +DEAL::average_hessian.norm(): 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 +DEAL::jump_gradient.norm(): 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 +DEAL::average_gradient.norm(): 0.125000 0.761802 0.125000 0.761802 0.125000 0.761802 0.125000 0.761802 0.125000 0.125000 0.125000 0.125000 +DEAL::average_hessian.norm(): 0.538675 1.41421 0.538675 1.41421 0.538675 1.41421 0.538675 1.41421 0.538675 0.538675 0.538675 0.538675 +DEAL::jump_gradient.norm(): 0.250000 0.500000 0.250000 0.500000 0.250000 0.500000 0.250000 0.500000 0.250000 0.250000 0.250000 0.250000 +DEAL::FE_DGQ<3>(0) +DEAL::average_gradient.norm(): 0.00000 +DEAL::average_hessian.norm(): 0.00000 +DEAL::jump_gradient.norm(): 0.00000 +DEAL::average_gradient.norm(): 0.00000 0.00000 +DEAL::average_hessian.norm(): 0.00000 0.00000 +DEAL::jump_gradient.norm(): 0.00000 0.00000 +DEAL::FE_DGQ<3>(1) +DEAL::average_gradient.norm(): 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 +DEAL::average_hessian.norm(): 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 +DEAL::jump_gradient.norm(): 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 +DEAL::average_gradient.norm(): 0.125000 0.405739 0.125000 0.405739 0.125000 0.405739 0.125000 0.405739 0.405739 0.125000 0.405739 0.125000 0.405739 0.125000 0.405739 0.125000 +DEAL::average_hessian.norm(): 0.538675 0.905750 0.538675 0.905750 0.538675 0.905750 0.538675 0.905750 0.905750 0.538675 0.905750 0.538675 0.905750 0.538675 0.905750 0.538675 +DEAL::jump_gradient.norm(): 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 0.811479 0.250000 -- 2.39.5