From: Peter Munch <peterrmuench@gmail.com>
Date: Fri, 10 Jul 2020 12:01:31 +0000 (+0200)
Subject: Implement Simplex::ScalarPolynomial and Simplex::PGauss also for 1D
X-Git-Tag: v9.3.0-rc1~1273^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F10690%2Fhead;p=dealii.git

Implement Simplex::ScalarPolynomial and Simplex::PGauss also for 1D

Reorder shape functions
---

diff --git a/include/deal.II/simplex/polynomials.h b/include/deal.II/simplex/polynomials.h
index 07c3d999e4..02f62ba434 100644
--- a/include/deal.II/simplex/polynomials.h
+++ b/include/deal.II/simplex/polynomials.h
@@ -40,8 +40,6 @@ namespace Simplex
   template <int dim>
   class ScalarPolynomial : public ScalarPolynomialsBase<dim>
   {
-    static_assert(dim == 2 || dim == 3, "Dimension not supported!");
-
   public:
     /**
      * Make the dimension available to the outside.
diff --git a/include/deal.II/simplex/quadrature_lib.h b/include/deal.II/simplex/quadrature_lib.h
index 82c9017d4c..7fc793d4eb 100644
--- a/include/deal.II/simplex/quadrature_lib.h
+++ b/include/deal.II/simplex/quadrature_lib.h
@@ -28,11 +28,14 @@ namespace Simplex
   /**
    * Integration rule for simplex entities.
    *
-   * Following number of quadrature points are currently supported:
+   * Following number of quadrature points are currently supported for 2D and
+   * 3D:
    *   - 2D: 1, 3, 7
    *   - 3D: 1, 4, 10
    *
-   *  @ingroup simplex
+   * For 1D, the quadrature rule degenerates to a `QGauss<1>(n_points)`.
+   *
+   * @ingroup simplex
    */
   template <int dim>
   class PGauss : public QSimplex<dim>
diff --git a/source/simplex/polynomials.cc b/source/simplex/polynomials.cc
index c1823cfdff..5c33c826dc 100644
--- a/source/simplex/polynomials.cc
+++ b/source/simplex/polynomials.cc
@@ -25,7 +25,14 @@ namespace Simplex
     unsigned int
     compute_n_polynomials(const unsigned int dim, const unsigned int degree)
     {
-      if (dim == 2)
+      if (dim == 1)
+        {
+          if (degree == 1)
+            return 2;
+          if (degree == 2)
+            return 3;
+        }
+      else if (dim == 2)
         {
           if (degree == 1)
             return 3;
@@ -60,22 +67,41 @@ namespace Simplex
   ScalarPolynomial<dim>::compute_value(const unsigned int i,
                                        const Point<dim> & p) const
   {
-    if (dim == 2)
+    if (dim == 1)
       {
         if (this->degree() == 1)
           {
             if (i == 0)
+              return 1.0 - p[0];
+            else if (i == 1)
               return p[0];
+          }
+        else if (this->degree() == 2)
+          {
+            if (i == 0)
+              return 2.0 * p[0] * p[0] - 3.0 * p[0] + 1;
             else if (i == 1)
-              return p[1];
+              return 2.0 * p[0] * p[0] - p[0];
             else if (i == 2)
+              return -4.0 * p[0] * p[0] + 4.0 * p[0];
+          }
+      }
+    else if (dim == 2)
+      {
+        if (this->degree() == 1)
+          {
+            if (i == 0)
               return 1.0 - p[0] - p[1];
+            else if (i == 1)
+              return p[0];
+            else if (i == 2)
+              return p[1];
           }
         else if (this->degree() == 2)
           {
-            const double t1 = p[0];
-            const double t2 = p[1];
-            const double t3 = 1.0 - p[0] - p[1];
+            const double t1 = 1.0 - p[0] - p[1];
+            const double t2 = p[0];
+            const double t3 = p[1];
 
             if (i == 0)
               return t1 * (2.0 * t1 - 1.0);
@@ -147,24 +173,47 @@ namespace Simplex
   {
     Tensor<1, dim> grad;
 
-    if (dim == 2)
+    if (dim == 1)
+      {
+        if (this->degree() == 1)
+          {
+            if (i == 0)
+              grad[0] = -1.0;
+            else if (i == 1)
+              grad[0] = 1.0;
+          }
+        else if (this->degree() == 2)
+          {
+            if (i == 0)
+              grad[0] = 4.0 * p[0] - 3.0;
+            else if (i == 1)
+              grad[0] = 4.0 * p[0] - 1.0;
+            else if (i == 2)
+              grad[0] = -8.0 * p[0] + 4.0;
+          }
+        else
+          {
+            Assert(false, ExcNotImplemented());
+          }
+      }
+    else if (dim == 2)
       {
         if (this->degree() == 1)
           {
             if (i == 0)
               {
-                grad[0] = +1.0;
-                grad[1] = +0.0;
+                grad[0] = -1.0;
+                grad[1] = -1.0;
               }
             else if (i == 1)
               {
-                grad[0] = +0.0;
-                grad[1] = +1.0;
+                grad[0] = +1.0;
+                grad[1] = +0.0;
               }
             else if (i == 2)
               {
-                grad[0] = -1.0;
-                grad[1] = -1.0;
+                grad[0] = +0.0;
+                grad[1] = +1.0;
               }
             else
               {
@@ -173,32 +222,32 @@ namespace Simplex
           }
         else if (this->degree() == 2)
           {
-            if (i == 2)
+            if (i == 0)
               {
                 grad[0] = -3.0 + 4.0 * (p[0] + p[1]);
                 grad[1] = -3.0 + 4.0 * (p[0] + p[1]);
               }
-            else if (i == 0)
+            else if (i == 1)
               {
                 grad[0] = 4.0 * p[0] - 1.0;
                 grad[1] = 0.0;
               }
-            else if (i == 1)
+            else if (i == 2)
               {
                 grad[0] = 0.0;
                 grad[1] = 4.0 * p[1] - 1.0;
               }
-            else if (i == 5)
+            else if (i == 3)
               {
                 grad[0] = 4.0 * (1.0 - 2.0 * p[0] - p[1]);
                 grad[1] = -4.0 * p[0];
               }
-            else if (i == 3)
+            else if (i == 4)
               {
                 grad[0] = 4.0 * p[1];
                 grad[1] = 4.0 * p[0];
               }
-            else if (i == 4)
+            else if (i == 5)
               {
                 grad[0] = -4.0 * p[1];
                 grad[1] = 4.0 * (1.0 - p[0] - 2.0 * p[1]);
@@ -436,6 +485,7 @@ namespace Simplex
     return std::make_unique<ScalarPolynomial<dim>>(*this);
   }
 
+  template class ScalarPolynomial<1>;
   template class ScalarPolynomial<2>;
   template class ScalarPolynomial<3>;
 
diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc
index 94bd3ad9c6..b0ffad033b 100644
--- a/source/simplex/quadrature_lib.cc
+++ b/source/simplex/quadrature_lib.cc
@@ -33,8 +33,14 @@ namespace Simplex
     : QSimplex<dim>(Quadrature<dim>())
   {
     // fill quadrature points and quadrature weights
+    if (dim == 1)
+      {
+        const dealii::QGauss<dim> quad(n_points);
 
-    if (dim == 2)
+        this->quadrature_points = quad.get_points();
+        this->weights           = quad.get_weights();
+      }
+    else if (dim == 2)
       {
         if (n_points == 1)
           {
@@ -140,6 +146,7 @@ namespace Simplex
 } // namespace Simplex
 
 
+template class Simplex::PGauss<1>;
 template class Simplex::PGauss<2>;
 template class Simplex::PGauss<3>;
 
diff --git a/tests/simplex/polynomials_01.cc b/tests/simplex/polynomials_01.cc
index 86c276ef9e..eb8c82cfcc 100644
--- a/tests/simplex/polynomials_01.cc
+++ b/tests/simplex/polynomials_01.cc
@@ -63,6 +63,16 @@ main()
 {
   initlog();
 
+  {
+    deallog.push("1d-1");
+    test<1>(1);
+    deallog.pop();
+  }
+  {
+    deallog.push("1d-2");
+    test<1>(2);
+    deallog.pop();
+  }
   {
     deallog.push("2d-1");
     test<2>(1);
diff --git a/tests/simplex/polynomials_01.output b/tests/simplex/polynomials_01.output
index fe28ac73ea..5d8db3cf7e 100644
--- a/tests/simplex/polynomials_01.output
+++ b/tests/simplex/polynomials_01.output
@@ -1,22 +1,32 @@
 
-DEAL:2d-1::0.211325 0.211325 0.577350 
-DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 
-DEAL:2d-1::0.788675 0.211325 2.77556e-17 
-DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 
-DEAL:2d-1::0.211325 0.788675 0.00000 
-DEAL:2d-1::1.00000 0.00000 0.00000 1.00000 -1.00000 -1.00000 
-DEAL:2d-2::-0.0872983 -0.0872983 0.425403 0.0508067 0.349193 0.349193 
-DEAL:2d-2::-0.549193 0.00000 0.00000 -0.549193 -2.09839 -2.09839 0.450807 0.450807 -0.450807 2.64758 2.64758 -0.450807 
-DEAL:2d-2::0.00000 -0.0872983 -0.0872983 0.225403 0.174597 0.774597 
-DEAL:2d-2::1.00000 0.00000 0.00000 -0.549193 -0.549193 -0.549193 0.450807 2.00000 -0.450807 1.09839 -0.450807 -2.00000 
-DEAL:2d-2::0.687298 -0.0872983 1.38778e-17 0.400000 -6.25620e-18 -4.92550e-17 
-DEAL:2d-2::2.54919 0.00000 0.00000 -0.549193 1.00000 1.00000 0.450807 3.54919 -0.450807 -0.450807 -3.54919 -3.54919 
-DEAL:2d-2::-0.0872983 0.00000 -0.0872983 0.225403 0.774597 0.174597 
-DEAL:2d-2::-0.549193 0.00000 0.00000 1.00000 -0.549193 -0.549193 2.00000 0.450807 -2.00000 -0.450807 1.09839 -0.450807 
-DEAL:2d-2::0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 
-DEAL:2d-2::1.00000 0.00000 0.00000 1.00000 1.00000 1.00000 2.00000 2.00000 -2.00000 -2.00000 -2.00000 -2.00000 
-DEAL:2d-2::-0.0872983 0.687298 0.00000 0.400000 0.00000 0.00000 
-DEAL:2d-2::-0.549193 0.00000 0.00000 2.54919 1.00000 1.00000 3.54919 0.450807 -3.54919 -3.54919 -0.450807 -0.450807 
+DEAL:1d-1::0.788675 0.211325 
+DEAL:1d-1::-1.00000 1.00000 
+DEAL:1d-1::0.211325 0.788675 
+DEAL:1d-1::-1.00000 1.00000 
+DEAL:1d-2::0.687298 -0.0872983 0.400000 
+DEAL:1d-2::-2.54919 -0.549193 3.09839 
+DEAL:1d-2::0.00000 0.00000 1.00000 
+DEAL:1d-2::-1.00000 1.00000 0.00000 
+DEAL:1d-2::-0.0872983 0.687298 0.400000 
+DEAL:1d-2::0.549193 2.54919 -3.09839 
+DEAL:2d-1::0.577350 0.211325 0.211325 
+DEAL:2d-1::-1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:2d-1::2.77556e-17 0.788675 0.211325 
+DEAL:2d-1::-1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:2d-1::0.00000 0.211325 0.788675 
+DEAL:2d-1::-1.00000 -1.00000 1.00000 0.00000 0.00000 1.00000 
+DEAL:2d-2::0.425403 -0.0872983 -0.0872983 0.349193 0.0508067 0.349193 
+DEAL:2d-2::-2.09839 -2.09839 -0.549193 0.00000 0.00000 -0.549193 2.64758 -0.450807 0.450807 0.450807 -0.450807 2.64758 
+DEAL:2d-2::-0.0872983 0.00000 -0.0872983 0.774597 0.225403 0.174597 
+DEAL:2d-2::-0.549193 -0.549193 1.00000 0.00000 0.00000 -0.549193 -0.450807 -2.00000 0.450807 2.00000 -0.450807 1.09839 
+DEAL:2d-2::1.38778e-17 0.687298 -0.0872983 -4.92550e-17 0.400000 -6.25620e-18 
+DEAL:2d-2::1.00000 1.00000 2.54919 0.00000 0.00000 -0.549193 -3.54919 -3.54919 0.450807 3.54919 -0.450807 -0.450807 
+DEAL:2d-2::-0.0872983 -0.0872983 0.00000 0.174597 0.225403 0.774597 
+DEAL:2d-2::-0.549193 -0.549193 -0.549193 0.00000 0.00000 1.00000 1.09839 -0.450807 2.00000 0.450807 -2.00000 -0.450807 
+DEAL:2d-2::0.00000 0.00000 0.00000 0.00000 1.00000 0.00000 
+DEAL:2d-2::1.00000 1.00000 1.00000 0.00000 0.00000 1.00000 -2.00000 -2.00000 2.00000 2.00000 -2.00000 -2.00000 
+DEAL:2d-2::0.00000 -0.0872983 0.687298 0.00000 0.400000 0.00000 
+DEAL:2d-2::1.00000 1.00000 -0.549193 0.00000 0.00000 2.54919 -0.450807 -0.450807 3.54919 0.450807 -3.54919 -3.54919 
 DEAL:3d-1::0.366025 0.211325 0.211325 0.211325 
 DEAL:3d-1::-1.00000 -1.00000 -1.00000 1.00000 0.00000 0.00000 0.00000 1.00000 0.00000 0.00000 0.00000 1.00000 
 DEAL:3d-2::0.214315 -0.0872983 -0.0872983 -0.0872983 0.298387 0.0508067 0.298387 0.298387 0.0508067 0.0508067 
diff --git a/tests/simplex/quadrature_lib_01.cc b/tests/simplex/quadrature_lib_01.cc
index a59a34fd6d..8e9f8c3002 100644
--- a/tests/simplex/quadrature_lib_01.cc
+++ b/tests/simplex/quadrature_lib_01.cc
@@ -42,9 +42,24 @@ main()
 {
   initlog();
 
+  {
+    deallog.push("1d-1");
+    test<1>(1 /*n_points*/);
+    deallog.pop();
+  }
+  {
+    deallog.push("1d-2");
+    test<1>(2);
+    deallog.pop();
+  }
+  {
+    deallog.push("1d-3");
+    test<1>(2);
+    deallog.pop();
+  }
   {
     deallog.push("2d-1");
-    test<2>(1 /*n_points*/);
+    test<2>(1);
     deallog.pop();
   }
   {
diff --git a/tests/simplex/quadrature_lib_01.output b/tests/simplex/quadrature_lib_01.output
index 03258812b2..24a8420bbd 100644
--- a/tests/simplex/quadrature_lib_01.output
+++ b/tests/simplex/quadrature_lib_01.output
@@ -1,4 +1,9 @@
 
+DEAL:1d-1::0.500000 1.00000 
+DEAL:1d-2::0.211325 0.500000 
+DEAL:1d-2::0.788675 0.500000 
+DEAL:1d-3::0.211325 0.500000 
+DEAL:1d-3::0.788675 0.500000 
 DEAL:2d-1::0.333333 0.333333 1.00000 
 DEAL:2d-3::0.666667 0.166667 0.166667 
 DEAL:2d-3::0.166667 0.666667 0.166667