--- /dev/null
+New: Added a function FETools::compute_nodal_quadrature()
+which computes the nodal Quadrature rule for an interpolatory
+FiniteElement.
+<br>
+(David Wells, 2023/07/05)
const unsigned int face,
FullMatrix<double> &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 <int dim, int spacedim = dim>
+ Quadrature<dim>
+ compute_nodal_quadrature(const FiniteElement<dim, spacedim> &fe);
/**
* Wrapper around
+ template <int dim, int spacedim>
+ Quadrature<dim>
+ compute_nodal_quadrature(const FiniteElement<dim, spacedim> &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<dim> q_gauss =
+ type.get_gauss_type_quadrature<dim>(fe.tensor_degree() + 1);
+ Triangulation<dim, spacedim> tria;
+ GridGenerator::reference_cell(tria, type);
+ const Mapping<dim, spacedim> &mapping =
+ type.template get_default_linear_mapping<dim, spacedim>();
+
+ auto cell = tria.begin_active();
+ FEValues<dim, spacedim> fe_values(mapping,
+ fe,
+ q_gauss,
+ update_values | update_JxW_values);
+ fe_values.reinit(cell);
+
+ const std::vector<Point<dim>> &nodal_quad_points =
+ fe.get_unit_support_points();
+ std::vector<double> 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
const Quadrature<deal_II_dimension> &,
const Quadrature<deal_II_dimension> &,
FullMatrix<double> &);
+
+ template Quadrature<deal_II_dimension>
+ compute_nodal_quadrature(
+ const FiniteElement<deal_II_dimension, deal_II_space_dimension> &);
#endif
#if deal_II_dimension == deal_II_space_dimension
--- /dev/null
+// ------------------------------------------------------------------------
+//
+// 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 <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_simplex_p.h>
+#include <deal.II/fe/fe_tools.h>
+#include <deal.II/fe/fe_wedge_p.h>
+
+#include "../tests.h"
+
+template <int dim>
+void
+print_quadrature(const FE_Q_Base<dim> &fe, std::ostream &out)
+{
+ out << fe.get_name() << std::endl;
+ const Quadrature<dim> 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<dim> lobatto = QGaussLobatto<dim>(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 <int dim>
+void
+print_quadrature(const FiniteElement<dim> &fe, std::ostream &out)
+{
+ out << fe.get_name() << std::endl;
+ const Quadrature<dim> 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());
+}
--- /dev/null
+
+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
#include <deal.II/dofs/dof_handler.h>
-#include <deal.II/fe/fe_pyramid_p.h>
#include <deal.II/fe/fe_simplex_p.h>
#include <deal.II/fe/fe_simplex_p_bubbles.h>
+#include <deal.II/fe/fe_tools.h>
#include <deal.II/fe/fe_values.h>
-#include <deal.II/fe/fe_wedge_p.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_tools.h>
#include "../tests.h"
-template <int dim, int spacedim = dim>
-Quadrature<dim>
-compute_nodal_quadrature(const FiniteElement<dim, spacedim> &fe)
-{
- Assert(fe.n_blocks() == 1, ExcNotImplemented());
- Assert(fe.n_components() == 1, ExcNotImplemented());
-
- const ReferenceCell type = fe.reference_cell();
-
- const Quadrature<dim> q_gauss =
- type.get_gauss_type_quadrature<dim>(fe.tensor_degree() + 1);
- Triangulation<dim, spacedim> tria;
- GridGenerator::reference_cell(tria, type);
- const Mapping<dim, spacedim> &mapping =
- type.template get_default_linear_mapping<dim, spacedim>();
-
- auto cell = tria.begin_active();
- FEValues<dim, spacedim> fe_values(mapping,
- fe,
- q_gauss,
- update_values | update_JxW_values);
- fe_values.reinit(cell);
-
- std::vector<Point<dim>> nodal_quad_points = fe.get_unit_support_points();
- std::vector<double> 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 <int dim, int spacedim = dim>
void
test_interpolate()
DoFHandler<dim, spacedim> dh(tria);
dh.distribute_dofs(fe);
deallog << "number of dofs = " << dh.n_dofs() << std::endl;
- const Quadrature<dim> nodal_quad = compute_nodal_quadrature(fe);
- const Quadrature<dim> cell_quad = QGaussSimplex<dim>(
+ const Quadrature<dim> nodal_quad =
+ FETools::compute_nodal_quadrature(fe);
+ const Quadrature<dim> cell_quad = QGaussSimplex<dim>(
std::max<unsigned int>(fe.tensor_degree() + 1, 2));
Vector<double> lumped_mass(dh.n_dofs());
projects / dealii.git / commitdiff