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 "unit_tests.hpp"
13 #include "mfem.hpp"
14
15 using namespace mfem;
16
17 namespace derefine
18 {
19
20 int dimension;
coeff(const Vector & x)21 double coeff(const Vector& x)
22 {
23 if (dimension == 2)
24 {
25 return sin(10.0*(x[0]+x[1]));
26 }
27 else
28 {
29 return sin(10.0*(x[0]+x[1]+x[2]));
30 }
31 }
32
33 TEST_CASE("Derefine")
34 {
35 for (dimension = 2; dimension <= 3; ++dimension)
36 {
37 for (int order = 0; order <= 2; ++order)
38 {
39 for (int map_type = FiniteElement::VALUE; map_type <= FiniteElement::INTEGRAL;
40 ++map_type)
41 {
42 const int ne = 8;
43 Mesh mesh;
44 if (dimension == 2)
45 {
46 mesh = Mesh::MakeCartesian2D(
47 ne, ne, Element::QUADRILATERAL, true, 1.0, 1.0);
48 }
49 else
50 {
51 mesh = Mesh::MakeCartesian3D(
52 ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
53 }
54
55 mesh.EnsureNCMesh();
56 mesh.SetCurvature(std::max(order,1), false, dimension, Ordering::byNODES);
57
58 L2_FECollection fec(order, dimension, BasisType::Positive, map_type);
59
60 FiniteElementSpace fespace(&mesh, &fec);
61
62 GridFunction x(&fespace);
63
64 FunctionCoefficient c(coeff);
65 x.ProjectCoefficient(c);
66
67 fespace.Update();
68 x.Update();
69
70 Array<Refinement> refinements;
71 refinements.Append(Refinement(1));
72 refinements.Append(Refinement(2));
73
74 int nonconformity_limit = 0; // 0 meaning allow unlimited ratio
75
76 // First refine two elements.
77 mesh.GeneralRefinement(refinements, 1, nonconformity_limit);
78
79 fespace.Update();
80 x.Update();
81
82 // Now refine one more element and then derefine it, comparing x before and after.
83 Vector diff(x);
84
85 refinements.DeleteAll();
86 refinements.Append(Refinement(2));
87 mesh.GeneralRefinement(refinements, 1, nonconformity_limit);
88
89 fespace.Update();
90 x.Update();
91
92 // Derefine by setting 0 error on the fine elements in coarse element 2.
93 Table coarse_to_fine_;
94 Table ref_type_to_matrix;
95 Array<int> coarse_to_ref_type;
96 Array<Geometry::Type> ref_type_to_geom;
97 const CoarseFineTransformations &rtrans = mesh.GetRefinementTransforms();
98 rtrans.GetCoarseToFineMap(mesh, coarse_to_fine_, coarse_to_ref_type,
99 ref_type_to_matrix, ref_type_to_geom);
100 Array<int> tabrow;
101
102 Vector local_err(mesh.GetNE());
103 double threshold = 1.0;
104 local_err = 2*threshold;
105 coarse_to_fine_.GetRow(2, tabrow);
106 for (int j = 0; j < tabrow.Size(); j++) { local_err(tabrow[j]) = 0.0; }
107 mesh.DerefineByError(local_err, threshold, 0, 1);
108
109 fespace.Update();
110 x.Update();
111
112 diff -= x;
113 REQUIRE(diff.Norml2() / x.Norml2() < 1e-11);
114 }
115 }
116 }
117 }
118
119 #ifdef MFEM_USE_MPI
120 TEST_CASE("ParDerefine", "[Parallel]")
121 {
122 for (dimension = 2; dimension <= 3; ++dimension)
123 {
124 for (int order = 0; order <= 2; ++order)
125 {
126 for (int map_type = FiniteElement::VALUE; map_type <= FiniteElement::INTEGRAL;
127 ++map_type)
128 {
129 const int ne = 8;
130 Mesh mesh;
131 if (dimension == 2)
132 {
133 mesh = Mesh::MakeCartesian2D(
134 ne, ne, Element::QUADRILATERAL, true, 1.0, 1.0);
135 }
136 else
137 {
138 mesh = Mesh::MakeCartesian3D(
139 ne, ne, ne, Element::HEXAHEDRON, 1.0, 1.0, 1.0);
140 }
141
142 mesh.EnsureNCMesh();
143 mesh.SetCurvature(std::max(order,1), false, dimension, Ordering::byNODES);
144
145 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
146
147 L2_FECollection fec(order, dimension, BasisType::Positive, map_type);
148 ParFiniteElementSpace fespace(pmesh, &fec);
149
150 ParGridFunction x(&fespace);
151
152 FunctionCoefficient c(coeff);
153 x.ProjectCoefficient(c);
154
155 fespace.Update();
156 x.Update();
157
158 // Refine two elements on each process and then derefine, comparing
159 // x before and after.
160 Vector diff(x);
161
162 Array<Refinement> refinements;
163 refinements.Append(Refinement(1));
164 refinements.Append(Refinement(2));
165
166 int nonconformity_limit = 0; // 0 meaning allow unlimited ratio
167
168 pmesh->GeneralRefinement(refinements, 1, nonconformity_limit);
169
170 fespace.Update();
171 x.Update();
172
173 // Derefine by setting 0 error on all fine elements.
174 Vector local_err(pmesh->GetNE());
175 double threshold = 1.0;
176 local_err = 0.0;
177 pmesh->DerefineByError(local_err, threshold, 0, 1);
178
179 fespace.Update();
180 x.Update();
181
182 diff -= x;
183 REQUIRE(diff.Norml2() / x.Norml2() < 1e-11);
184
185 delete pmesh;
186 }
187 }
188 }
189 }
190 #endif
191 }
192