From: Dominik Still <dominik.still@tum.de>
Date: Fri, 5 Jul 2024 11:59:44 +0000 (+0200)
Subject: Add support for integrations on simplex unit support points
X-Git-Tag: v9.6.0-rc1~127^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=7ba096de7b5dadf863fa3e3963510eded69750ad;p=dealii.git

Add support for integrations on simplex unit support points
---

diff --git a/include/deal.II/matrix_free/util.h b/include/deal.II/matrix_free/util.h
index a78eb78ffb..532740af5e 100644
--- a/include/deal.II/matrix_free/util.h
+++ b/include/deal.II/matrix_free/util.h
@@ -22,6 +22,8 @@
 #include <deal.II/base/quadrature.h>
 #include <deal.II/base/quadrature_lib.h>
 
+#include <deal.II/fe/fe_simplex_p.h>
+
 #include <deal.II/grid/reference_cell.h>
 
 #include <deal.II/hp/q_collection.h>
@@ -56,6 +58,15 @@ namespace internal
               return {ReferenceCells::get_simplex<dim>(),
                       dealii::hp::QCollection<dim - 1>(
                         QWitherdenVincentSimplex<dim - 1>(i))};
+
+          for (unsigned int i = 1; i <= 3; ++i)
+            {
+              FE_SimplexP<dim> fe(i);
+              if (quad == Quadrature<dim>(fe.get_unit_support_points()))
+                return {ReferenceCells::get_simplex<dim>(),
+                        dealii::hp::QCollection<dim - 1>(Quadrature<dim - 1>(
+                          fe.get_unit_face_support_points()))};
+            }
         }
 
       if (dim == 3)
@@ -138,6 +149,22 @@ namespace internal
                   return {Quadrature<dim - 1>(),
                           QWitherdenVincentSimplex<dim - 1>(i)};
               }
+
+          for (unsigned int i = 1; i <= 3; ++i)
+            {
+              FE_SimplexP<dim> fe(i);
+              if (quad == Quadrature<dim>(fe.get_unit_support_points()))
+                {
+                  if (dim == 2)
+                    return {Quadrature<dim - 1>(
+                              fe.get_unit_face_support_points()), // line!
+                            Quadrature<dim - 1>()};
+                  else
+                    return {Quadrature<dim - 1>(),
+                            Quadrature<dim - 1>(
+                              fe.get_unit_face_support_points())};
+                }
+            }
         }
 
       if (dim == 3)