From d425c7d85266e9e543c439319ca43f15ac65a617 Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Thu, 30 Mar 2023 14:41:09 -0400
Subject: [PATCH] Fix nodal and Gauss-type quadratures for dim = 0.

---
 include/deal.II/grid/reference_cell.h         |  2 +-
 source/base/quadrature_lib.cc                 |  2 +
 source/grid/reference_cell.cc                 |  7 +-
 tests/grid/reference_cell_type_03_q.cc        | 10 +++
 tests/grid/reference_cell_type_03_q.output    |  1 +
 tests/simplex/nodal_type_quadrature_01.cc     | 87 +++++++++++++++++++
 tests/simplex/nodal_type_quadrature_01.output | 33 +++++++
 7 files changed, 140 insertions(+), 2 deletions(-)
 create mode 100644 tests/simplex/nodal_type_quadrature_01.cc
 create mode 100644 tests/simplex/nodal_type_quadrature_01.output

diff --git a/include/deal.II/grid/reference_cell.h b/include/deal.II/grid/reference_cell.h
index 8392edacb1..2a70b4d5a0 100644
--- a/include/deal.II/grid/reference_cell.h
+++ b/include/deal.II/grid/reference_cell.h
@@ -1226,7 +1226,7 @@ ReferenceCell::vertex(const unsigned int v) const
       case 0:
         {
           if (*this == ReferenceCells::Vertex)
-            return Point<dim>(0);
+            return Point<dim>();
           break;
         }
       case 1:
diff --git a/source/base/quadrature_lib.cc b/source/base/quadrature_lib.cc
index 3e28a1e9a0..9218b04ef0 100644
--- a/source/base/quadrature_lib.cc
+++ b/source/base/quadrature_lib.cc
@@ -2228,9 +2228,11 @@ template class QGaussSimplex<0>;
 template class QGaussSimplex<1>;
 template class QGaussSimplex<2>;
 template class QGaussSimplex<3>;
+template class QGaussWedge<0>;
 template class QGaussWedge<1>;
 template class QGaussWedge<2>;
 template class QGaussWedge<3>;
+template class QGaussPyramid<0>;
 template class QGaussPyramid<1>;
 template class QGaussPyramid<2>;
 template class QGaussPyramid<3>;
diff --git a/source/grid/reference_cell.cc b/source/grid/reference_cell.cc
index d7e7b0e9b3..2eb2aa2d51 100644
--- a/source/grid/reference_cell.cc
+++ b/source/grid/reference_cell.cc
@@ -848,7 +848,12 @@ operator>>(std::istream &in, ReferenceCell &reference_cell)
   return in;
 }
 
-
+// explicitly instantiate dimension 0 quadrature in addition to the standard
+// dimensions
+template Quadrature<0>
+ReferenceCell::get_gauss_type_quadrature(const unsigned n_points_1D) const;
+template const Quadrature<0> &
+ReferenceCell::get_nodal_type_quadrature() const;
 
 #include "reference_cell.inst"
 
diff --git a/tests/grid/reference_cell_type_03_q.cc b/tests/grid/reference_cell_type_03_q.cc
index 755f4cd225..d4975490f1 100644
--- a/tests/grid/reference_cell_type_03_q.cc
+++ b/tests/grid/reference_cell_type_03_q.cc
@@ -66,6 +66,16 @@ main()
 {
   initlog();
 
+  {
+    deallog.push("0D");
+    // It doesn't make sense to integrate in 0D, but make sure that
+    // get_gauss_type_quadrature() still works
+    deallog << "0D quadrature size: "
+            << ReferenceCells::Vertex.get_gauss_type_quadrature<0>(1).size()
+            << std::endl;
+    deallog.pop();
+  }
+
   {
     deallog.push("1D");
     test<1>(ReferenceCells::Line);
diff --git a/tests/grid/reference_cell_type_03_q.output b/tests/grid/reference_cell_type_03_q.output
index 42bbf99968..dd920473ba 100644
--- a/tests/grid/reference_cell_type_03_q.output
+++ b/tests/grid/reference_cell_type_03_q.output
@@ -1,4 +1,5 @@
 
+DEAL:0D::0D quadrature size: 1
 DEAL:1D::ReferenceCell: Line
 DEAL:1D::  computed barycenter = 0.500000
 DEAL:1D::  self-reported barycenter = 0.500000
diff --git a/tests/simplex/nodal_type_quadrature_01.cc b/tests/simplex/nodal_type_quadrature_01.cc
new file mode 100644
index 0000000000..d6c09a6243
--- /dev/null
+++ b/tests/simplex/nodal_type_quadrature_01.cc
@@ -0,0 +1,87 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 - 2023 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.
+//
+// ---------------------------------------------------------------------
+
+// Test ReferenceCell::get_nodal_type_quadrature()
+
+#include <deal.II/base/quadrature.h>
+
+#include <deal.II/grid/reference_cell.h>
+
+#include "../tests.h"
+
+using namespace dealii;
+
+template <int dim>
+void
+test(const ReferenceCell reference_cell)
+{
+  const Quadrature<dim> &quad =
+    reference_cell.template get_nodal_type_quadrature<dim>();
+
+  for (unsigned int q = 0; q < quad.size(); ++q)
+    {
+      // There is nothing to print for dim = 0 in the point
+      if (dim > 0)
+        deallog << quad.point(q) << ' ';
+      deallog << quad.weight(q) << ' ';
+      deallog << std::endl;
+    }
+}
+
+int
+main()
+{
+  initlog();
+
+  {
+    deallog.push("0d");
+    deallog.pop();
+  }
+  {
+    deallog.push("1d");
+    test<1>(ReferenceCells::Line);
+    deallog.pop();
+  }
+  {
+    deallog.push("2d-1");
+    test<2>(ReferenceCells::Triangle);
+    deallog.pop();
+  }
+  {
+    deallog.push("2d-2");
+    test<2>(ReferenceCells::Quadrilateral);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d-1");
+    test<3>(ReferenceCells::Tetrahedron);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d-2");
+    test<3>(ReferenceCells::Pyramid);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d-2");
+    test<3>(ReferenceCells::Wedge);
+    deallog.pop();
+  }
+  {
+    deallog.push("3d-3");
+    test<3>(ReferenceCells::Hexahedron);
+    deallog.pop();
+  }
+}
diff --git a/tests/simplex/nodal_type_quadrature_01.output b/tests/simplex/nodal_type_quadrature_01.output
new file mode 100644
index 0000000000..82bb892f5a
--- /dev/null
+++ b/tests/simplex/nodal_type_quadrature_01.output
@@ -0,0 +1,33 @@
+
+DEAL:1d::0.00000 inf 
+DEAL:1d::1.00000 inf 
+DEAL:2d-1::0.00000 0.00000 inf 
+DEAL:2d-1::1.00000 0.00000 inf 
+DEAL:2d-1::0.00000 1.00000 inf 
+DEAL:2d-2::0.00000 0.00000 inf 
+DEAL:2d-2::1.00000 0.00000 inf 
+DEAL:2d-2::0.00000 1.00000 inf 
+DEAL:2d-2::1.00000 1.00000 inf 
+DEAL:3d-1::0.00000 0.00000 0.00000 inf 
+DEAL:3d-1::1.00000 0.00000 0.00000 inf 
+DEAL:3d-1::0.00000 1.00000 0.00000 inf 
+DEAL:3d-1::0.00000 0.00000 1.00000 inf 
+DEAL:3d-2::-1.00000 -1.00000 0.00000 inf 
+DEAL:3d-2::1.00000 -1.00000 0.00000 inf 
+DEAL:3d-2::-1.00000 1.00000 0.00000 inf 
+DEAL:3d-2::1.00000 1.00000 0.00000 inf 
+DEAL:3d-2::0.00000 0.00000 1.00000 inf 
+DEAL:3d-2::0.00000 0.00000 0.00000 inf 
+DEAL:3d-2::1.00000 0.00000 0.00000 inf 
+DEAL:3d-2::0.00000 1.00000 0.00000 inf 
+DEAL:3d-2::0.00000 0.00000 1.00000 inf 
+DEAL:3d-2::1.00000 0.00000 1.00000 inf 
+DEAL:3d-2::0.00000 1.00000 1.00000 inf 
+DEAL:3d-3::0.00000 0.00000 0.00000 inf 
+DEAL:3d-3::1.00000 0.00000 0.00000 inf 
+DEAL:3d-3::0.00000 1.00000 0.00000 inf 
+DEAL:3d-3::1.00000 1.00000 0.00000 inf 
+DEAL:3d-3::0.00000 0.00000 1.00000 inf 
+DEAL:3d-3::1.00000 0.00000 1.00000 inf 
+DEAL:3d-3::0.00000 1.00000 1.00000 inf 
+DEAL:3d-3::1.00000 1.00000 1.00000 inf 
-- 
2.39.5