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