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 <memory>
16 #include <array>
17
18 using namespace mfem;
19
20 #if defined(MFEM_USE_MPI)
21
22 namespace testhelper
23 {
SmoothSolutionX(const mfem::Vector & x)24 double SmoothSolutionX(const mfem::Vector& x)
25 {
26 return x(0);
27 }
28
SmoothSolutionY(const mfem::Vector & x)29 double SmoothSolutionY(const mfem::Vector& x)
30 {
31 return x(1);
32 }
33
SmoothSolutionZ(const mfem::Vector & x)34 double SmoothSolutionZ(const mfem::Vector& x)
35 {
36 return x(2);
37 }
38
NonsmoothSolutionX(const mfem::Vector & x)39 double NonsmoothSolutionX(const mfem::Vector& x)
40 {
41 return std::abs(x(0)-0.5);
42 }
43
NonsmoothSolutionY(const mfem::Vector & x)44 double NonsmoothSolutionY(const mfem::Vector& x)
45 {
46 return std::abs(x(1)-0.5);
47 }
48
NonsmoothSolutionZ(const mfem::Vector & x)49 double NonsmoothSolutionZ(const mfem::Vector& x)
50 {
51 return std::abs(x(2)-0.5);
52 }
53 }
54
55 TEST_CASE("Kelly Error Estimator on 2D NCMesh",
56 "[NCMesh], [Parallel]")
57 {
58 // Setup
59 const auto order = GENERATE(1, 3, 5);
60 Mesh mesh = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL);
61
62 // Make the mesh NC
63 mesh.EnsureNCMesh();
64 {
65 Array<int> elements_to_refine(1);
66 elements_to_refine[0] = 1;
67 mesh.GeneralRefinement(elements_to_refine, 1, 0);
68 }
69
70 auto pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
71 mesh.Clear();
72
73 H1_FECollection fe_coll(order, pmesh->Dimension());
74 ParFiniteElementSpace fespace(pmesh, &fe_coll);
75
76 SECTION("Perfect Approximation X")
77 {
78 FunctionCoefficient u_analytic(testhelper::SmoothSolutionX);
79 ParGridFunction u_gf(&fespace);
80 u_gf.ProjectCoefficient(u_analytic);
81
82 L2_FECollection flux_fec(order, pmesh->Dimension());
83 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
84 DiffusionIntegrator di;
85 KellyErrorEstimator estimator(di, u_gf, flux_fes);
86
87 auto &local_errors = estimator.GetLocalErrors();
88 for (int i=0; i<local_errors.Size(); i++)
89 {
90 REQUIRE(local_errors(i) == MFEM_Approx(0.0));
91 }
92 REQUIRE(estimator.GetTotalError() == MFEM_Approx(0.0));
93 }
94
95 SECTION("Perfect Approximation Y")
96 {
97 FunctionCoefficient u_analytic(testhelper::SmoothSolutionY);
98 ParGridFunction u_gf(&fespace);
99 u_gf.ProjectCoefficient(u_analytic);
100
101 L2_FECollection flux_fec(order, pmesh->Dimension());
102 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
103 DiffusionIntegrator di;
104 KellyErrorEstimator estimator(di, u_gf, flux_fes);
105
106 auto &local_errors = estimator.GetLocalErrors();
107 for (int i=0; i<local_errors.Size(); i++)
108 {
109 REQUIRE(local_errors(i) == MFEM_Approx(0.0));
110 }
111 REQUIRE(estimator.GetTotalError() == MFEM_Approx(0.0));
112 }
113
114 SECTION("Nonsmooth Approximation X")
115 {
116 FunctionCoefficient u_analytic(testhelper::NonsmoothSolutionX);
117 ParGridFunction u_gf(&fespace);
118 u_gf.ProjectCoefficient(u_analytic);
119
120 L2_FECollection flux_fec(order, pmesh->Dimension());
121 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
122 DiffusionIntegrator di;
123 KellyErrorEstimator estimator(di, u_gf, flux_fes);
124
125 auto &local_errors = estimator.GetLocalErrors();
126 for (int i=0; i<local_errors.Size(); i++)
127 {
128 REQUIRE(local_errors(i) >= 0.0);
129 }
130 REQUIRE(estimator.GetTotalError() > 0.0);
131 }
132
133 SECTION("Nonsmooth Approximation Y")
134 {
135 FunctionCoefficient u_analytic(testhelper::NonsmoothSolutionY);
136 ParGridFunction u_gf(&fespace);
137 u_gf.ProjectCoefficient(u_analytic);
138
139 L2_FECollection flux_fec(order, pmesh->Dimension());
140 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
141 DiffusionIntegrator di;
142 KellyErrorEstimator estimator(di, u_gf, flux_fes);
143
144 auto &local_errors = estimator.GetLocalErrors();
145 for (int i=0; i<local_errors.Size(); i++)
146 {
147 REQUIRE(local_errors(i) >= MFEM_Approx(0.0));
148 }
149 REQUIRE(estimator.GetTotalError() > 0.0);
150 }
151
152 delete pmesh;
153 }
154
155 TEST_CASE("Kelly Error Estimator on 2D NCMesh embedded in 3D",
156 "[NCMesh], [Parallel]")
157 {
158 // Setup
159 const auto order = GENERATE(1, 3, 5);
160
161 // Manually construct embedded mesh
162 std::array<double, 4*3> vertices =
163 {
164 0.0,0.0,0.0,
165 0.0,1.0,0.0,
166 1.0,1.0,0.0,
167 1.0,0.0,0.0
168 };
169
170 std::array<int, 4> element_indices =
171 {
172 0,1,2,3
173 };
174
175 std::array<int, 1> element_attributes =
176 {
177 1
178 };
179
180 std::array<int, 8> boundary_indices =
181 {
182 0,1,
183 1,2,
184 2,3,
185 3,0
186 };
187
188 std::array<int, 4> boundary_attributes =
189 {
190 1,
191 1,
192 1,
193 1
194 };
195
196 auto mesh = new Mesh(
197 vertices.data(), 4,
198 element_indices.data(), Geometry::SQUARE,
199 element_attributes.data(), 1,
200 boundary_indices.data(), Geometry::SEGMENT,
201 boundary_attributes.data(), 4,
202 2, 3
203 );
204 mesh->UniformRefinement();
205 mesh->Finalize();
206
207 // Make the mesh NC
208 mesh->EnsureNCMesh();
209 {
210 Array<int> elements_to_refine(1);
211 elements_to_refine[0] = 1;
212 mesh->GeneralRefinement(elements_to_refine, 1, 0);
213 }
214
215 auto pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
216 delete mesh;
217
218 H1_FECollection fe_coll(order, pmesh->Dimension());
219 ParFiniteElementSpace fespace(pmesh, &fe_coll);
220
221 SECTION("Perfect Approximation X")
222 {
223 FunctionCoefficient u_analytic(testhelper::SmoothSolutionX);
224 ParGridFunction u_gf(&fespace);
225 u_gf.ProjectCoefficient(u_analytic);
226
227 L2_FECollection flux_fec(order, pmesh->Dimension());
228 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
229 DiffusionIntegrator di;
230 KellyErrorEstimator estimator(di, u_gf, flux_fes);
231
232 auto &local_errors = estimator.GetLocalErrors();
233 for (int i=0; i<local_errors.Size(); i++)
234 {
235 REQUIRE(local_errors(i) == MFEM_Approx(0.0));
236 }
237 REQUIRE(estimator.GetTotalError() == MFEM_Approx(0.0));
238 }
239
240 SECTION("Perfect Approximation Y")
241 {
242 FunctionCoefficient u_analytic(testhelper::SmoothSolutionY);
243 ParGridFunction u_gf(&fespace);
244 u_gf.ProjectCoefficient(u_analytic);
245
246 L2_FECollection flux_fec(order, pmesh->Dimension());
247 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
248 DiffusionIntegrator di;
249 KellyErrorEstimator estimator(di, u_gf, flux_fes);
250
251 auto &local_errors = estimator.GetLocalErrors();
252 for (int i=0; i<local_errors.Size(); i++)
253 {
254 REQUIRE(local_errors(i) == MFEM_Approx(0.0));
255 }
256 REQUIRE(estimator.GetTotalError() == MFEM_Approx(0.0));
257 }
258
259 SECTION("Nonsmooth Approximation X")
260 {
261 FunctionCoefficient u_analytic(testhelper::NonsmoothSolutionX);
262 ParGridFunction u_gf(&fespace);
263 u_gf.ProjectCoefficient(u_analytic);
264
265 L2_FECollection flux_fec(order, pmesh->Dimension());
266 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
267 DiffusionIntegrator di;
268 KellyErrorEstimator estimator(di, u_gf, flux_fes);
269
270 auto &local_errors = estimator.GetLocalErrors();
271 for (int i=0; i<local_errors.Size(); i++)
272 {
273 REQUIRE(local_errors(i) >= 0.0);
274 }
275 REQUIRE(estimator.GetTotalError() > 0.0);
276 }
277
278 SECTION("Nonsmooth Approximation Y")
279 {
280 FunctionCoefficient u_analytic(testhelper::NonsmoothSolutionY);
281 ParGridFunction u_gf(&fespace);
282 u_gf.ProjectCoefficient(u_analytic);
283
284 L2_FECollection flux_fec(order, pmesh->Dimension());
285 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
286 DiffusionIntegrator di;
287 KellyErrorEstimator estimator(di, u_gf, flux_fes);
288
289 auto &local_errors = estimator.GetLocalErrors();
290 for (int i=0; i<local_errors.Size(); i++)
291 {
292 REQUIRE(local_errors(i) >= MFEM_Approx(0.0));
293 }
294 REQUIRE(estimator.GetTotalError() > 0.0);
295 }
296
297 delete pmesh;
298 }
299
300 TEST_CASE("Kelly Error Estimator on 3D NCMesh",
301 "[NCMesh], [Parallel]")
302 {
303 // Setup
304 const auto order = GENERATE(1, 3, 5);
305 Mesh mesh = Mesh::MakeCartesian3D(2, 2, 2, Element::HEXAHEDRON);
306
307 // Make the mesh NC
308 mesh.EnsureNCMesh();
309 {
310 Array<int> elements_to_refine(1);
311 elements_to_refine[0] = 1;
312 mesh.GeneralRefinement(elements_to_refine, 1, 0);
313 }
314
315 auto pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
316 mesh.Clear();
317
318 H1_FECollection fe_coll(order, pmesh->Dimension());
319 ParFiniteElementSpace fespace(pmesh, &fe_coll);
320
321 SECTION("Perfect Approximation X")
322 {
323 FunctionCoefficient u_analytic(testhelper::SmoothSolutionX);
324 ParGridFunction u_gf(&fespace);
325 u_gf.ProjectCoefficient(u_analytic);
326
327 L2_FECollection flux_fec(order, pmesh->Dimension());
328 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
329 DiffusionIntegrator di;
330 KellyErrorEstimator estimator(di, u_gf, flux_fes);
331
332 auto &local_errors = estimator.GetLocalErrors();
333 for (int i=0; i<local_errors.Size(); i++)
334 {
335 REQUIRE(local_errors(i) == MFEM_Approx(0.0));
336 }
337 REQUIRE(estimator.GetTotalError() == MFEM_Approx(0.0));
338 }
339
340 SECTION("Perfect Approximation Y")
341 {
342 FunctionCoefficient u_analytic(testhelper::SmoothSolutionY);
343 ParGridFunction u_gf(&fespace);
344 u_gf.ProjectCoefficient(u_analytic);
345
346 L2_FECollection flux_fec(order, pmesh->Dimension());
347 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
348 DiffusionIntegrator di;
349 KellyErrorEstimator estimator(di, u_gf, flux_fes);
350
351 auto &local_errors = estimator.GetLocalErrors();
352 for (int i=0; i<local_errors.Size(); i++)
353 {
354 REQUIRE(local_errors(i) == MFEM_Approx(0.0));
355 }
356 REQUIRE(estimator.GetTotalError() == MFEM_Approx(0.0));
357 }
358
359 SECTION("Perfect Approximation Z")
360 {
361 FunctionCoefficient u_analytic(testhelper::SmoothSolutionZ);
362 ParGridFunction u_gf(&fespace);
363 u_gf.ProjectCoefficient(u_analytic);
364
365 L2_FECollection flux_fec(order, pmesh->Dimension());
366 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
367 DiffusionIntegrator di;
368 KellyErrorEstimator estimator(di, u_gf, flux_fes);
369
370 auto &local_errors = estimator.GetLocalErrors();
371 for (int i=0; i<local_errors.Size(); i++)
372 {
373 REQUIRE(local_errors(i) == MFEM_Approx(0.0));
374 }
375 REQUIRE(estimator.GetTotalError() == MFEM_Approx(0.0));
376 }
377
378 SECTION("Nonsmooth Approximation X")
379 {
380 FunctionCoefficient u_analytic(testhelper::NonsmoothSolutionX);
381 ParGridFunction u_gf(&fespace);
382 u_gf.ProjectCoefficient(u_analytic);
383
384 L2_FECollection flux_fec(order, pmesh->Dimension());
385 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
386 DiffusionIntegrator di;
387 KellyErrorEstimator estimator(di, u_gf, flux_fes);
388
389 auto &local_errors = estimator.GetLocalErrors();
390 for (int i=0; i<local_errors.Size(); i++)
391 {
392 REQUIRE(local_errors(i) >= 0.0);
393 }
394 REQUIRE(estimator.GetTotalError() > 0.0);
395 }
396
397 SECTION("Nonsmooth Approximation Y")
398 {
399 FunctionCoefficient u_analytic(testhelper::NonsmoothSolutionY);
400 ParGridFunction u_gf(&fespace);
401 u_gf.ProjectCoefficient(u_analytic);
402
403 L2_FECollection flux_fec(order, pmesh->Dimension());
404 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
405 DiffusionIntegrator di;
406 KellyErrorEstimator estimator(di, u_gf, flux_fes);
407
408 auto &local_errors = estimator.GetLocalErrors();
409 for (int i=0; i<local_errors.Size(); i++)
410 {
411 REQUIRE(local_errors(i) >= MFEM_Approx(0.0));
412 }
413 REQUIRE(estimator.GetTotalError() > 0.0);
414 }
415
416 SECTION("Nonsmooth Approximation Z")
417 {
418 FunctionCoefficient u_analytic(testhelper::NonsmoothSolutionZ);
419 ParGridFunction u_gf(&fespace);
420 u_gf.ProjectCoefficient(u_analytic);
421
422 L2_FECollection flux_fec(order, pmesh->Dimension());
423 ParFiniteElementSpace flux_fes(pmesh, &flux_fec, pmesh->SpaceDimension());
424 DiffusionIntegrator di;
425 KellyErrorEstimator estimator(di, u_gf, flux_fes);
426
427 auto &local_errors = estimator.GetLocalErrors();
428 for (int i=0; i<local_errors.Size(); i++)
429 {
430 REQUIRE(local_errors(i) >= 0.0);
431 }
432 REQUIRE(estimator.GetTotalError() > 0.0);
433 }
434
435 delete pmesh;
436 }
437
438 #endif
439