From: David Wells Date: Fri, 5 Jul 2024 18:09:44 +0000 (-0400) Subject: Add FETools::compute_nodal_quadrature(). X-Git-Tag: v9.6.0-rc1~126^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F17214%2Fhead;p=dealii.git Add FETools::compute_nodal_quadrature(). --- diff --git a/doc/news/changes/minor/20240705DavidWells b/doc/news/changes/minor/20240705DavidWells new file mode 100644 index 0000000000..043a06dc66 --- /dev/null +++ b/doc/news/changes/minor/20240705DavidWells @@ -0,0 +1,5 @@ +New: Added a function FETools::compute_nodal_quadrature() +which computes the nodal Quadrature rule for an interpolatory +FiniteElement. +
+(David Wells, 2023/07/05) diff --git a/include/deal.II/fe/fe_tools.h b/include/deal.II/fe/fe_tools.h index bdddae6fe7..6897a9c1d5 100644 --- a/include/deal.II/fe/fe_tools.h +++ b/include/deal.II/fe/fe_tools.h @@ -586,7 +586,41 @@ namespace FETools const unsigned int face, FullMatrix &X); - + /** + * @brief Compute the nodal quadrature rule associated with an element. + * + * Many common finite elements (e.g., FE_Q) are interpolatory - i.e., they are + * multidimensional generalizations of Lagrange interpolating polynomials and + * are collectively referred to as @ref GlossLagrange "Lagrange elements". + * + * Since these elements define a polynomial interpolation, they implicitly + * also define a Quadrature rule. For example, the weights of QTrapezoid are + * the integrals of multilinear shape functions and the quadrature points are + * the vertices of the reference hypercube - i.e., QTrapezoid is the nodal + * rule corresponding to FE_Q(1). Some elements, such as FE_SimplexP(2), do + * not correspond to useful nodal quadrature rules since the integrals of some + * shape functions (in this case, the ones with support points at vertices) + * are zero. Other elements, like FE_RaviartThomas, are noninterpolatory and + * therefore also do not correspond to nodal quadrature rules. + * + * Given a FiniteElement with support points, this function determines the + * associated quadrature rule by setting the weights equal to the integrals of + * its shape functions (even if they are zero or negative) and points equal to + * the support points of the FiniteElement. + * + * @note If @p fe is an FE_Q or FE_DGQ with the default support points (i.e., + * it was constructed with a polynomial order argument and not with a custom + * list of support points) then the corresponding nodal quadrature rule has + * the same points and weights as QGaussLobatto (but they will generally be in + * a different order). + * + * @note This function is only implemented for scalar (i.e., one block and one + * component) elements which have + * @ref GlossSupport "support points". + */ + template + Quadrature + compute_nodal_quadrature(const FiniteElement &fe); /** * Wrapper around diff --git a/include/deal.II/fe/fe_tools.templates.h b/include/deal.II/fe/fe_tools.templates.h index f499d71169..e7636e5370 100644 --- a/include/deal.II/fe/fe_tools.templates.h +++ b/include/deal.II/fe/fe_tools.templates.h @@ -2947,6 +2947,46 @@ namespace FETools + template + Quadrature + compute_nodal_quadrature(const FiniteElement &fe) + { + Assert(fe.has_support_points(), ExcNotImplemented()); + Assert(fe.n_blocks() == 1, ExcNotImplemented()); + Assert(fe.n_components() == 1, ExcNotImplemented()); + + const ReferenceCell type = fe.reference_cell(); + + const Quadrature q_gauss = + type.get_gauss_type_quadrature(fe.tensor_degree() + 1); + Triangulation tria; + GridGenerator::reference_cell(tria, type); + const Mapping &mapping = + type.template get_default_linear_mapping(); + + auto cell = tria.begin_active(); + FEValues fe_values(mapping, + fe, + q_gauss, + update_values | update_JxW_values); + fe_values.reinit(cell); + + const std::vector> &nodal_quad_points = + fe.get_unit_support_points(); + std::vector nodal_quad_weights(nodal_quad_points.size()); + Assert(nodal_quad_points.size() > 0, ExcNotImplemented()); + for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) + { + double integral = 0.0; + for (unsigned int q = 0; q < q_gauss.size(); ++q) + integral += fe_values.shape_value(i, q) * fe_values.JxW(q); + nodal_quad_weights[i] = integral; + } + return {nodal_quad_points, nodal_quad_weights}; + } + + + namespace internal { namespace FEToolsConvertHelper diff --git a/source/fe/fe_tools.inst.in b/source/fe/fe_tools.inst.in index 7071c7e7cd..a2bff7ba44 100644 --- a/source/fe/fe_tools.inst.in +++ b/source/fe/fe_tools.inst.in @@ -170,6 +170,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS) const Quadrature &, const Quadrature &, FullMatrix &); + + template Quadrature + compute_nodal_quadrature( + const FiniteElement &); #endif #if deal_II_dimension == deal_II_space_dimension diff --git a/tests/fe/fe_tools_nodal_quadrature.cc b/tests/fe/fe_tools_nodal_quadrature.cc new file mode 100644 index 0000000000..2555a74dbf --- /dev/null +++ b/tests/fe/fe_tools_nodal_quadrature.cc @@ -0,0 +1,94 @@ +// ------------------------------------------------------------------------ +// +// SPDX-License-Identifier: LGPL-2.1-or-later +// Copyright (C) 2024 by the deal.II authors +// +// This file is part of the deal.II library. +// +// Part of the source code is dual licensed under Apache-2.0 WITH +// LLVM-exception OR LGPL-2.1-or-later. Detailed license information +// governing the source code and code contributions can be found in +// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II. +// +// ------------------------------------------------------------------------ + + +// Test FETools::compute_nodal_quadrature() + +#include + +#include +#include +#include +#include +#include + +#include "../tests.h" + +template +void +print_quadrature(const FE_Q_Base &fe, std::ostream &out) +{ + out << fe.get_name() << std::endl; + const Quadrature nodal_quadrature = + FETools::compute_nodal_quadrature(fe); + + out << "points:" << std::endl; + for (const auto &point : nodal_quadrature.get_points()) + out << " " << point << std::endl; + out << "weights:" << std::endl; + for (const auto &weight : nodal_quadrature.get_weights()) + out << " " << weight << std::endl; + + // Check that we are close enough to QGaussLobatto + const Quadrature lobatto = QGaussLobatto(fe.tensor_degree() + 1); + + for (unsigned int q = 0; q < lobatto.size(); ++q) + { + AssertThrow((lobatto.point(q) - nodal_quadrature.point(q)).norm() < 1e-6, + ExcInternalError()); + AssertThrow(std::abs(lobatto.weight(q) - nodal_quadrature.weight(q)) < + 1e-6, + ExcInternalError()); + } +} + +template +void +print_quadrature(const FiniteElement &fe, std::ostream &out) +{ + out << fe.get_name() << std::endl; + const Quadrature nodal_quadrature = + FETools::compute_nodal_quadrature(fe); + + out << "points:" << std::endl; + for (const auto &point : nodal_quadrature.get_points()) + out << " " << point << std::endl; + out << "weights:" << std::endl; + for (const auto &weight : nodal_quadrature.get_weights()) + out << " " << weight << std::endl; +} + +int +main() +{ + initlog(); + + print_quadrature(FE_SimplexP<1>(1), deallog.get_file_stream()); + deallog.get_file_stream() << std::endl; + print_quadrature(FE_SimplexP<1>(2), deallog.get_file_stream()); + deallog.get_file_stream() << std::endl; + print_quadrature(FE_SimplexP<2>(1), deallog.get_file_stream()); + deallog.get_file_stream() << std::endl; + print_quadrature(FE_SimplexP<2>(2), deallog.get_file_stream()); + deallog.get_file_stream() << std::endl; + + print_quadrature(FE_WedgeP<3>(1), deallog.get_file_stream()); + deallog.get_file_stream() << std::endl; + + print_quadrature(FE_DGQ<2>(1), deallog.get_file_stream()); + deallog.get_file_stream() << std::endl; + print_quadrature(FE_DGQ<2>(2), deallog.get_file_stream()); + deallog.get_file_stream() << std::endl; + print_quadrature(FE_Q<3>(1), deallog.get_file_stream()); +} diff --git a/tests/fe/fe_tools_nodal_quadrature.output b/tests/fe/fe_tools_nodal_quadrature.output new file mode 100644 index 0000000000..301a387531 --- /dev/null +++ b/tests/fe/fe_tools_nodal_quadrature.output @@ -0,0 +1,114 @@ + +FE_SimplexP<1>(1) +points: + 0.00000 + 1.00000 +weights: + 0.500000 + 0.500000 + +FE_SimplexP<1>(2) +points: + 0.00000 + 1.00000 + 0.500000 +weights: + 0.166667 + 0.166667 + 0.666667 + +FE_SimplexP<2>(1) +points: + 0.00000 0.00000 + 1.00000 0.00000 + 0.00000 1.00000 +weights: + 0.166667 + 0.166667 + 0.166667 + +FE_SimplexP<2>(2) +points: + 0.00000 0.00000 + 1.00000 0.00000 + 0.00000 1.00000 + 0.500000 0.00000 + 0.500000 0.500000 + 0.00000 0.500000 +weights: + 5.63785e-18 + -6.50521e-18 + -3.25261e-18 + 0.166667 + 0.166667 + 0.166667 + +FE_WedgeP<3>(1) +points: + 0.00000 0.00000 0.00000 + 1.00000 0.00000 0.00000 + 0.00000 1.00000 0.00000 + 0.00000 0.00000 1.00000 + 1.00000 0.00000 1.00000 + 0.00000 1.00000 1.00000 +weights: + 0.0833333 + 0.0833333 + 0.0833333 + 0.0833333 + 0.0833333 + 0.0833333 + +FE_DGQ<2>(1) +points: + 0.00000 0.00000 + 1.00000 0.00000 + 0.00000 1.00000 + 1.00000 1.00000 +weights: + 0.250000 + 0.250000 + 0.250000 + 0.250000 + +FE_DGQ<2>(2) +points: + 0.00000 0.00000 + 0.500000 0.00000 + 1.00000 0.00000 + 0.00000 0.500000 + 0.500000 0.500000 + 1.00000 0.500000 + 0.00000 1.00000 + 0.500000 1.00000 + 1.00000 1.00000 +weights: + 0.0277778 + 0.111111 + 0.0277778 + 0.111111 + 0.444444 + 0.111111 + 0.0277778 + 0.111111 + 0.0277778 + +FE_Q<3>(1) +points: + 0.00000 0.00000 0.00000 + 1.00000 0.00000 0.00000 + 0.00000 1.00000 0.00000 + 1.00000 1.00000 0.00000 + 0.00000 0.00000 1.00000 + 1.00000 0.00000 1.00000 + 0.00000 1.00000 1.00000 + 1.00000 1.00000 1.00000 +weights: + 0.125000 + 0.125000 + 0.125000 + 0.125000 + 0.125000 + 0.125000 + 0.125000 + 0.125000 diff --git a/tests/simplex/fe_p_bubbles_02.cc b/tests/simplex/fe_p_bubbles_02.cc index ea09f4e7ca..11f3a53ffb 100644 --- a/tests/simplex/fe_p_bubbles_02.cc +++ b/tests/simplex/fe_p_bubbles_02.cc @@ -20,11 +20,10 @@ #include -#include #include #include +#include #include -#include #include #include @@ -37,44 +36,6 @@ #include "../tests.h" -template -Quadrature -compute_nodal_quadrature(const FiniteElement &fe) -{ - Assert(fe.n_blocks() == 1, ExcNotImplemented()); - Assert(fe.n_components() == 1, ExcNotImplemented()); - - const ReferenceCell type = fe.reference_cell(); - - const Quadrature q_gauss = - type.get_gauss_type_quadrature(fe.tensor_degree() + 1); - Triangulation tria; - GridGenerator::reference_cell(tria, type); - const Mapping &mapping = - type.template get_default_linear_mapping(); - - auto cell = tria.begin_active(); - FEValues fe_values(mapping, - fe, - q_gauss, - update_values | update_JxW_values); - fe_values.reinit(cell); - - std::vector> nodal_quad_points = fe.get_unit_support_points(); - std::vector nodal_quad_weights(nodal_quad_points.size()); - Assert(nodal_quad_points.size() > 0, ExcNotImplemented()); - for (unsigned int i = 0; i < fe.n_dofs_per_cell(); ++i) - { - double integral = 0.0; - for (unsigned int q = 0; q < q_gauss.size(); ++q) - integral += fe_values.shape_value(i, q) * fe_values.JxW(q); - nodal_quad_weights[i] = integral; - } - return {nodal_quad_points, nodal_quad_weights}; -} - - - template void test_interpolate() @@ -179,8 +140,9 @@ test_lumped_project() DoFHandler dh(tria); dh.distribute_dofs(fe); deallog << "number of dofs = " << dh.n_dofs() << std::endl; - const Quadrature nodal_quad = compute_nodal_quadrature(fe); - const Quadrature cell_quad = QGaussSimplex( + const Quadrature nodal_quad = + FETools::compute_nodal_quadrature(fe); + const Quadrature cell_quad = QGaussSimplex( std::max(fe.tensor_degree() + 1, 2)); Vector lumped_mass(dh.n_dofs());