From: Simon Sticko Date: Wed, 8 Jan 2020 13:37:47 +0000 (+0100) Subject: Implement jump_hessian and jump_3rd_derivative on FEInterfaceValues. X-Git-Tag: v9.2.0-rc1~683^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=dac138ac42243a0fc476f9193b4a6929b361d95b;p=dealii.git Implement jump_hessian and jump_3rd_derivative on FEInterfaceValues. In the same way as the other jump_* functions. --- diff --git a/include/deal.II/fe/fe_interface_values.h b/include/deal.II/fe/fe_interface_values.h index 4b3b173425..b6923b0ddc 100644 --- a/include/deal.II/fe/fe_interface_values.h +++ b/include/deal.II/fe/fe_interface_values.h @@ -371,6 +371,34 @@ public: const unsigned int q_point, const unsigned int component = 0) const; + /** + * Return the jump in the Hessian $[\nabla^2 u] = \nabla^2 u_{\text{cell0}} - + * \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> + jump_hessian(const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component = 0) const; + + /** + * Return the jump in the third derivative $[\nabla^3 u] = \nabla^3 + * u_{\text{cell0}} - \nabla^3 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^3 u] = \nabla^3 u_{\text{cell0}}$. + */ + Tensor<3, dim> + jump_3rd_derivative(const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component = 0) const; + /** * @} */ @@ -872,6 +900,67 @@ FEInterfaceValues::jump_gradient( } + +template +Tensor<2, dim> +FEInterfaceValues::jump_hessian( + const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component) const +{ + const auto dof_pair = dofmap[interface_dof_index]; + + if (at_boundary()) + return get_fe_face_values(0).shape_hessian_component(dof_pair[0], + q_point, + component); + + Tensor<2, dim> value; + + if (dof_pair[0] != numbers::invalid_unsigned_int) + value += 1.0 * get_fe_face_values(0).shape_hessian_component(dof_pair[0], + q_point, + component); + if (dof_pair[1] != numbers::invalid_unsigned_int) + value += -1.0 * get_fe_face_values(1).shape_hessian_component(dof_pair[1], + q_point, + component); + + return value; +} + + +template +Tensor<3, dim> +FEInterfaceValues::jump_3rd_derivative( + const unsigned int interface_dof_index, + const unsigned int q_point, + const unsigned int component) const +{ + const auto dof_pair = dofmap[interface_dof_index]; + + if (at_boundary()) + return get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0], + q_point, + component); + + Tensor<3, dim> value; + + if (dof_pair[0] != numbers::invalid_unsigned_int) + value += + 1.0 * get_fe_face_values(0).shape_3rd_derivative_component(dof_pair[0], + q_point, + component); + if (dof_pair[1] != numbers::invalid_unsigned_int) + value += + -1.0 * get_fe_face_values(1).shape_3rd_derivative_component(dof_pair[1], + q_point, + component); + + return value; +} + + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/tests/feinterface/fe_interface_values_09.cc b/tests/feinterface/fe_interface_values_09.cc new file mode 100644 index 0000000000..55766fbcf4 --- /dev/null +++ b/tests/feinterface/fe_interface_values_09.cc @@ -0,0 +1,150 @@ +// --------------------------------------------------------------------- +// +// 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 jump_hessian and jump_third_derivative 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 +print_norm_of_average_over_quadrature_points(const FEInterfaceValues &fiv) +{ + const unsigned int n_dofs = fiv.n_current_interface_dofs(); + Vector cell_vector(n_dofs); + + cell_vector = 0.0; + for (unsigned int qpoint = 0; qpoint < fiv.n_quadrature_points; ++qpoint) + for (unsigned int i = 0; i < n_dofs; ++i) + cell_vector(i) += + fiv.jump_hessian(i, qpoint).norm() * fiv.get_JxW_values()[qpoint]; + deallog << "jump_hessian.norm(): " << cell_vector; + + cell_vector = 0.0; + for (unsigned int qpoint = 0; qpoint < fiv.n_quadrature_points; ++qpoint) + for (unsigned int i = 0; i < n_dofs; ++i) + cell_vector(i) += fiv.jump_3rd_derivative(i, qpoint).norm() * + fiv.get_JxW_values()[qpoint]; + deallog << "jump_3rd_derivative.norm(): " << cell_vector; +} + + + +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); + const UpdateFlags update_flags = + update_hessians | update_3rd_derivatives | update_JxW_values; + + FEInterfaceValues fiv(mapping, + fe, + QGauss(fe.degree + 1), + update_flags); + + auto cell = dofh.begin(); + + // Print jump in Hessian and third derivative for a face at the boundary. + { + unsigned int face = 0; + Assert(cell->at_boundary(face), ExcInternalError()); + + fiv.reinit(cell, face); + + print_norm_of_average_over_quadrature_points(fiv); + } + + // Print jump in Hessian and third derivative for a face between two cells. + { + const unsigned int face = 1; + Assert(!cell->at_boundary(face), ExcInternalError()); + + fiv.reinit(cell, + face, + numbers::invalid_unsigned_int, + cell->neighbor(face), + cell->neighbor_of_neighbor(face), + numbers::invalid_unsigned_int); + + print_norm_of_average_over_quadrature_points(fiv); + } +} + + + +int +main() +{ + initlog(); + // Test the lowest order cg and dg elements which have a non-zero third + // derivative. + test<2>(FE_Q<2>(2)); + test<2>(FE_DGQ<2>(2)); + test<3>(FE_Q<3>(1)); + test<3>(FE_DGQ<3>(1)); +} diff --git a/tests/feinterface/fe_interface_values_09.output b/tests/feinterface/fe_interface_values_09.output new file mode 100644 index 0000000000..904fbfa46f --- /dev/null +++ b/tests/feinterface/fe_interface_values_09.output @@ -0,0 +1,21 @@ + +DEAL::FE_Q<2>(2) +DEAL::jump_hessian.norm(): 7.17398 2.12446 7.17398 2.12446 12.5704 4.36931 7.68564 7.68564 13.4538 +DEAL::jump_3rd_derivative.norm(): 23.1831 11.8202 23.1831 11.8202 44.4667 20.3528 34.2248 34.2248 63.5828 +DEAL::jump_hessian.norm(): 2.12446 11.0742 2.12446 11.0742 4.36931 14.6059 7.68564 7.68564 13.4538 2.12446 2.12446 4.36931 7.68564 7.68564 13.4538 +DEAL::jump_3rd_derivative.norm(): 11.8202 41.5692 11.8202 41.5692 20.3528 83.1384 34.2248 34.2248 63.5828 11.8202 11.8202 20.3528 34.2248 34.2248 63.5828 +DEAL::FE_DGQ<2>(2) +DEAL::jump_hessian.norm(): 7.17398 7.68564 2.12446 12.5704 13.4538 4.36931 7.17398 7.68564 2.12446 +DEAL::jump_3rd_derivative.norm(): 23.1831 34.2248 11.8202 44.4667 63.5828 20.3528 23.1831 34.2248 11.8202 +DEAL::jump_hessian.norm(): 2.12446 7.68564 7.17398 4.36931 13.4538 12.5704 2.12446 7.68564 7.17398 7.17398 7.68564 2.12446 12.5704 13.4538 4.36931 7.17398 7.68564 2.12446 +DEAL::jump_3rd_derivative.norm(): 11.8202 34.2248 23.1831 20.3528 63.5828 44.4667 11.8202 34.2248 23.1831 23.1831 34.2248 11.8202 44.4667 63.5828 20.3528 23.1831 34.2248 11.8202 +DEAL::FE_Q<3>(1) +DEAL::jump_hessian.norm(): 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 +DEAL::jump_3rd_derivative.norm(): 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 +DEAL::jump_hessian.norm(): 1.07735 2.15470 1.07735 2.15470 1.07735 2.15470 1.07735 2.15470 1.07735 1.07735 1.07735 1.07735 +DEAL::jump_3rd_derivative.norm(): 2.44949 4.89898 2.44949 4.89898 2.44949 4.89898 2.44949 4.89898 2.44949 2.44949 2.44949 2.44949 +DEAL::FE_DGQ<3>(1) +DEAL::jump_hessian.norm(): 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 +DEAL::jump_3rd_derivative.norm(): 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 +DEAL::jump_hessian.norm(): 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 1.81150 1.07735 +DEAL::jump_3rd_derivative.norm(): 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949 2.44949