libMesh
cell_tet4.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 
21 // Local includes
22 #include "libmesh/side.h"
23 #include "libmesh/cell_tet4.h"
24 #include "libmesh/edge_edge2.h"
25 #include "libmesh/face_tri3.h"
26 #include "libmesh/tensor_value.h"
27 
28 namespace libMesh
29 {
30 
31 
32 
33 // ------------------------------------------------------------
34 // Tet4 class static member initializations
35 const unsigned int Tet4::side_nodes_map[4][3] =
36  {
37  {0, 2, 1}, // Side 0
38  {0, 1, 3}, // Side 1
39  {1, 2, 3}, // Side 2
40  {2, 0, 3} // Side 3
41  };
42 
43 const unsigned int Tet4::edge_nodes_map[6][2] =
44  {
45  {0, 1}, // Edge 0
46  {1, 2}, // Edge 1
47  {0, 2}, // Edge 2
48  {0, 3}, // Edge 3
49  {1, 3}, // Edge 4
50  {2, 3} // Edge 5
51  };
52 
53 
54 // ------------------------------------------------------------
55 // Tet4 class member functions
56 
57 bool Tet4::is_vertex(const unsigned int) const
58 {
59  return true;
60 }
61 
62 bool Tet4::is_edge(const unsigned int) const
63 {
64  return false;
65 }
66 
67 bool Tet4::is_face(const unsigned int) const
68 {
69  return false;
70 }
71 
72 bool Tet4::is_node_on_edge(const unsigned int n,
73  const unsigned int e) const
74 {
75  libmesh_assert_less (e, n_edges());
76  for (unsigned int i = 0; i != 2; ++i)
77  if (edge_nodes_map[e][i] == n)
78  return true;
79  return false;
80 }
81 
82 
83 
84 
85 #ifdef LIBMESH_ENABLE_AMR
86 
87 // This function only works if LIBMESH_ENABLE_AMR...
88 bool Tet4::is_child_on_side(const unsigned int c,
89  const unsigned int s) const
90 {
91  // OK, for the Tet4, this is pretty obvious... it is sets of nodes
92  // not equal to the current node. But if we want this algorithm to
93  // be generic and work for Tet10 also it helps to do it this way.
94  const unsigned int nodes_opposite[4][3] =
95  {
96  {1,2,3}, // nodes opposite node 0
97  {0,2,3}, // nodes opposite node 1
98  {0,1,3}, // nodes opposite node 2
99  {0,1,2} // nodes opposite node 3
100  };
101 
102  // Call the base class helper function
103  return Tet::is_child_on_side_helper(c, s, nodes_opposite);
104 }
105 
106 #else
107 
108 bool Tet4::is_child_on_side(const unsigned int /*c*/,
109  const unsigned int /*s*/) const
110 {
111  libmesh_not_implemented();
112  return false;
113 }
114 
115 #endif //LIBMESH_ENABLE_AMR
116 
117 
118 
119 
120 bool Tet4::is_node_on_side(const unsigned int n,
121  const unsigned int s) const
122 {
123  libmesh_assert_less (s, n_sides());
124  for (unsigned int i = 0; i != 3; ++i)
125  if (side_nodes_map[s][i] == n)
126  return true;
127  return false;
128 }
129 
130 UniquePtr<Elem> Tet4::build_side_ptr (const unsigned int i,
131  bool proxy)
132 {
133  libmesh_assert_less (i, this->n_sides());
134 
135  if (proxy)
136  return UniquePtr<Elem>(new Side<Tri3,Tet4>(this,i));
137 
138  else
139  {
140  Elem * face = new Tri3;
141  face->subdomain_id() = this->subdomain_id();
142 
143  for (unsigned n=0; n<face->n_nodes(); ++n)
144  face->set_node(n) = this->node_ptr(Tet4::side_nodes_map[i][n]);
145 
146  return UniquePtr<Elem>(face);
147  }
148 
149  libmesh_error_msg("We'll never get here!");
150  return UniquePtr<Elem>();
151 }
152 
153 
154 UniquePtr<Elem> Tet4::build_edge_ptr (const unsigned int i)
155 {
156  libmesh_assert_less (i, this->n_edges());
157 
158  return UniquePtr<Elem>(new SideEdge<Edge2,Tet4>(this,i));
159 }
160 
161 
162 void Tet4::connectivity(const unsigned int libmesh_dbg_var(sc),
163  const IOPackage iop,
164  std::vector<dof_id_type> & conn) const
165 {
166  libmesh_assert(_nodes);
167  libmesh_assert_less (sc, this->n_sub_elem());
168  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
169 
170 
171  switch (iop)
172  {
173  case TECPLOT:
174  {
175  conn.resize(8);
176  conn[0] = this->node_id(0)+1;
177  conn[1] = this->node_id(1)+1;
178  conn[2] = this->node_id(2)+1;
179  conn[3] = this->node_id(2)+1;
180  conn[4] = this->node_id(3)+1;
181  conn[5] = this->node_id(3)+1;
182  conn[6] = this->node_id(3)+1;
183  conn[7] = this->node_id(3)+1;
184  return;
185  }
186 
187  case VTK:
188  {
189  conn.resize(4);
190  conn[0] = this->node_id(0);
191  conn[1] = this->node_id(1);
192  conn[2] = this->node_id(2);
193  conn[3] = this->node_id(3);
194  return;
195  }
196 
197  default:
198  libmesh_error_msg("Unsupported IO package " << iop);
199  }
200 }
201 
202 
203 
204 #ifdef LIBMESH_ENABLE_AMR
205 
206 const float Tet4::_embedding_matrix[8][4][4] =
207  {
208  // embedding matrix for child 0
209  {
210  // 0 1 2 3
211  {1.0, 0.0, 0.0, 0.0}, // 0
212  {0.5, 0.5, 0.0, 0.0}, // 1
213  {0.5, 0.0, 0.5, 0.0}, // 2
214  {0.5, 0.0, 0.0, 0.5} // 3
215  },
216 
217  // embedding matrix for child 1
218  {
219  // 0 1 2 3
220  {0.5, 0.5, 0.0, 0.0}, // 0
221  {0.0, 1.0, 0.0, 0.0}, // 1
222  {0.0, 0.5, 0.5, 0.0}, // 2
223  {0.0, 0.5, 0.0, 0.5} // 3
224  },
225 
226  // embedding matrix for child 2
227  {
228  // 0 1 2 3
229  {0.5, 0.0, 0.5, 0.0}, // 0
230  {0.0, 0.5, 0.5, 0.0}, // 1
231  {0.0, 0.0, 1.0, 0.0}, // 2
232  {0.0, 0.0, 0.5, 0.5} // 3
233  },
234 
235  // embedding matrix for child 3
236  {
237  // 0 1 2 3
238  {0.5, 0.0, 0.0, 0.5}, // 0
239  {0.0, 0.5, 0.0, 0.5}, // 1
240  {0.0, 0.0, 0.5, 0.5}, // 2
241  {0.0, 0.0, 0.0, 1.0} // 3
242  },
243 
244  // embedding matrix for child 4
245  {
246  // 0 1 2 3
247  {0.5, 0.5, 0.0, 0.0}, // 0
248  {0.0, 0.5, 0.0, 0.5}, // 1
249  {0.5, 0.0, 0.5, 0.0}, // 2
250  {0.5, 0.0, 0.0, 0.5} // 3
251  },
252 
253  // embedding matrix for child 5
254  {
255  // 0 1 2 3
256  {0.5, 0.5, 0.0, 0.0}, // 0
257  {0.0, 0.5, 0.5, 0.0}, // 1
258  {0.5, 0.0, 0.5, 0.0}, // 2
259  {0.0, 0.5, 0.0, 0.5} // 3
260  },
261 
262  // embedding matrix for child 6
263  {
264  // 0 1 2 3
265  {0.5, 0.0, 0.5, 0.0}, // 0
266  {0.0, 0.5, 0.5, 0.0}, // 1
267  {0.0, 0.0, 0.5, 0.5}, // 2
268  {0.0, 0.5, 0.0, 0.5} // 3
269  },
270 
271  // embedding matrix for child 7
272  {
273  // 0 1 2 3
274  {0.5, 0.0, 0.5, 0.0}, // 0
275  {0.0, 0.5, 0.0, 0.5}, // 1
276  {0.0, 0.0, 0.5, 0.5}, // 2
277  {0.5, 0.0, 0.0, 0.5} // 3
278  }
279  };
280 
281 #endif // #ifdef LIBMESH_ENABLE_AMR
282 
283 
284 
285 
286 
287 Real Tet4::volume () const
288 {
289  // The volume of a tetrahedron is 1/6 the box product formed
290  // by its base and apex vectors
291  Point a = point(3) - point(0);
292 
293  // b is the vector pointing from 0 to 1
294  Point b = point(1) - point(0);
295 
296  // c is the vector pointing from 0 to 2
297  Point c = point(2) - point(0);
298 
299  return triple_product(a, b, c) / 6.;
300 }
301 
302 
303 
304 
305 std::pair<Real, Real> Tet4::min_and_max_angle() const
306 {
307  Point n[4];
308 
309  // Compute the outward normal vectors on each face
310  n[0] = (this->point(2) - this->point(0)).cross(this->point(1) - this->point(0));
311  n[1] = (this->point(1) - this->point(0)).cross(this->point(3) - this->point(0));
312  n[2] = (this->point(2) - this->point(1)).cross(this->point(3) - this->point(1));
313  n[3] = (this->point(0) - this->point(2)).cross(this->point(3) - this->point(2));
314 
315  Real dihedral_angles[6]; // 01, 02, 03, 12, 13, 23
316 
317  // Compute dihedral angles
318  for (unsigned int k=0,i=0; i<4; ++i)
319  for (unsigned int j=i+1; j<4; ++j,k+=1)
320  dihedral_angles[k] = std::acos(n[i]*n[j] / n[i].norm() / n[j].norm()); // return value is between 0 and PI
321 
322  // Return max/min dihedral angles
323  return std::make_pair(*std::min_element(dihedral_angles, dihedral_angles+6),
324  *std::max_element(dihedral_angles, dihedral_angles+6));
325 
326 }
327 
328 
329 
330 dof_id_type Tet4::key () const
331 {
332  return this->compute_key(this->node_id(0),
333  this->node_id(1),
334  this->node_id(2),
335  this->node_id(3));
336 }
337 
338 
339 
340 bool Tet4::contains_point (const Point & p, Real tol) const
341 {
342  // See the response by Tony Noe on this thread.
343  // http://bbs.dartmouth.edu/~fangq/MATH/download/source/Point_in_Tetrahedron.htm
344  Point
345  col1 = point(1) - point(0),
346  col2 = point(2) - point(0),
347  col3 = point(3) - point(0);
348 
349  Point r;
350 
351  libmesh_try
352  {
353  RealTensorValue(col1(0), col2(0), col3(0),
354  col1(1), col2(1), col3(1),
355  col1(2), col2(2), col3(2)).solve(p - point(0), r);
356  }
357  libmesh_catch (ConvergenceFailure &)
358  {
359  // The Jacobian was singular, therefore the Tet had
360  // zero volume. In this degenerate case, it is not
361  // possible for the Tet to contain any points.
362  return false;
363  }
364 
365  // The point p is inside the tetrahedron if r1 + r2 + r3 < 1 and
366  // 0 <= ri <= 1 for i = 1,2,3.
367  return
368  r(0) > -tol &&
369  r(1) > -tol &&
370  r(2) > -tol &&
371  r(0) + r(1) + r(2) < 1.0 + tol;
372 }
373 
374 
375 
376 #ifdef LIBMESH_ENABLE_AMR
377 float Tet4::embedding_matrix (const unsigned int i,
378  const unsigned int j,
379  const unsigned int k) const
380 {
381  // Choose an optimal diagonal, if one has not already been selected
382  this->choose_diagonal();
383 
384  // Permuted j and k indices
385  unsigned int
386  jp=j,
387  kp=k;
388 
389  if ((i>3) && (this->_diagonal_selection!=DIAG_02_13))
390  {
391  // Just the enum value cast to an unsigned int...
392  const unsigned ds = static_cast<unsigned int>(this->_diagonal_selection);
393 
394  // Permute j, k:
395  // ds==1 ds==2
396  // 0 -> 1 0 -> 2
397  // 1 -> 2 1 -> 0
398  // 2 -> 0 2 -> 1
399  if (jp != 3)
400  jp = (jp+ds)%3;
401 
402  if (kp != 3)
403  kp = (kp+ds)%3;
404  }
405 
406  // Debugging
407  // libMesh::err << "Selected diagonal " << _diagonal_selection << std::endl;
408  // libMesh::err << "j=" << j << std::endl;
409  // libMesh::err << "k=" << k << std::endl;
410  // libMesh::err << "jp=" << jp << std::endl;
411  // libMesh::err << "kp=" << kp << std::endl;
412 
413  // Call embedding matrix with permuted indices
414  return this->_embedding_matrix[i][jp][kp];
415 }
416 
417 
418 
419 // void Tet4::reselect_diagonal (const Diagonal diag)
420 // {
421 // /* Make sure that the element has just been refined. */
422 // libmesh_assert(_children);
423 // libmesh_assert_equal_to (n_children(), 8);
424 // libmesh_assert_equal_to (_children[0]->refinement_flag(), JUST_REFINED);
425 // libmesh_assert_equal_to (_children[1]->refinement_flag(), JUST_REFINED);
426 // libmesh_assert_equal_to (_children[2]->refinement_flag(), JUST_REFINED);
427 // libmesh_assert_equal_to (_children[3]->refinement_flag(), JUST_REFINED);
428 // libmesh_assert_equal_to (_children[4]->refinement_flag(), JUST_REFINED);
429 // libmesh_assert_equal_to (_children[5]->refinement_flag(), JUST_REFINED);
430 // libmesh_assert_equal_to (_children[6]->refinement_flag(), JUST_REFINED);
431 // libmesh_assert_equal_to (_children[7]->refinement_flag(), JUST_REFINED);
432 //
433 // /* Check whether anything has to be changed. */
434 // if (_diagonal_selection!=diag)
435 // {
436 // /* Set new diagonal selection. */
437 // _diagonal_selection = diag;
438 //
439 // /* The first four children do not have to be changed. For the
440 // others, only the nodes have to be changed. Note that we have
441 // to keep track of the nodes ourselves since there is no \p
442 // MeshRefinement object with a valid \p _new_nodes_map
443 // available. */
444 // for (unsigned int c=4; c<this->n_children(); c++)
445 // {
446 // Elem * child = this->child_ptr(c);
447 // for (unsigned int nc=0; nc<child->n_nodes(); nc++)
448 // {
449 // /* Unassign the current node. */
450 // child->set_node(nc) = libmesh_nullptr;
451 //
452 // /* We have to find the correct new node now. We know
453 // that it exists somewhere. We make use of the fact
454 // that the embedding matrix for these children consists
455 // of entries 0.0 and 0.5 only. Also, we make use of
456 // the properties of the embedding matrix for the first
457 // (unchanged) four children, which allow us to use a
458 // simple mechanism to find the required node. */
459 //
460 //
461 // unsigned int first_05_in_embedding_matrix = libMesh::invalid_uint;
462 // for (unsigned int n=0; n<this->n_nodes(); n++)
463 // {
464 // if (this->embedding_matrix(c,nc,n) != 0.0)
465 // {
466 // /* It must be 0.5 then. Check whether it's the
467 // first or second time that we get a 0.5
468 // value. */
469 // if (first_05_in_embedding_matrix==libMesh::invalid_uint)
470 // {
471 // /* First time, so just remember this position. */
472 // first_05_in_embedding_matrix = n;
473 // }
474 // else
475 // {
476 // /* Second time, so we know now which node to
477 // use. */
478 // child->set_node(nc) = this->child_ptr(n)->node_ptr(first_05_in_embedding_matrix);
479 // }
480 //
481 // }
482 // }
483 //
484 // /* Make sure that a node has been found. */
485 // libmesh_assert(child->node_ptr(nc));
486 // }
487 // }
488 // }
489 // }
490 
491 
492 
493 // void Tet4::reselect_optimal_diagonal (const Diagonal exclude_this)
494 // {
495 // Real diag_01_23 = (this->point(0)+this->point(1)-this->point(2)-this->point(3)).norm_sq();
496 // Real diag_02_13 = (this->point(0)-this->point(1)+this->point(2)-this->point(3)).norm_sq();
497 // Real diag_03_12 = (this->point(0)-this->point(1)-this->point(2)+this->point(3)).norm_sq();
498 //
499 // Diagonal use_this = INVALID_DIAG;
500 // switch (exclude_this)
501 // {
502 // case DIAG_01_23:
503 // use_this = DIAG_02_13;
504 // if (diag_03_12 < diag_02_13)
505 // {
506 // use_this = DIAG_03_12;
507 // }
508 // break;
509 //
510 // case DIAG_02_13:
511 // use_this = DIAG_03_12;
512 // if (diag_01_23 < diag_03_12)
513 // {
514 // use_this = DIAG_01_23;
515 // }
516 // break;
517 //
518 // case DIAG_03_12:
519 // use_this = DIAG_02_13;
520 // if (diag_01_23 < diag_02_13)
521 // {
522 // use_this = DIAG_01_23;
523 // }
524 // break;
525 //
526 // default:
527 // use_this = DIAG_02_13;
528 // if (diag_01_23 < diag_02_13 || diag_03_12 < diag_02_13)
529 // {
530 // if (diag_01_23 < diag_03_12)
531 // {
532 // use_this = DIAG_01_23;
533 // }
534 // else
535 // {
536 // use_this = DIAG_03_12;
537 // }
538 // }
539 // break;
540 // }
541 //
542 // reselect_diagonal (use_this);
543 // }
544 #endif // #ifdef LIBMESH_ENABLE_AMR
545 
546 } // namespace libMesh
bool is_child_on_side_helper(const unsigned int c, const unsigned int s, const unsigned int checked_nodes[][3]) const
Called by descendant classes with appropriate data to determine if child c is on side s...
Definition: cell_tet.C:103
The libMesh namespace provides an interface to certain functionality in the library.
libmesh_assert(j)
TensorValue< Real > RealTensorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1051
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
uint8_t dof_id_type
Definition: id_types.h:64