MFEM  v4.1.0 Finite element discretization library
tetrahedron.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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
15
16 namespace mfem
17 {
18
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
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
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
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
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  {
83  if (e1 == 2 && e2 == 1) { break; }
84  // if (e1 == 3 && e2 == 4) break;
85  mfem_error("Error in Tetrahedron::CreateRefinementFlag(...) #1");
86  break;
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;
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
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
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
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
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
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
327 {
328  v.SetSize(4);
329  for (int i = 0; i < 4; i++)
330  {
331  v[i] = indices[i];
332  }
333 }
334
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);
345  return tet;
346 }
347
348 }
void Init(int ind1, int ind2, int ind3, int ind4, int attr=1, int ref_flag=0)
Initialize the vertex indices and the attribute of a Tetrahedron.
Definition: tetrahedron.cpp:43
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual void SetVertices(const int *ind)
Set the vertices according to the given input.
virtual Element * Duplicate(Mesh *m) const
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:153
double b
Definition: lissajous.cpp:42
void CreateRefinementFlag(int refinement_edges[2], int type, int flag=0)
Definition: tetrahedron.cpp:71
Data type tetrahedron element.
Definition: tetrahedron.hpp:22
void Swap(Array< T > &, Array< T > &)
Definition: array.hpp:592
int attribute
Element&#39;s attribute (specifying material property, etc).
Definition: element.hpp:33
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
virtual int * GetVertices()
Definition: tetrahedron.hpp:94
double a
Definition: lissajous.cpp:41
virtual int NeedRefinement(HashTable< Hashed2 > &v_to_v) const
Return 1 if the element needs refinement in order to get conforming mesh.
static void GetPointMatrix(unsigned transform, DenseMatrix &pm)
Calculate point matrix corresponding to a chain of transformations.
void SetAttribute(const int attr)
Set element&#39;s attribute.
Definition: element.hpp:58
void GetMarkedFace(const int face, int *fv)
MemAlloc< Tetrahedron, 1024 > TetMemory
Definition: mesh.hpp:173
void SetRefinementFlag(int rf)
Definition: tetrahedron.hpp:68
void ParseRefinementFlag(int refinement_edges[2], int &type, int &flag)
Definition: tetrahedron.cpp:55
int FindId(int p1, int p2) const
Find id of item whose parents are p1, p2... Return -1 if it doesn&#39;t exist.
Definition: hash.hpp:370
virtual void MarkEdge(const DSTable &v_to_v, const int *length)
Abstract data type element.
Definition: element.hpp:28