From b7e0a2cbc0a9e0ffdcbbd6d1988a08b89d3fb85e Mon Sep 17 00:00:00 2001
From: David Wells <drwells@email.unc.edu>
Date: Fri, 22 Jan 2021 09:10:52 -0500
Subject: [PATCH] Add another simplex quadrature rule.

---
 source/simplex/quadrature_lib.cc | 118 ++++++++++++++++++++++++++++++-
 1 file changed, 117 insertions(+), 1 deletion(-)

diff --git a/source/simplex/quadrature_lib.cc b/source/simplex/quadrature_lib.cc
index 6e8b9832ea..55d8c9eb69 100644
--- a/source/simplex/quadrature_lib.cc
+++ b/source/simplex/quadrature_lib.cc
@@ -82,6 +82,44 @@ namespace Simplex
             this->weights.emplace_back(q12 * 0.132394152789);
             this->weights.emplace_back(q12 * 0.132394152789);
           }
+        else if (n_points_1D == 4)
+          {
+            const double q12 = 0.5;
+
+            // clang-format off
+            this->quadrature_points.emplace_back(0.06493051315916486, 0.06493051315916486);
+            this->quadrature_points.emplace_back(0.06493051315916486, 0.8701389736816703);
+            this->quadrature_points.emplace_back(0.8701389736816703, 0.06493051315916486);
+            this->quadrature_points.emplace_back(0.2845755842491703, 0.517039939069323);
+            this->quadrature_points.emplace_back(0.3135591843849315, 0.04386347179237249);
+            this->quadrature_points.emplace_back(0.19838447668150672, 0.2845755842491703);
+            this->quadrature_points.emplace_back(0.642577343822696, 0.3135591843849315);
+            this->quadrature_points.emplace_back(0.517039939069323, 0.19838447668150672);
+            this->quadrature_points.emplace_back(0.04386347179237249, 0.642577343822696);
+            this->quadrature_points.emplace_back(0.517039939069323, 0.2845755842491703);
+            this->quadrature_points.emplace_back(0.04386347179237249, 0.3135591843849315);
+            this->quadrature_points.emplace_back(0.19838447668150672, 0.517039939069323);
+            this->quadrature_points.emplace_back(0.642577343822696, 0.04386347179237249);
+            this->quadrature_points.emplace_back(0.2845755842491703, 0.19838447668150672);
+            this->quadrature_points.emplace_back(0.3135591843849315, 0.642577343822696);
+            // clang-format on
+
+            this->weights.emplace_back(q12 * 0.053077801790232484);
+            this->weights.emplace_back(q12 * 0.053077801790232484);
+            this->weights.emplace_back(q12 * 0.053077801790232484);
+            this->weights.emplace_back(q12 * 0.07085308369213365);
+            this->weights.emplace_back(q12 * 0.069274682079417);
+            this->weights.emplace_back(q12 * 0.07085308369213365);
+            this->weights.emplace_back(q12 * 0.069274682079417);
+            this->weights.emplace_back(q12 * 0.07085308369213365);
+            this->weights.emplace_back(q12 * 0.069274682079417);
+            this->weights.emplace_back(q12 * 0.07085308369213365);
+            this->weights.emplace_back(q12 * 0.069274682079417);
+            this->weights.emplace_back(q12 * 0.07085308369213365);
+            this->weights.emplace_back(q12 * 0.069274682079417);
+            this->weights.emplace_back(q12 * 0.07085308369213365);
+            this->weights.emplace_back(q12 * 0.069274682079417);
+          }
       }
     else if (dim == 3)
       {
@@ -136,13 +174,91 @@ namespace Simplex
             this->weights.emplace_back(0.0214899534130631 * Q16);
             this->weights.emplace_back(0.0214899534130631 * Q16);
           }
+        else if (n_points_1D == 4)
+          {
+            const double Q16 = 1.0 / 6.0;
+
+            // clang-format off
+            this->quadrature_points.emplace_back(0.25, 0.25, 0.25);
+            this->quadrature_points.emplace_back(0.3157011497782028, 0.3157011497782028, 0.3157011497782028);
+            this->quadrature_points.emplace_back(0.3157011497782028, 0.3157011497782028, 0.05289655066539156);
+            this->quadrature_points.emplace_back(0.3157011497782028, 0.05289655066539156, 0.3157011497782028);
+            this->quadrature_points.emplace_back(0.05289655066539156, 0.3157011497782028, 0.3157011497782028);
+            this->quadrature_points.emplace_back(0.05048982259839635, 0.05048982259839635, 0.44951017740160365);
+            this->quadrature_points.emplace_back(0.05048982259839635, 0.44951017740160365, 0.05048982259839635);
+            this->quadrature_points.emplace_back(0.44951017740160365, 0.05048982259839635, 0.05048982259839635);
+            this->quadrature_points.emplace_back(0.05048982259839635, 0.44951017740160365, 0.44951017740160365);
+            this->quadrature_points.emplace_back(0.44951017740160365, 0.05048982259839635, 0.44951017740160365);
+            this->quadrature_points.emplace_back(0.44951017740160365, 0.44951017740160365, 0.05048982259839635);
+            this->quadrature_points.emplace_back(0.188833831026001, 0.188833831026001, 0.575171637587);
+            this->quadrature_points.emplace_back(0.02126547254148314, 0.02126547254148314, 0.8108302410985486);
+            this->quadrature_points.emplace_back(0.188833831026001, 0.575171637587, 0.188833831026001);
+            this->quadrature_points.emplace_back(0.02126547254148314, 0.8108302410985486, 0.02126547254148314);
+            this->quadrature_points.emplace_back(0.575171637587, 0.188833831026001, 0.188833831026001);
+            this->quadrature_points.emplace_back(0.8108302410985486, 0.02126547254148314, 0.02126547254148314);
+            this->quadrature_points.emplace_back(0.188833831026001, 0.575171637587, 0.047160700360998065);
+            this->quadrature_points.emplace_back(0.02126547254148314, 0.8108302410985486, 0.1466388138184851);
+            this->quadrature_points.emplace_back(0.575171637587, 0.188833831026001, 0.047160700360998065);
+            this->quadrature_points.emplace_back(0.8108302410985486, 0.02126547254148314, 0.1466388138184851);
+            this->quadrature_points.emplace_back(0.575171637587, 0.047160700360998065, 0.188833831026001);
+            this->quadrature_points.emplace_back(0.8108302410985486, 0.1466388138184851, 0.02126547254148314);
+            this->quadrature_points.emplace_back(0.188833831026001, 0.188833831026001, 0.047160700360998065);
+            this->quadrature_points.emplace_back(0.02126547254148314, 0.02126547254148314, 0.1466388138184851);
+            this->quadrature_points.emplace_back(0.188833831026001, 0.047160700360998065, 0.188833831026001);
+            this->quadrature_points.emplace_back(0.02126547254148314, 0.1466388138184851, 0.02126547254148314);
+            this->quadrature_points.emplace_back(0.047160700360998065, 0.188833831026001, 0.188833831026001);
+            this->quadrature_points.emplace_back(0.1466388138184851, 0.02126547254148314, 0.02126547254148314);
+            this->quadrature_points.emplace_back(0.188833831026001, 0.047160700360998065, 0.575171637587);
+            this->quadrature_points.emplace_back(0.02126547254148314, 0.1466388138184851, 0.8108302410985486);
+            this->quadrature_points.emplace_back(0.047160700360998065, 0.188833831026001, 0.575171637587);
+            this->quadrature_points.emplace_back(0.1466388138184851, 0.02126547254148314, 0.8108302410985486);
+            this->quadrature_points.emplace_back(0.047160700360998065, 0.575171637587, 0.188833831026001);
+            this->quadrature_points.emplace_back(0.1466388138184851, 0.8108302410985486, 0.02126547254148314);
+            // clang-format on
+
+            this->weights.emplace_back(Q16 * 0.09548528946413085);
+            this->weights.emplace_back(Q16 * 0.04232958120996703);
+            this->weights.emplace_back(Q16 * 0.04232958120996703);
+            this->weights.emplace_back(Q16 * 0.04232958120996703);
+            this->weights.emplace_back(Q16 * 0.04232958120996703);
+            this->weights.emplace_back(Q16 * 0.03189692783285758);
+            this->weights.emplace_back(Q16 * 0.03189692783285758);
+            this->weights.emplace_back(Q16 * 0.03189692783285758);
+            this->weights.emplace_back(Q16 * 0.03189692783285758);
+            this->weights.emplace_back(Q16 * 0.03189692783285758);
+            this->weights.emplace_back(Q16 * 0.03189692783285758);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+            this->weights.emplace_back(Q16 * 0.03720713072833462);
+            this->weights.emplace_back(Q16 * 0.008110770829903342);
+          }
       }
 
     AssertDimension(this->quadrature_points.size(), this->weights.size());
     Assert(this->quadrature_points.size() > 0,
            ExcNotImplemented(
              "Simplex::QGauss is currently only implemented for "
-             "n_points_1D = 1, 2, and 3, while you are asking for "
+             "n_points_1D = 1, 2, 3, and 4 while you are asking for "
              "n_points_1D = " +
              Utilities::to_string(n_points_1D)));
   }
-- 
2.39.5