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