1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "mfem.hpp"
13 using namespace mfem;
14
15 #include "unit_tests.hpp"
16
17 // You typically want to start by testing things one object at a time.
18 TEST_CASE("Integration rule container with no refinement", "[IntegrationRules]")
19 {
20 // This code is automatically re-executed for all of the sections.
21 IntegrationRules my_intrules(0, Quadrature1D::GaussLegendre);
22 IntegrationRule single_point_rule(1);
23 const IntegrationRule *ir;
24
25 // The tests will be reported in these sections.
26 // Each REQUIRE counts as an assertion.
27 // true = pass, false = fail
28 SECTION("can set int rules in empty container")
29 {
30 single_point_rule[0].weight = 1.0;
31 my_intrules.SetOwnRules(0);
32
33 my_intrules.Set(Geometry::SEGMENT, 0, single_point_rule);
34 ir = &my_intrules.Get(Geometry::SEGMENT, 0);
35
36 REQUIRE(ir->Size() == 1);
37 REQUIRE(ir->IntPoint(0).weight == 1.0);
38 }
39
40 SECTION("user set int rules really are owned by the user")
41 {
42 // Set up a custom int rule and put it in the container
43 single_point_rule[0].weight = 1.0;
44 my_intrules.SetOwnRules(0);
45 my_intrules.Set(Geometry::SEGMENT, 0, single_point_rule);
46
47 // Alter the int rule
48 single_point_rule[0].weight = 2.0;
49
50 // Test ownership by making sure that the int rule has changed
51 ir = &my_intrules.Get(Geometry::SEGMENT, 0);
52 REQUIRE(ir->IntPoint(0).weight == 2.0);
53 }
54
55 SECTION("point int rules 0, 1 accessible")
56 {
57 ir = &my_intrules.Get(Geometry::POINT, 0);
58 REQUIRE(ir->Size() == 1);
59
60 ir = &my_intrules.Get(Geometry::POINT, 1);
61 REQUIRE(ir->Size() == 1);
62 }
63
64 // Can't really unit test these because it will crash due to
65 // a null pointer dereference if it doesn't work. This will
66 // force the crash in the unit tests though.
67 SECTION("resize of the SEGMENT int rule array")
68 {
69 ir = &my_intrules.Get(Geometry::SEGMENT, 100);
70 REQUIRE(true);
71 }
72
73 SECTION("setting the integration point index works")
74 {
75 ir = &my_intrules.Get(Geometry::CUBE, 5);
76 for (int i = 0; i < ir->Size(); i++)
77 {
78 REQUIRE(ir->IntPoint(i).index == i);
79 }
80 }
81
82 SECTION("intrules up to order 16 accessible")
83 {
84 for (int order = 0; order <= 16; order ++)
85 {
86 // Do this in reverse the usual order to make sure that
87 // the higher dimension cases are causing their constituent
88 // lower dimension cases to lazily create properly.
89 my_intrules.Get(Geometry::CUBE, order);
90 my_intrules.Get(Geometry::TETRAHEDRON, order);
91 my_intrules.Get(Geometry::SQUARE, order);
92 my_intrules.Get(Geometry::TRIANGLE, order);
93 my_intrules.Get(Geometry::SEGMENT, order);
94 }
95 REQUIRE(true);
96 }
97 }
98
99
poly2d(const IntegrationPoint & ip,int m,int n)100 double poly2d(const IntegrationPoint &ip, int m, int n)
101 {
102 return pow(ip.x, m)*pow(ip.y, n);
103 // exact integral over the reference triangle is
104 // m!n!/(m+n+2)! = 1/binom(m+n,m)/(m+n+1)/(m+n+2)
105 }
106
apoly2d(const IntegrationPoint & ip,int i,int j,int k)107 double apoly2d(const IntegrationPoint &ip, int i, int j, int k)
108 {
109 return pow(1. - ip.x - ip.y, i)*pow(ip.x, j)*pow(ip.y, k);
110 // exact integral over the reference triangle is (with p = i+j+k)
111 // i!j!k!/(p+2)! = 1/binom(p,i+j)/binom(i+j,i)/(p+1)/(p+2)
112 }
113
poly3d(const IntegrationPoint & ip,int l,int m,int n)114 double poly3d(const IntegrationPoint &ip, int l, int m, int n)
115 {
116 return pow(ip.x, l)*pow(ip.y, m)*pow(ip.z, n);
117 // exact integral over the reference tetrahedron is (with p = l+m+n)
118 // l!m!n!/(p+3)! = 1/binom(p,l+m)/binom(l+m,l)/(p+1)/(p+2)/(p+3)
119 }
120
121 TEST_CASE("Simplex integration rules", "[SimplexRules]")
122 {
123 // This code is automatically re-executed for all of the sections.
124 IntegrationRules my_intrules(0, Quadrature1D::GaussLegendre);
125 IntegrationRule single_point_rule(1);
126
127 const int maxn = 32;
128 int binom[maxn+1][maxn+1];
129 for (int n = 0; n <= maxn; n++)
130 {
131 binom[n][0] = binom[n][n] = 1;
132 for (int k = 1; k < n; k++)
133 {
134 binom[n][k] = binom[n-1][k] + binom[n-1][k-1];
135 }
136 }
137
138 SECTION("low triangle integration error on reference element for f=x^m y^n, where m+n <= p")
139 {
140 for (int order = 0; order <= 25; order++)
141 {
142 const IntegrationRule &ir = IntRules.Get(Geometry::TRIANGLE, order);
143
144 // using the monomial basis: x^m y^n, 0 <= m+n <= order
145 for (int p = 0; p <= order; p++)
146 {
147 for (int m = p; m >= 0; m--)
148 {
149 int n = p - m;
150
151 double integral = 0.0;
152 for (int i = 0; i < ir.GetNPoints(); i++)
153 {
154 const IntegrationPoint &ip = ir.IntPoint(i);
155 integral += ip.weight*poly2d(ip, m, n);
156 }
157
158 double exact = 1.0/binom[p][m]/(p + 1)/(p + 2);
159 double relerr = 1. - integral/exact;
160
161 // If a test fails any INFO statements preceding the REQUIRE are displayed
162 INFO("p=" << p << ", m=" << m << ", n=" << n);
163 REQUIRE(fabs(relerr) < 1e-11);
164 }
165 }
166 }
167 }
168
169 SECTION("low tet integration error on reference element for f=x^l y^m z^n, where l+m+n <= p")
170 {
171 for (int order = 0; order <= 21; order++)
172 {
173 const IntegrationRule &ir = IntRules.Get(Geometry::TETRAHEDRON, order);
174
175 for (int p = 0; p <= order; p++)
176 {
177 for (int l = p; l >= 0; l--)
178 {
179 for (int m = p - l; m >= 0; m--)
180 {
181 int n = p - l - m;
182
183 double integral = 0.0;
184 for (int i = 0; i < ir.GetNPoints(); i++)
185 {
186 const IntegrationPoint &ip = ir.IntPoint(i);
187 integral += ip.weight*poly3d(ip, l, m, n);
188 }
189
190 double exact = 1.0/binom[p][l+m]/binom[l+m][l]/(p+1)/(p+2)/(p+3);
191 double relerr = 1. - integral/exact;
192
193 // If a test fails any INFO statements preceding the REQUIRE are displayed
194 INFO("p=" << p << ", l=" << l << ", m=" << m << ", n=" << n);
195 REQUIRE(fabs(relerr) < 1e-11);
196 }
197 }
198 }
199 }
200 }
201 }
202