From f9e551f002937206704f0f7dfbcb1c3557f6afec Mon Sep 17 00:00:00 2001 From: Brandon Talamini Date: Sun, 3 Dec 2023 07:20:53 -0800 Subject: [PATCH] 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. --- optimism/QuadratureRule.py | 54 ++++++++++++++++++++++++++++ optimism/test/test_QuadratureRule.py | 17 ++++----- 2 files changed, 63 insertions(+), 8 deletions(-) diff --git a/optimism/QuadratureRule.py b/optimism/QuadratureRule.py index 5e9e577..d14948b 100644 --- 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 index 9d4d645..d0c89fb 100644 --- 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()