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 #include "unit_tests.hpp"
14
15 #include <iostream>
16 #include <cmath>
17
18 using namespace mfem;
19
20 /**
21 * Utility function to generate IntegerationPoints, based on param ip
22 * that are outside the unit interval. Results are placed in output
23 * parameter arr.
24 */
GetRelatedIntegrationPoints(const IntegrationPoint & ip,int dim,Array<IntegrationPoint> & arr)25 void GetRelatedIntegrationPoints(const IntegrationPoint& ip, int dim,
26 Array<IntegrationPoint>& arr)
27 {
28 IntegrationPoint pt = ip;
29 int idx = 0;
30
31 switch (dim)
32 {
33 case 1:
34 arr.SetSize(3);
35
36 pt.x = ip.x; arr[idx++] = pt;
37 pt.x = -ip.x; arr[idx++] = pt;
38 pt.x = 1+ip.x; arr[idx++] = pt;
39 break;
40 case 2:
41 arr.SetSize(7);
42
43 pt.Set2( ip.x, ip.y); arr[idx++] = pt;
44 pt.Set2( -ip.x, ip.y); arr[idx++] = pt;
45 pt.Set2( ip.x, -ip.y); arr[idx++] = pt;
46 pt.Set2( -ip.x, -ip.y); arr[idx++] = pt;
47 pt.Set2(1+ip.x, ip.y); arr[idx++] = pt;
48 pt.Set2( ip.x, 1+ip.y); arr[idx++] = pt;
49 pt.Set2(1+ip.x, 1+ip.y); arr[idx++] = pt;
50 break;
51 case 3:
52 arr.SetSize(15);
53
54 pt.Set3( ip.x, ip.y, ip.z ); arr[idx++] = pt;
55 pt.Set3( -ip.x, ip.y, ip.z ); arr[idx++] = pt;
56 pt.Set3( ip.x, -ip.y, ip.z ); arr[idx++] = pt;
57 pt.Set3( -ip.x, -ip.y, ip.z ); arr[idx++] = pt;
58 pt.Set3( ip.x, ip.y, -ip.z ); arr[idx++] = pt;
59 pt.Set3( -ip.x, ip.y, -ip.z ); arr[idx++] = pt;
60 pt.Set3( ip.x, -ip.y, -ip.z ); arr[idx++] = pt;
61 pt.Set3( -ip.x, -ip.y, -ip.z ); arr[idx++] = pt;
62 pt.Set3(1+ip.x, ip.y, ip.z ); arr[idx++] = pt;
63 pt.Set3( ip.x, 1+ip.y, ip.z ); arr[idx++] = pt;
64 pt.Set3(1+ip.x, 1+ip.y, ip.z ); arr[idx++] = pt;
65 pt.Set3( ip.x, ip.y, 1+ip.z ); arr[idx++] = pt;
66 pt.Set3(1+ip.x, ip.y, 1+ip.z ); arr[idx++] = pt;
67 pt.Set3( ip.x, 1+ip.y, 1+ip.z ); arr[idx++] = pt;
68 pt.Set3(1+ip.x, 1+ip.y, 1+ip.z ); arr[idx++] = pt;
69 break;
70 }
71 }
72
73 /**
74 * Tests fe->CalcShape() over a grid of IntegrationPoints
75 * of resolution res. Also tests at integration points
76 * that are outside the element.
77 */
TestCalcShape(FiniteElement * fe,int res,double tol=1e-12)78 void TestCalcShape(FiniteElement* fe, int res, double tol=1e-12)
79 {
80 int dim = fe->GetDim();
81
82 Vector weights( fe->GetDof() );
83
84 // Get a uniform grid or integration points
85 RefinedGeometry* ref = GlobGeometryRefiner.Refine( fe->GetGeomType(), res);
86 const IntegrationRule& intRule = ref->RefPts;
87
88 int npoints = intRule.GetNPoints();
89 for (int i=0; i < npoints; ++i)
90 {
91 // Get the current integration point from intRule
92 IntegrationPoint pt = intRule.IntPoint(i);
93
94 // Get several variants of this integration point
95 // some of which are inside the element and some are outside
96 Array<IntegrationPoint> ipArr;
97 GetRelatedIntegrationPoints( pt, dim, ipArr );
98
99 // For each such integration point check that the weights
100 // from CalcShape() sum to one
101 for (int j=0; j < ipArr.Size(); ++j)
102 {
103 IntegrationPoint& ip = ipArr[j];
104 fe->CalcShape(ip, weights);
105 REQUIRE(weights.Sum() == MFEM_Approx(1., tol, tol));
106 }
107 }
108 }
109
110
111 TEST_CASE("CalcShape for several Lagrange FiniteElement instances",
112 "[Lagrange1DFiniteElement]"
113 "[BiLinear2DFiniteElement]"
114 "[BiQuad2DFiniteElement]"
115 "[LagrangeHexFiniteElement]")
116 {
117 int maxOrder = 5;
118 int resolution = 10;
119
120 SECTION("Lagrange1DFiniteElement")
121 {
122 for (int order =1; order <= maxOrder; ++order)
123 {
124 std::cout << "Testing Lagrange1DFiniteElement::CalcShape() "
125 << "for order " << order << std::endl;
126 Lagrange1DFiniteElement fe(order);
127 TestCalcShape(&fe, resolution);
128 }
129 }
130
131 SECTION("BiLinear2DFiniteElement")
132 {
133 std::cout << "Testing BiLinear2DFiniteElement::CalcShape()" << std::endl;
134 BiLinear2DFiniteElement fe;
135 TestCalcShape(&fe, resolution);
136 }
137
138 SECTION("BiQuad2DFiniteElement")
139 {
140 std::cout << "Testing BiQuad2DFiniteElement::CalcShape()" << std::endl;
141 BiQuad2DFiniteElement fe;
142 TestCalcShape(&fe, resolution);
143 }
144
145
146 SECTION("LagrangeHexFiniteElement")
147 {
148 std::cout << "Testing LagrangeHexFiniteElement::CalcShape() "
149 << "for order 2" << std::endl;
150
151 // Comments for LagrangeHexFiniteElement state
152 // that only degree 2 is functional for this class
153 LagrangeHexFiniteElement fe(2);
154 TestCalcShape(&fe, resolution);
155 }
156 }
157
158 TEST_CASE("CalcShape for several H1 FiniteElement instances",
159 "[H1_SegmentElement]"
160 "[H1_TriangleElement]"
161 "[H1_QuadrilateralElement]"
162 "[H1_TetrahedronElement]"
163 "[H1_HexahedronElement]"
164 "[H1_WedgeElement]")
165 {
166 int maxOrder = 5;
167 int resolution = 10;
168
169 SECTION("H1_SegmentElement")
170 {
171 for (int order =1; order <= maxOrder; ++order)
172 {
173 std::cout << "Testing H1_SegmentElement::CalcShape() "
174 << "for order " << order << std::endl;
175 H1_SegmentElement fe(order);
176 TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
177 }
178 }
179
180 SECTION("H1_TriangleElement")
181 {
182 for (int order =1; order <= maxOrder; ++order)
183 {
184 std::cout << "Testing H1_TriangleElement::CalcShape() "
185 << "for order " << order << std::endl;
186 H1_TriangleElement fe(order);
187 TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
188 }
189 }
190
191 SECTION("H1_QuadrilateralElement")
192 {
193 for (int order =1; order <= maxOrder; ++order)
194 {
195 std::cout << "Testing H1_QuadrilateralElement::CalcShape() "
196 << "for order " << order << std::endl;
197 H1_QuadrilateralElement fe(order);
198 TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
199 }
200 }
201
202 SECTION("H1_TetrahedronElement")
203 {
204 for (int order =1; order <= maxOrder; ++order)
205 {
206 std::cout << "Testing H1_TetrahedronElement::CalcShape() "
207 << "for order " << order << std::endl;
208 H1_TetrahedronElement fe(order);
209 TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
210 }
211 }
212
213 SECTION("H1_HexahedronElement")
214 {
215 for (int order =1; order <= maxOrder; ++order)
216 {
217 std::cout << "Testing H1_HexahedronElement::CalcShape() "
218 << "for order " << order << std::endl;
219 H1_HexahedronElement fe(order);
220 TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
221 }
222 }
223
224 SECTION("H1_WedgeElement")
225 {
226 for (int order =1; order <= maxOrder; ++order)
227 {
228 std::cout << "Testing H1_WedgeElement::CalcShape() "
229 << "for order " << order << std::endl;
230 H1_WedgeElement fe(order);
231 TestCalcShape(&fe, resolution, 2e-11*std::pow(10, order));
232 }
233 }
234
235 }
236