Add a higher order quadrature rule (exact up to p=10)

Also: now a ValueError is emitted if a user asks for a quadrature rule precision that is above what we have implemented. Before the error was a bit obscure. Improved the unit testing of the quadrature rules by making subtests for every polynomial degree check. This way, if the test fails, you can figure out which quadrature rules are broken.
sandialabs · Dec 3, 2023 · fbf927a · fbf927a
1 parent d81daee
commit fbf927a
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()