Merge pull request #73 from btalamini/feature/quadrature_precision_10

Add a higher order quadrature rule (exact up to p=10)
sandialabs · Dec 3, 2023 · 56b43b8 · 56b43b8
2 parents d00ffd1 + f9e551f
commit 56b43b8
Show file tree

Hide file tree

Showing 2 changed files with 63 additions and 8 deletions.
diff --git a/optimism/QuadratureRule.py b/optimism/QuadratureRule.py
@@ -131,6 +131,60 @@ def create_quadrature_rule_on_triangle(degree):
                        4.14255378091870E-02,
                        4.14255378091870E-02,
                        4.14255378091870E-02])
+    elif degree <= 10:
+        xi = onp.array([[0.33333333333333333E+00, 0.33333333333333333E+00],
+                        [0.4269134091050342E-02,  0.49786543295447483E+00],
+                        [0.49786543295447483E+00, 0.4269134091050342E-02],
+                        [0.49786543295447483E+00, 0.49786543295447483E+00],
+                        [0.14397510054188759E+00, 0.42801244972905617E+00],
+                        [0.42801244972905617E+00, 0.14397510054188759E+00],
+                        [0.42801244972905617E+00, 0.42801244972905617E+00],
+                        [0.6304871745135507E+00,  0.18475641274322457E+00],
+                        [0.18475641274322457E+00, 0.6304871745135507E+00],
+                        [0.18475641274322457E+00, 0.18475641274322457E+00],
+                        [0.9590375628566448E+00,  0.20481218571677562E-01],
+                        [0.20481218571677562E-01, 0.9590375628566448E+00],
+                        [0.20481218571677562E-01, 0.20481218571677562E-01],
+                        [0.3500298989727196E-01,  0.1365735762560334E+00],
+                        [0.1365735762560334E+00,  0.8284234338466947E+00],
+                        [0.8284234338466947E+00,  0.3500298989727196E-01],
+                        [0.1365735762560334E+00,  0.3500298989727196E-01],
+                        [0.8284234338466947E+00,  0.1365735762560334E+00],
+                        [0.3500298989727196E-01,  0.8284234338466947E+00],
+                        [0.37549070258442674E-01, 0.3327436005886386E+00],
+                        [0.3327436005886386E+00,  0.6297073291529187E+00],
+                        [0.6297073291529187E+00,  0.37549070258442674E-01],
+                        [0.3327436005886386E+00,  0.37549070258442674E-01],
+                        [0.6297073291529187E+00,  0.3327436005886386E+00],
+                        [0.37549070258442674E-01, 0.6297073291529187E+00]])
+
+        w = onp.array([0.4176169990259819E-01,
+                       0.36149252960283717E-02,
+                       0.36149252960283717E-02,
+                       0.36149252960283717E-02,
+                       0.3724608896049025E-01,
+                       0.3724608896049025E-01,
+                       0.3724608896049025E-01,
+                       0.39323236701554264E-01,
+                       0.39323236701554264E-01,
+                       0.39323236701554264E-01,
+                       0.3464161543553752E-02,
+                       0.3464161543553752E-02,
+                       0.3464161543553752E-02,
+                       0.147591601673897E-01,
+                       0.147591601673897E-01,
+                       0.147591601673897E-01,
+                       0.147591601673897E-01,
+                       0.147591601673897E-01,
+                       0.147591601673897E-01,
+                       0.1978968359803062E-01,
+                       0.1978968359803062E-01,
+                       0.1978968359803062E-01,
+                       0.1978968359803062E-01,
+                       0.1978968359803062E-01,
+                       0.1978968359803062E-01])
+    else:
+        raise ValueError("Quadrature of precision this high is not implemented.")
 
     return QuadratureRule(xi, w)
 

diff --git a/optimism/test/test_QuadratureRule.py b/optimism/test/test_QuadratureRule.py
@@ -43,7 +43,7 @@ def are_positive_weights(QuadratureRuleFactory, degree):
 
 class TestQuadratureRules(TestFixture.TestFixture):
     endpoints = (0.0, 1.0)     # better if the quadrature rule module provided this
-    max_degree_2D = 6
+    max_degree_2D = 10
     max_degree_1D = 25
 
 
@@ -87,13 +87,14 @@ def test_triangle_quadrature_points_in_domain(self):
 
     def test_triangle_quadrature_exactness(self):
         for degree in range(self.max_degree_2D + 1):
-            qr = QuadratureRule.create_quadrature_rule_on_triangle(degree)
-            for i in range(degree + 1):
-                for j in range(degree + 1 - i):
-                    monomial = qr.xigauss[:,0]**i * qr.xigauss[:,1]**j
-                    quadratureAnswer = np.sum(monomial * qr.wgauss)
-                    exactAnswer = integrate_2D_monomial_on_triangle(i, j)
-                    self.assertNear(quadratureAnswer, exactAnswer, 14)
+            with self.subTest(i=degree):
+                qr = QuadratureRule.create_quadrature_rule_on_triangle(degree)
+                for i in range(degree + 1):
+                    for j in range(degree + 1 - i):
+                        monomial = qr.xigauss[:,0]**i * qr.xigauss[:,1]**j
+                        quadratureAnswer = np.sum(monomial * qr.wgauss)
+                        exactAnswer = integrate_2D_monomial_on_triangle(i, j)
+                        self.assertNear(quadratureAnswer, exactAnswer, 14)
 
 if __name__ == '__main__':
     unittest.main()