]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add FETools::compute_nodal_quadrature(). 17214/head
authorDavid Wells <drwells@email.unc.edu>
Fri, 5 Jul 2024 18:09:44 +0000 (14:09 -0400)
committerDavid Wells <drwells@email.unc.edu>
Sat, 6 Jul 2024 12:40:00 +0000 (08:40 -0400)
doc/news/changes/minor/20240705DavidWells [new file with mode: 0644]
include/deal.II/fe/fe_tools.h
include/deal.II/fe/fe_tools.templates.h
source/fe/fe_tools.inst.in
tests/fe/fe_tools_nodal_quadrature.cc [new file with mode: 0644]
tests/fe/fe_tools_nodal_quadrature.output [new file with mode: 0644]
tests/simplex/fe_p_bubbles_02.cc

diff --git a/doc/news/changes/minor/20240705DavidWells b/doc/news/changes/minor/20240705DavidWells
new file mode 100644 (file)
index 0000000..043a06d
--- /dev/null
@@ -0,0 +1,5 @@
+New: Added a function FETools::compute_nodal_quadrature()
+which computes the nodal Quadrature rule for an interpolatory
+FiniteElement.
+<br>
+(David Wells, 2023/07/05)
index bdddae6fe71f304a2753da7d00628305bd4fc450..6897a9c1d58c2933f2136f2d4b2652daed87380a 100644 (file)
@@ -586,7 +586,41 @@ namespace FETools
     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
index f499d71169fba915153613c3bf5fb4874023dd67..e7636e5370dfd6cb5f3086bb7eaf0c87a7f99a7c 100644 (file)
@@ -2947,6 +2947,46 @@ namespace FETools
 
 
 
+  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
index 7071c7e7cdf5e692b9214d876e1cbc0bde4a9963..a2bff7ba44059c4b985132a44dfc8634338f761e 100644 (file)
@@ -170,6 +170,10 @@ for (deal_II_dimension : DIMENSIONS; deal_II_space_dimension : SPACE_DIMENSIONS)
         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
diff --git a/tests/fe/fe_tools_nodal_quadrature.cc b/tests/fe/fe_tools_nodal_quadrature.cc
new file mode 100644 (file)
index 0000000..2555a74
--- /dev/null
@@ -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 <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());
+}
diff --git a/tests/fe/fe_tools_nodal_quadrature.output b/tests/fe/fe_tools_nodal_quadrature.output
new file mode 100644 (file)
index 0000000..301a387
--- /dev/null
@@ -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
index ea09f4e7cacacb687c8f9911883322e12e8505d7..11f3a53ffb3a87e78dd42e726c14f01783294391 100644 (file)
 
 #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()
@@ -179,8 +140,9 @@ test_lumped_project()
           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());

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.