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 // Implementation of class Tetrahedron
13 
14 #include "mesh_headers.hpp"
15 
16 namespace mfem
17 {
18 
Tetrahedron(const int * ind,int attr)19 Tetrahedron::Tetrahedron(const int *ind, int attr)
20    : Element(Geometry::TETRAHEDRON)
21 {
22    attribute = attr;
23    for (int i = 0; i < 4; i++)
24    {
25       indices[i] = ind[i];
26    }
27    refinement_flag = 0;
28    transform = 0;
29 }
30 
Tetrahedron(int ind1,int ind2,int ind3,int ind4,int attr)31 Tetrahedron::Tetrahedron(int ind1, int ind2, int ind3, int ind4, int attr)
32    : Element(Geometry::TETRAHEDRON)
33 {
34    attribute  = attr;
35    indices[0] = ind1;
36    indices[1] = ind2;
37    indices[2] = ind3;
38    indices[3] = ind4;
39    refinement_flag = 0;
40    transform = 0;
41 }
42 
Init(int ind1,int ind2,int ind3,int ind4,int attr,int ref_flag)43 void Tetrahedron::Init(int ind1, int ind2, int ind3, int ind4, int attr,
44                        int ref_flag)
45 {
46    attribute  = attr;
47    indices[0] = ind1;
48    indices[1] = ind2;
49    indices[2] = ind3;
50    indices[3] = ind4;
51    refinement_flag = ref_flag;
52    transform = 0;
53 }
54 
ParseRefinementFlag(int refinement_edges[2],int & type,int & flag)55 void Tetrahedron::ParseRefinementFlag(int refinement_edges[2], int &type,
56                                       int &flag)
57 {
58    int i, f = refinement_flag;
59 
60    MFEM_VERIFY(f != 0, "tetrahedron is not marked");
61 
62    for (i = 0; i < 2; i++)
63    {
64       refinement_edges[i] = f & 7;
65       f = f >> 3;
66    }
67    type = f & 7;
68    flag = (f >> 3);
69 }
70 
CreateRefinementFlag(int refinement_edges[2],int type,int flag)71 void Tetrahedron::CreateRefinementFlag(int refinement_edges[2], int type,
72                                        int flag)
73 {
74    // Check for correct type
75 #ifdef MFEM_DEBUG
76    int e1, e2;
77    e1 = refinement_edges[0];
78    e2 = refinement_edges[1];
79    // if (e1 > e2)  e1 = e2, e2 = refinement_edges[0];
80    switch (type)
81    {
82       case Tetrahedron::TYPE_PU:
83          if (e1 == 2 && e2 == 1) { break; }
84          // if (e1 == 3 && e2 == 4) break;
85          mfem_error("Error in Tetrahedron::CreateRefinementFlag(...) #1");
86          break;
87       case Tetrahedron::TYPE_A:
88          if (e1 == 3 && e2 == 1) { break; }
89          if (e1 == 2 && e2 == 4) { break; }
90          // if (flag == 0)  // flag is assumed to be the generation
91          //   if (e2 == 5)
92          //      if (e1 >= 1 && e1 <= 5) break; // type is actually O or M
93          //                                     //   ==>  ok for generation = 0
94          mfem_error("Error in Tetrahedron::CreateRefinementFlag(...) #2");
95          break;
96       case Tetrahedron::TYPE_PF:
97          if (flag > 0)  // PF is ok only for generation > 0
98          {
99             if (e1 == 2 && e2 == 1) { break; }
100             // if (e1 == 3 && e2 == 4) break;
101          }
102          mfem_error("Error in Tetrahedron::CreateRefinementFlag(...) #3");
103          break;
104       case Tetrahedron::TYPE_O:
105          if (flag == 0 && e1 == 5 && e2 == 5)
106          {
107             break;
108          }
109          mfem_error("Error in Tetrahedron::CreateRefinementFlag(...) #4");
110          break;
111       case Tetrahedron::TYPE_M:
112          if (flag == 0)
113          {
114             if (e1 == 5 && e2 == 1) { break; }
115             if (e1 == 2 && e2 == 5) { break; }
116          }
117          mfem_error("Error in Tetrahedron::CreateRefinementFlag(...) #5");
118          break;
119       default:
120          mfem_error("Error in Tetrahedron::CreateRefinementFlag(...) #6");
121          break;
122    }
123 #endif
124 
125    refinement_flag = flag;
126    refinement_flag <<= 3;
127 
128    refinement_flag |= type;
129    refinement_flag <<= 3;
130 
131    refinement_flag |= refinement_edges[1];
132    refinement_flag <<= 3;
133 
134    refinement_flag |= refinement_edges[0];
135 }
136 
GetMarkedFace(const int face,int * fv)137 void Tetrahedron::GetMarkedFace(const int face, int *fv)
138 {
139    int re[2], type, flag, *tv = this->indices;
140    ParseRefinementFlag(re, type, flag);
141    switch (face)
142    {
143       case 0:
144          switch (re[1])
145          {
146             case 1: fv[0] = tv[1]; fv[1] = tv[2]; fv[2] = tv[3]; break;
147             case 4: fv[0] = tv[3]; fv[1] = tv[1]; fv[2] = tv[2]; break;
148             case 5: fv[0] = tv[2]; fv[1] = tv[3]; fv[2] = tv[1]; break;
149          }
150          break;
151       case 1:
152          switch (re[0])
153          {
154             case 2: fv[0] = tv[2]; fv[1] = tv[0]; fv[2] = tv[3]; break;
155             case 3: fv[0] = tv[0]; fv[1] = tv[3]; fv[2] = tv[2]; break;
156             case 5: fv[0] = tv[3]; fv[1] = tv[2]; fv[2] = tv[0]; break;
157          }
158          break;
159       case 2:
160          fv[0] = tv[0]; fv[1] = tv[1]; fv[2] = tv[3];
161          break;
162       case 3:
163          fv[0] = tv[1]; fv[1] = tv[0]; fv[2] = tv[2];
164          break;
165    }
166 }
167 
NeedRefinement(HashTable<Hashed2> & v_to_v) const168 int Tetrahedron::NeedRefinement(HashTable<Hashed2> &v_to_v) const
169 {
170    if (v_to_v.FindId(indices[0], indices[1]) != -1) { return 1; }
171    if (v_to_v.FindId(indices[1], indices[2]) != -1) { return 1; }
172    if (v_to_v.FindId(indices[2], indices[0]) != -1) { return 1; }
173    if (v_to_v.FindId(indices[0], indices[3]) != -1) { return 1; }
174    if (v_to_v.FindId(indices[1], indices[3]) != -1) { return 1; }
175    if (v_to_v.FindId(indices[2], indices[3]) != -1) { return 1; }
176    return 0;
177 }
178 
SetVertices(const int * ind)179 void Tetrahedron::SetVertices(const int *ind)
180 {
181    for (int i = 0; i < 4; i++)
182    {
183       indices[i] = ind[i];
184    }
185 }
186 
MarkEdge(const DSTable & v_to_v,const int * length)187 void Tetrahedron::MarkEdge(const DSTable &v_to_v, const int *length)
188 {
189    int ind[4], i, j, l, L, type;
190 
191    // determine the longest edge
192    L = length[v_to_v(indices[0], indices[1])]; j = 0;
193    if ((l = length[v_to_v(indices[1], indices[2])]) > L) { L = l; j = 1; }
194    if ((l = length[v_to_v(indices[2], indices[0])]) > L) { L = l; j = 2; }
195    if ((l = length[v_to_v(indices[0], indices[3])]) > L) { L = l; j = 3; }
196    if ((l = length[v_to_v(indices[1], indices[3])]) > L) { L = l; j = 4; }
197    if ((l = length[v_to_v(indices[2], indices[3])]) > L) { j = 5; }
198 
199    for (i = 0; i < 4; i++)
200    {
201       ind[i] = indices[i];
202    }
203 
204    switch (j)
205    {
206       case 1:
207          indices[0] = ind[1]; indices[1] = ind[2];
208          indices[2] = ind[0]; indices[3] = ind[3];
209          break;
210       case 2:
211          indices[0] = ind[2]; indices[1] = ind[0];
212          indices[2] = ind[1]; indices[3] = ind[3];
213          break;
214       case 3:
215          indices[0] = ind[3]; indices[1] = ind[0];
216          indices[2] = ind[2]; indices[3] = ind[1];
217          break;
218       case 4:
219          indices[0] = ind[1]; indices[1] = ind[3];
220          indices[2] = ind[2]; indices[3] = ind[0];
221          break;
222       case 5:
223          indices[0] = ind[2]; indices[1] = ind[3];
224          indices[2] = ind[0]; indices[3] = ind[1];
225          break;
226    }
227 
228    // Determine the two longest edges for the other two faces and
229    // store them in ind[0] and ind[1]
230    ind[0] = 2; ind[1] = 1;
231    L = length[v_to_v(indices[0], indices[2])];
232    if ((l = length[v_to_v(indices[0], indices[3])]) > L) { L = l; ind[0] = 3; }
233    if ((l = length[v_to_v(indices[2], indices[3])]) > L) { ind[0] = 5; }
234 
235    L = length[v_to_v(indices[1], indices[2])];
236    if ((l = length[v_to_v(indices[1], indices[3])]) > L) { L = l; ind[1] = 4; }
237    if ((l = length[v_to_v(indices[2], indices[3])]) > L) { ind[1] = 5; }
238 
239    j = 0;
240    switch (ind[0])
241    {
242       case 2:
243          switch (ind[1])
244          {
245             case 1:  type = Tetrahedron::TYPE_PU; break;
246             case 4:  type = Tetrahedron::TYPE_A;  break;
247             case 5:
248             default: type = Tetrahedron::TYPE_M;
249          }
250          break;
251       case 3:
252          switch (ind[1])
253          {
254             case 1:  type = Tetrahedron::TYPE_A;  break;
255             case 4:  type = Tetrahedron::TYPE_PU;
256                j = 1; ind[0] = 2; ind[1] = 1; break;
257             case 5:
258             default: type = Tetrahedron::TYPE_M;
259                j = 1; ind[0] = 5; ind[1] = 1;
260          }
261          break;
262       case 5:
263       default:
264          switch (ind[1])
265          {
266             case 1:  type = Tetrahedron::TYPE_M;  break;
267             case 4:  type = Tetrahedron::TYPE_M;
268                j = 1; ind[0] = 2; ind[1] = 5; break;
269             case 5:
270             default: type = Tetrahedron::TYPE_O;
271          }
272    }
273 
274    if (j)
275    {
276       mfem::Swap(indices[0], indices[1]);
277       mfem::Swap(indices[2], indices[3]);
278    }
279 
280    CreateRefinementFlag(ind, type);
281 }
282 
283 // static method
GetPointMatrix(unsigned transform,DenseMatrix & pm)284 void Tetrahedron::GetPointMatrix(unsigned transform, DenseMatrix &pm)
285 {
286    double *a = &pm(0,0), *b = &pm(0,1), *c = &pm(0,2), *d = &pm(0,3);
287 
288    // initialize to identity
289    a[0] = 0.0, a[1] = 0.0, a[2] = 0.0;
290    b[0] = 1.0, b[1] = 0.0, b[2] = 0.0;
291    c[0] = 0.0, c[1] = 1.0, c[2] = 0.0;
292    d[0] = 0.0, d[1] = 0.0, d[2] = 1.0;
293 
294    int chain[12], n = 0;
295    while (transform)
296    {
297       chain[n++] = (transform & 7) - 1;
298       transform >>= 3;
299    }
300 
301    /* The transformations and orientations here match the six cases in
302       Mesh::Bisection for tetrahedra. */
303    while (n)
304    {
305 #define ASGN(a, b) (a[0] = b[0], a[1] = b[1], a[2] = b[2])
306 #define SWAP(a, b) for (int i = 0; i < 3; i++) { std::swap(a[i], b[i]); }
307 #define AVG(a, b, c) for (int i = 0; i < 3; i++) { a[i] = (b[i]+c[i])*0.5; }
308 
309       double e[3];
310       AVG(e, a, b);
311       switch (chain[--n])
312       {
313          case 0: ASGN(b, c); ASGN(c, d); break;
314          case 1: ASGN(a, c); ASGN(c, d); break;
315          case 2: ASGN(b, a); ASGN(a, d); break;
316          case 3: ASGN(a, b); ASGN(b, d); break;
317          case 4: SWAP(a, c); ASGN(b, d); break;
318          case 5: SWAP(b, c); ASGN(a, d); break;
319          default:
320             MFEM_ABORT("Invalid transform.");
321       }
322       ASGN(d, e);
323    }
324 }
325 
GetVertices(Array<int> & v) const326 void Tetrahedron::GetVertices(Array<int> &v) const
327 {
328    v.SetSize(4);
329    for (int i = 0; i < 4; i++)
330    {
331       v[i] = indices[i];
332    }
333 }
334 
Duplicate(Mesh * m) const335 Element *Tetrahedron::Duplicate(Mesh *m) const
336 {
337 #ifdef MFEM_USE_MEMALLOC
338    Tetrahedron *tet = m->TetMemory.Alloc();
339 #else
340    Tetrahedron *tet = new Tetrahedron;
341 #endif
342    tet->SetVertices(indices);
343    tet->SetAttribute(attribute);
344    tet->SetRefinementFlag(refinement_flag);
345    return tet;
346 }
347 
348 }
349