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