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 namespace mfem
16 {
17
18 namespace build_dof_to_arrays
19 {
20
21 static double a_ = M_PI;
22 static double b_ = M_PI / sqrt(2.0);
23 static double c_ = M_PI / 2.0;
24
25 enum MeshType
26 {
27 SEGMENT = 0,
28 QUADRILATERAL = 1,
29 TRIANGLE2A = 2,
30 TRIANGLE2B = 3,
31 TRIANGLE2C = 4,
32 TRIANGLE4 = 5,
33 MIXED2D = 6,
34 HEXAHEDRON = 7,
35 HEXAHEDRON2A = 8,
36 HEXAHEDRON2B = 9,
37 HEXAHEDRON2C = 10,
38 HEXAHEDRON2D = 11,
39 WEDGE2 = 12,
40 TETRAHEDRA = 13,
41 WEDGE4 = 14,
42 MIXED3D6 = 15,
43 MIXED3D8 = 16
44 };
45
46 Mesh * GetMesh(MeshType type);
47
48 enum BasisType
49 {
50 H1 = 0, ND = 1, RT = 2, L2 = 3
51 };
52
53 TEST_CASE("Build Dof To Arrays",
54 "[BuildDofToArrays]"
55 "[FiniteElementSpace]")
56 {
57 int order = 3;
58
59 Array<int> dofs;
60
61 for (int mt = (int)MeshType::SEGMENT;
62 mt <= (int)MeshType::MIXED3D8; mt++)
63 {
64 Mesh *mesh = GetMesh((MeshType)mt);
65 int dim = mesh->Dimension();
66 if (dim < 3 ||
67 mt == MeshType::HEXAHEDRON ||
68 mt == MeshType::WEDGE2 ||
69 mt == MeshType::TETRAHEDRA ||
70 mt == MeshType::WEDGE4 ||
71 mt == MeshType::MIXED3D8 )
72 {
73 mesh->UniformRefinement();
74 }
75 if (dim == 3)
76 {
77 mesh->ReorientTetMesh();
78 }
79
80 for (int bt = (int)BasisType::H1; bt <= (int)BasisType::L2; bt++)
81 {
82 if (dim == 1 && bt == (int)BasisType::ND) { continue; }
83 if (dim == 1 && bt == (int)BasisType::RT) { continue; }
84 // FIX THIS: ND Wedge elements are not yet implemented
85 if ((mt == (int)MeshType::WEDGE2 ||
86 mt == (int)MeshType::MIXED3D6 ||
87 mt == (int)MeshType::MIXED3D8) &&
88 bt == (int)BasisType::ND) { continue; }
89
90 int num_elem_fails = 0;
91 int num_rang_fails = 0;
92 int num_ldof_fails = 0;
93
to_string(mt)94 SECTION("Mesh Type: " + std::to_string(mt) +
95 ", Basis Type: " + std::to_string(bt))
96 {
97 FiniteElementCollection * fec = NULL;
98 if (bt == (int)BasisType::H1)
99 {
100 fec = new H1_FECollection(order, dim);
101 }
102 else if (bt == (int)BasisType::ND)
103 {
104 fec = new ND_FECollection(order, dim);
105 }
106 else if (bt == (int)BasisType::RT)
107 {
108 fec = new RT_FECollection(order-1, dim);
109 }
110 else
111 {
112 fec = new L2_FECollection(order-1, dim);
113 }
114 FiniteElementSpace fespace(mesh, fec);
115 int size = fespace.GetTrueVSize();
116
117 fespace.BuildDofToArrays();
118
119 for (int i = 0; i<size; i++)
120 {
121 int e = fespace.GetElementForDof(i);
122 int l = fespace.GetLocalDofForDof(i);
123
124 if (e < 0 || e >= mesh->GetNE()) { num_elem_fails++; }
125
126 fespace.GetElementDofs(e, dofs);
127
128 if (l < 0 || l >= dofs.Size()) { num_rang_fails++; }
129
130 int ldof = (dofs[l] >= 0) ? dofs[l] : (-1 - dofs[l]);
131
132 if (i != ldof) { num_ldof_fails++; }
133 }
134 delete fec;
135 }
136 REQUIRE(num_elem_fails == 0);
137 REQUIRE(num_rang_fails == 0);
138 REQUIRE(num_ldof_fails == 0);
139 }
140 }
141 }
142
143 #ifdef MFEM_USE_MPI
144 #
145 TEST_CASE("Build Dof To Arrays (Parallel)",
146 "[BuildDofToArrays]"
147 "[ParFiniteElementSpace]"
148 "[Parallel]")
149 {
150 int num_procs;
151 MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
152
153 int my_rank;
154 MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
155
156 int order = 3;
157
158 Array<int> dofs;
159
160 for (int mt = (int)MeshType::SEGMENT;
161 mt <= (int)MeshType::MIXED3D8; mt++)
162 {
163 Mesh *mesh = GetMesh((MeshType)mt);
164 int dim = mesh->Dimension();
165 if (dim < 3 ||
166 mt == MeshType::HEXAHEDRON ||
167 mt == MeshType::WEDGE2 ||
168 mt == MeshType::TETRAHEDRA ||
169 mt == MeshType::WEDGE4 ||
170 mt == MeshType::MIXED3D8 )
171 {
172 mesh->UniformRefinement();
173 }
174 while (mesh->GetNE() < num_procs)
175 {
176 mesh->UniformRefinement();
177 }
178 ParMesh pmesh(MPI_COMM_WORLD, *mesh);
179 delete mesh;
180 if (dim == 3)
181 {
182 pmesh.ReorientTetMesh();
183 }
184
185 for (int bt = (int)BasisType::H1; bt <= (int)BasisType::L2; bt++)
186 {
187 if (dim == 1 && bt == (int)BasisType::ND) { continue; }
188 if (dim == 1 && bt == (int)BasisType::RT) { continue; }
189 // FIX THIS: ND Wedge elements are not yet implemented
190 if ((mt == (int)MeshType::WEDGE2 ||
191 mt == (int)MeshType::MIXED3D6 ||
192 mt == (int)MeshType::MIXED3D8) &&
193 bt == (int)BasisType::ND) { continue; }
194
195 int num_elem_fails = 0;
196 int num_rang_fails = 0;
197 int num_ldof_fails = 0;
198
to_string(mt)199 SECTION("Mesh Type: " + std::to_string(mt) +
200 ", Basis Type: " + std::to_string(bt))
201 {
202 FiniteElementCollection * fec = NULL;
203 if (bt == (int)BasisType::H1)
204 {
205 fec = new H1_FECollection(order, dim);
206 }
207 else if (bt == (int)BasisType::ND)
208 {
209 fec = new ND_FECollection(order, dim);
210 }
211 else if (bt == (int)BasisType::RT)
212 {
213 fec = new RT_FECollection(order-1, dim);
214 }
215 else
216 {
217 fec = new L2_FECollection(order-1, dim);
218 }
219 ParFiniteElementSpace fespace(&pmesh, fec);
220 HYPRE_Int size = fespace.GetTrueVSize();
221
222 fespace.BuildDofToArrays();
223
224 for (int i = 0; i<size; i++)
225 {
226 int e = fespace.GetElementForDof(i);
227 int l = fespace.GetLocalDofForDof(i);
228
229 if (e < 0 || e >= pmesh.GetNE()) { num_elem_fails++; }
230
231 fespace.GetElementDofs(e, dofs);
232
233 if (l < 0 || l >= dofs.Size()) { num_rang_fails++; }
234
235 int ldof = (dofs[l] >= 0) ? dofs[l] : (-1 - dofs[l]);
236
237 if (i != ldof) { num_ldof_fails++; }
238 }
239 delete fec;
240 }
241 REQUIRE(num_elem_fails == 0);
242 REQUIRE(num_rang_fails == 0);
243 REQUIRE(num_ldof_fails == 0);
244 }
245 }
246 }
247 #endif // MFEM_USE_MPI
248
GetMesh(MeshType type)249 Mesh * GetMesh(MeshType type)
250 {
251 Mesh * mesh = NULL;
252
253 switch (type)
254 {
255 case SEGMENT:
256 mesh = new Mesh(1, 2, 1);
257 mesh->AddVertex(0.0);
258 mesh->AddVertex(a_);
259
260 mesh->AddSegment(0, 1);
261
262 mesh->AddBdrPoint(0);
263 mesh->AddBdrPoint(1);
264 break;
265 case QUADRILATERAL:
266 mesh = new Mesh(2, 4, 1);
267 mesh->AddVertex(0.0, 0.0);
268 mesh->AddVertex(a_, 0.0);
269 mesh->AddVertex(a_, b_);
270 mesh->AddVertex(0.0, b_);
271
272 mesh->AddQuad(0, 1, 2, 3);
273 break;
274 case TRIANGLE2A:
275 mesh = new Mesh(2, 4, 2);
276 mesh->AddVertex(0.0, 0.0);
277 mesh->AddVertex(a_, 0.0);
278 mesh->AddVertex(a_, b_);
279 mesh->AddVertex(0.0, b_);
280
281 mesh->AddTriangle(0, 1, 2);
282 mesh->AddTriangle(2, 3, 0);
283 break;
284 case TRIANGLE2B:
285 mesh = new Mesh(2, 4, 2);
286 mesh->AddVertex(0.0, 0.0);
287 mesh->AddVertex(a_, 0.0);
288 mesh->AddVertex(a_, b_);
289 mesh->AddVertex(0.0, b_);
290
291 mesh->AddTriangle(1, 2, 0);
292 mesh->AddTriangle(3, 0, 2);
293 break;
294 case TRIANGLE2C:
295 mesh = new Mesh(2, 4, 2);
296 mesh->AddVertex(0.0, 0.0);
297 mesh->AddVertex(a_, 0.0);
298 mesh->AddVertex(a_, b_);
299 mesh->AddVertex(0.0, b_);
300
301 mesh->AddTriangle(2, 0, 1);
302 mesh->AddTriangle(0, 2, 3);
303 break;
304 case TRIANGLE4:
305 mesh = new Mesh(2, 5, 4);
306 mesh->AddVertex(0.0, 0.0);
307 mesh->AddVertex(a_, 0.0);
308 mesh->AddVertex(a_, b_);
309 mesh->AddVertex(0.0, b_);
310 mesh->AddVertex(0.5 * a_, 0.5 * b_);
311
312 mesh->AddTriangle(0, 1, 4);
313 mesh->AddTriangle(1, 2, 4);
314 mesh->AddTriangle(2, 3, 4);
315 mesh->AddTriangle(3, 0, 4);
316 break;
317 case MIXED2D:
318 mesh = new Mesh(2, 6, 4);
319 mesh->AddVertex(0.0, 0.0);
320 mesh->AddVertex(a_, 0.0);
321 mesh->AddVertex(a_, b_);
322 mesh->AddVertex(0.0, b_);
323 mesh->AddVertex(0.5 * b_, 0.5 * b_);
324 mesh->AddVertex(a_ - 0.5 * b_, 0.5 * b_);
325
326 mesh->AddQuad(0, 1, 5, 4);
327 mesh->AddTriangle(1, 2, 5);
328 mesh->AddQuad(2, 3, 4, 5);
329 mesh->AddTriangle(3, 0, 4);
330 break;
331 case HEXAHEDRON:
332 mesh = new Mesh(3, 8, 1);
333 mesh->AddVertex(0.0, 0.0, 0.0);
334 mesh->AddVertex(a_, 0.0, 0.0);
335 mesh->AddVertex(a_, b_, 0.0);
336 mesh->AddVertex(0.0, b_, 0.0);
337 mesh->AddVertex(0.0, 0.0, c_);
338 mesh->AddVertex(a_, 0.0, c_);
339 mesh->AddVertex(a_, b_, c_);
340 mesh->AddVertex(0.0, b_, c_);
341
342 mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
343 break;
344 case HEXAHEDRON2A:
345 case HEXAHEDRON2B:
346 case HEXAHEDRON2C:
347 case HEXAHEDRON2D:
348 mesh = new Mesh(3, 12, 2);
349 mesh->AddVertex(0.0, 0.0, 0.0);
350 mesh->AddVertex(0.5 * a_, 0.0, 0.0);
351 mesh->AddVertex(a_, 0.0, 0.0);
352 mesh->AddVertex(a_, b_, 0.0);
353 mesh->AddVertex(0.5 * a_, b_, 0.0);
354 mesh->AddVertex(0.0, b_, 0.0);
355 mesh->AddVertex(0.0, 0.0, c_);
356 mesh->AddVertex(0.5 * a_, 0.0, c_);
357 mesh->AddVertex(a_, 0.0, c_);
358 mesh->AddVertex(a_, b_, c_);
359 mesh->AddVertex(0.5 * a_, b_, c_);
360 mesh->AddVertex(0.0,b_, c_);
361
362 mesh->AddHex(0, 5, 11, 6, 1, 4, 10, 7);
363
364 switch (type)
365 {
366 case HEXAHEDRON2A: // Face Orientation 1
367 mesh->AddHex(4, 10, 7, 1, 3, 9, 8, 2);
368 break;
369 case HEXAHEDRON2B: // Face Orientation 3
370 mesh->AddHex(10, 7, 1, 4, 9, 8, 2, 3);
371 break;
372 case HEXAHEDRON2C: // Face Orientation 5
373 mesh->AddHex(7, 1, 4, 10, 8, 2, 3, 9);
374 break;
375 case HEXAHEDRON2D: // Face Orientation 7
376 mesh->AddHex(1, 4, 10, 7, 2, 3, 9, 8);
377 break;
378 default:
379 // Cannot happen
380 break;
381 }
382 break;
383 case WEDGE2:
384 mesh = new Mesh(3, 8, 2);
385 mesh->AddVertex(0.0, 0.0, 0.0);
386 mesh->AddVertex(a_, 0.0, 0.0);
387 mesh->AddVertex(a_, b_, 0.0);
388 mesh->AddVertex(0.0, b_, 0.0);
389 mesh->AddVertex(0.0, 0.0, c_);
390 mesh->AddVertex(a_, 0.0, c_);
391 mesh->AddVertex(a_, b_, c_);
392 mesh->AddVertex(0.0, b_, c_);
393
394 mesh->AddWedge(0, 1, 2, 4, 5, 6);
395 mesh->AddWedge(0, 2, 3, 4, 6, 7);
396 break;
397 case TETRAHEDRA:
398 mesh = new Mesh(3, 8, 5);
399 mesh->AddVertex(0.0, 0.0, 0.0);
400 mesh->AddVertex(a_, 0.0, 0.0);
401 mesh->AddVertex(a_, b_, 0.0);
402 mesh->AddVertex(0.0, b_, 0.0);
403 mesh->AddVertex(0.0, 0.0, c_);
404 mesh->AddVertex(a_, 0.0, c_);
405 mesh->AddVertex(a_, b_, c_);
406 mesh->AddVertex(0.0, b_, c_);
407
408 mesh->AddTet(0, 2, 7, 5);
409 mesh->AddTet(6, 7, 2, 5);
410 mesh->AddTet(4, 7, 5, 0);
411 mesh->AddTet(1, 0, 5, 2);
412 mesh->AddTet(3, 7, 0, 2);
413 break;
414 case WEDGE4:
415 mesh = new Mesh(3, 10, 4);
416 mesh->AddVertex(0.0, 0.0, 0.0);
417 mesh->AddVertex(a_, 0.0, 0.0);
418 mesh->AddVertex(a_, b_, 0.0);
419 mesh->AddVertex(0.0, b_, 0.0);
420 mesh->AddVertex(0.5 * a_, 0.5 * b_, 0.0);
421 mesh->AddVertex(0.0, 0.0, c_);
422 mesh->AddVertex(a_, 0.0, c_);
423 mesh->AddVertex(a_, b_, c_);
424 mesh->AddVertex(0.0, b_, c_);
425 mesh->AddVertex(0.5 * a_, 0.5 * b_, c_);
426
427 mesh->AddWedge(0, 1, 4, 5, 6, 9);
428 mesh->AddWedge(1, 2, 4, 6, 7, 9);
429 mesh->AddWedge(2, 3, 4, 7, 8, 9);
430 mesh->AddWedge(3, 0, 4, 8, 5, 9);
431 break;
432 case MIXED3D6:
433 mesh = new Mesh(3, 12, 6);
434 mesh->AddVertex(0.0, 0.0, 0.0);
435 mesh->AddVertex(a_, 0.0, 0.0);
436 mesh->AddVertex(a_, b_, 0.0);
437 mesh->AddVertex(0.0, b_, 0.0);
438 mesh->AddVertex(0.5 * c_, 0.5 * c_, 0.5 * c_);
439 mesh->AddVertex(a_ - 0.5 * c_, 0.5 * c_, 0.5 * c_);
440 mesh->AddVertex(a_ - 0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
441 mesh->AddVertex(0.5 * c_, b_ - 0.5 * c_, 0.5 * c_);
442 mesh->AddVertex(0.0, 0.0, c_);
443 mesh->AddVertex(a_, 0.0, c_);
444 mesh->AddVertex(a_, b_, c_);
445 mesh->AddVertex(0.0, b_, c_);
446
447 mesh->AddHex(0, 1, 2, 3, 4, 5, 6, 7);
448 mesh->AddWedge(0, 4, 8, 1, 5, 9);
449 mesh->AddWedge(1, 5, 9, 2, 6, 10);
450 mesh->AddWedge(2, 6, 10, 3, 7, 11);
451 mesh->AddWedge(3, 7, 11, 0, 4, 8);
452 mesh->AddHex(4, 5, 6, 7, 8, 9, 10, 11);
453 break;
454 case MIXED3D8:
455 mesh = new Mesh(3, 10, 8);
456 mesh->AddVertex(0.0, 0.0, 0.0);
457 mesh->AddVertex(a_, 0.0, 0.0);
458 mesh->AddVertex(a_, b_, 0.0);
459 mesh->AddVertex(0.0, b_, 0.0);
460 mesh->AddVertex(0.25 * a_, 0.5 * b_, 0.5 * c_);
461 mesh->AddVertex(0.75 * a_, 0.5 * b_, 0.5 * c_);
462 mesh->AddVertex(0.0, 0.0, c_);
463 mesh->AddVertex(a_, 0.0, c_);
464 mesh->AddVertex(a_, b_, c_);
465 mesh->AddVertex(0.0, b_, c_);
466
467 mesh->AddWedge(0, 3, 4, 1, 2, 5);
468 mesh->AddWedge(3, 9, 4, 2, 8, 5);
469 mesh->AddWedge(9, 6, 4, 8, 7, 5);
470 mesh->AddWedge(6, 0, 4, 7, 1, 5);
471 mesh->AddTet(0, 3, 9, 4);
472 mesh->AddTet(0, 9, 6, 4);
473 mesh->AddTet(1, 7, 2, 5);
474 mesh->AddTet(8, 2, 7, 5);
475 break;
476 }
477 mesh->FinalizeTopology();
478
479 return mesh;
480 }
481
482 } // namespace build_dof_to_arrays
483
484 } // namespace mfem
485