cell_tet4.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
27 namespace libMesh
28 {
29 
30 
31 
32 // ------------------------------------------------------------
33 // Tet4 class static member initializations
34 const unsigned int Tet4::side_nodes_map[4][3] =
35 {
36  {0, 2, 1}, // Side 0
37  {0, 1, 3}, // Side 1
38  {1, 2, 3}, // Side 2
39  {2, 0, 3} // Side 3
40 };
41 
42 const unsigned int Tet4::edge_nodes_map[6][2] =
43 {
44  {0, 1}, // Side 0
45  {1, 2}, // Side 1
46  {0, 2}, // Side 2
47  {0, 3}, // Side 3
48  {1, 3}, // Side 4
49  {2, 3} // Side 5
50 };
51 
52 
53 // ------------------------------------------------------------
54 // Tet4 class member functions
55 
56 bool Tet4::is_vertex(const unsigned int) const
57 {
58  return true;
59 }
60 
61 bool Tet4::is_edge(const unsigned int) const
62 {
63  return false;
64 }
65 
66 bool Tet4::is_face(const unsigned int) const
67 {
68  return false;
69 }
70 
71 bool Tet4::is_node_on_edge(const unsigned int n,
72  const unsigned int e) const
73 {
74  libmesh_assert_less (e, n_edges());
75  for (unsigned int i = 0; i != 2; ++i)
76  if (edge_nodes_map[e][i] == n)
77  return true;
78  return false;
79 }
80 
81 
82 
83 
84 #ifdef LIBMESH_ENABLE_AMR
85 
86 // This function only works if LIBMESH_ENABLE_AMR...
87 bool Tet4::is_child_on_side(const unsigned int c,
88  const unsigned int s) const
89 {
90  // OK, for the Tet4, this is pretty obvious... it is sets of nodes
91  // not equal to the current node. But if we want this algorithm to
92  // be generic and work for Tet10 also it helps to do it this way.
93  const unsigned int nodes_opposite[4][3] =
94  {
95  {1,2,3}, // nodes opposite node 0
96  {0,2,3}, // nodes opposite node 1
97  {0,1,3}, // nodes opposite node 2
98  {0,1,2} // nodes opposite node 3
99  };
100 
101  // Call the base class helper function
102  return Tet::is_child_on_side_helper(c, s, nodes_opposite);
103 }
104 
105 #else
106 
107 bool Tet4::is_child_on_side(const unsigned int /*c*/,
108  const unsigned int /*s*/) const
109 {
110  libmesh_not_implemented();
111  return false;
112 }
113 
114 #endif //LIBMESH_ENABLE_AMR
115 
116 
117 
118 
119 bool Tet4::is_node_on_side(const unsigned int n,
120  const unsigned int s) const
121 {
122  libmesh_assert_less (s, n_sides());
123  for (unsigned int i = 0; i != 3; ++i)
124  if (side_nodes_map[s][i] == n)
125  return true;
126  return false;
127 }
128 
129 AutoPtr<Elem> Tet4::build_side (const unsigned int i,
130  bool proxy) const
131 {
132  libmesh_assert_less (i, this->n_sides());
133 
134  if (proxy)
135  {
136  AutoPtr<Elem> ap(new Side<Tri3,Tet4>(this,i));
137  return ap;
138  }
139 
140  else
141  {
142  AutoPtr<Elem> face(new Tri3);
143  face->subdomain_id() = this->subdomain_id();
144 
145  switch (i)
146  {
147  case 0:
148  {
149  face->set_node(0) = this->get_node(0);
150  face->set_node(1) = this->get_node(2);
151  face->set_node(2) = this->get_node(1);
152 
153  return face;
154  }
155  case 1:
156  {
157  face->set_node(0) = this->get_node(0);
158  face->set_node(1) = this->get_node(1);
159  face->set_node(2) = this->get_node(3);
160 
161  return face;
162  }
163  case 2:
164  {
165  face->set_node(0) = this->get_node(1);
166  face->set_node(1) = this->get_node(2);
167  face->set_node(2) = this->get_node(3);
168 
169  return face;
170  }
171  case 3:
172  {
173  face->set_node(0) = this->get_node(2);
174  face->set_node(1) = this->get_node(0);
175  face->set_node(2) = this->get_node(3);
176 
177  return face;
178  }
179  default:
180  {
181  libmesh_error();
182  }
183  }
184  }
185 
186  // We'll never get here.
187  libmesh_error();
188  AutoPtr<Elem> ap(NULL); return ap;
189 }
190 
191 
192 AutoPtr<Elem> Tet4::build_edge (const unsigned int i) const
193 {
194  libmesh_assert_less (i, this->n_edges());
195 
196  return AutoPtr<Elem>(new SideEdge<Edge2,Tet4>(this,i));
197 }
198 
199 
200 void Tet4::connectivity(const unsigned int libmesh_dbg_var(sc),
201  const IOPackage iop,
202  std::vector<dof_id_type>& conn) const
203 {
205  libmesh_assert_less (sc, this->n_sub_elem());
206  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
207 
208 
209  switch (iop)
210  {
211  case TECPLOT:
212  {
213  conn.resize(8);
214  conn[0] = this->node(0)+1;
215  conn[1] = this->node(1)+1;
216  conn[2] = this->node(2)+1;
217  conn[3] = this->node(2)+1;
218  conn[4] = this->node(3)+1;
219  conn[5] = this->node(3)+1;
220  conn[6] = this->node(3)+1;
221  conn[7] = this->node(3)+1;
222  return;
223  }
224 
225  case VTK:
226  {
227  conn.resize(4);
228  conn[0] = this->node(0);
229  conn[1] = this->node(1);
230  conn[2] = this->node(2);
231  conn[3] = this->node(3);
232  return;
233  }
234 
235  default:
236  libmesh_error();
237  }
238 
239  libmesh_error();
240 }
241 
242 
243 
244 #ifdef LIBMESH_ENABLE_AMR
245 
246 const float Tet4::_embedding_matrix[8][4][4] =
247  {
248  // embedding matrix for child 0
249  {
250  // 0 1 2 3
251  {1.0, 0.0, 0.0, 0.0}, // 0
252  {0.5, 0.5, 0.0, 0.0}, // 1
253  {0.5, 0.0, 0.5, 0.0}, // 2
254  {0.5, 0.0, 0.0, 0.5} // 3
255  },
256 
257  // embedding matrix for child 1
258  {
259  // 0 1 2 3
260  {0.5, 0.5, 0.0, 0.0}, // 0
261  {0.0, 1.0, 0.0, 0.0}, // 1
262  {0.0, 0.5, 0.5, 0.0}, // 2
263  {0.0, 0.5, 0.0, 0.5} // 3
264  },
265 
266  // embedding matrix for child 2
267  {
268  // 0 1 2 3
269  {0.5, 0.0, 0.5, 0.0}, // 0
270  {0.0, 0.5, 0.5, 0.0}, // 1
271  {0.0, 0.0, 1.0, 0.0}, // 2
272  {0.0, 0.0, 0.5, 0.5} // 3
273  },
274 
275  // embedding matrix for child 3
276  {
277  // 0 1 2 3
278  {0.5, 0.0, 0.0, 0.5}, // 0
279  {0.0, 0.5, 0.0, 0.5}, // 1
280  {0.0, 0.0, 0.5, 0.5}, // 2
281  {0.0, 0.0, 0.0, 1.0} // 3
282  },
283 
284  // embedding matrix for child 4
285  {
286  // 0 1 2 3
287  {0.5, 0.5, 0.0, 0.0}, // 0
288  {0.0, 0.5, 0.0, 0.5}, // 1
289  {0.5, 0.0, 0.5, 0.0}, // 2
290  {0.5, 0.0, 0.0, 0.5} // 3
291  },
292 
293  // embedding matrix for child 5
294  {
295  // 0 1 2 3
296  {0.5, 0.5, 0.0, 0.0}, // 0
297  {0.0, 0.5, 0.5, 0.0}, // 1
298  {0.5, 0.0, 0.5, 0.0}, // 2
299  {0.0, 0.5, 0.0, 0.5} // 3
300  },
301 
302  // embedding matrix for child 6
303  {
304  // 0 1 2 3
305  {0.5, 0.0, 0.5, 0.0}, // 0
306  {0.0, 0.5, 0.5, 0.0}, // 1
307  {0.0, 0.0, 0.5, 0.5}, // 2
308  {0.0, 0.5, 0.0, 0.5} // 3
309  },
310 
311  // embedding matrix for child 7
312  {
313  // 0 1 2 3
314  {0.5, 0.0, 0.5, 0.0}, // 0
315  {0.0, 0.5, 0.0, 0.5}, // 1
316  {0.0, 0.0, 0.5, 0.5}, // 2
317  {0.5, 0.0, 0.0, 0.5} // 3
318  }
319  };
320 
321 #endif // #ifdef LIBMESH_ENABLE_AMR
322 
323 
324 
325 
326 
328 {
329  // The volume of a tetrahedron is 1/6 the box product formed
330  // by its base and apex vectors
331  Point a ( *this->get_node(3) - *this->get_node(0) );
332 
333  // b is the vector pointing from 0 to 1
334  Point b ( *this->get_node(1) - *this->get_node(0) );
335 
336  // c is the vector pointing from 0 to 2
337  Point c ( *this->get_node(2) - *this->get_node(0) );
338 
339  return (1.0 / 6.0) * (a * (b.cross(c)));
340 }
341 
342 
343 
344 
345 std::pair<Real, Real> Tet4::min_and_max_angle() const
346 {
347  Point n[4];
348 
349  // Compute the outward normal vectors on each face
350  n[0] = (this->point(2) - this->point(0)).cross(this->point(1) - this->point(0));
351  n[1] = (this->point(1) - this->point(0)).cross(this->point(3) - this->point(0));
352  n[2] = (this->point(2) - this->point(1)).cross(this->point(3) - this->point(1));
353  n[3] = (this->point(0) - this->point(2)).cross(this->point(3) - this->point(2));
354 
355  Real dihedral_angles[6]; // 01, 02, 03, 12, 13, 23
356 
357  // Compute dihedral angles
358  for (unsigned int k=0,i=0; i<4; ++i)
359  for (unsigned int j=i+1; j<4; ++j,k+=1)
360  dihedral_angles[k] = std::acos(n[i]*n[j] / n[i].size() / n[j].size()); // return value is between 0 and PI
361 
362  // Return max/min dihedral angles
363  return std::make_pair(*std::min_element(dihedral_angles, dihedral_angles+6),
364  *std::max_element(dihedral_angles, dihedral_angles+6));
365 
366 }
367 
368 
369 
370 #ifdef LIBMESH_ENABLE_AMR
371 float Tet4::embedding_matrix (const unsigned int i,
372  const unsigned int j,
373  const unsigned int k) const
374 {
375  // Choose an optimal diagonal, if one has not already been selected
376  this->choose_diagonal();
377 
378  // Permuted j and k indices
379  unsigned int
380  jp=j,
381  kp=k;
382 
383  if ((i>3) && (this->_diagonal_selection!=DIAG_02_13))
384  {
385  // Just the enum value cast to an unsigned int...
386  const unsigned ds = static_cast<unsigned int>(this->_diagonal_selection);
387 
388  // Permute j, k:
389  // ds==1 ds==2
390  // 0 -> 1 0 -> 2
391  // 1 -> 2 1 -> 0
392  // 2 -> 0 2 -> 1
393  if (jp != 3)
394  jp = (jp+ds)%3;
395 
396  if (kp != 3)
397  kp = (kp+ds)%3;
398  }
399 
400  // Debugging
401  // libMesh::err << "Selected diagonal " << _diagonal_selection << std::endl;
402  // libMesh::err << "j=" << j << std::endl;
403  // libMesh::err << "k=" << k << std::endl;
404  // libMesh::err << "jp=" << jp << std::endl;
405  // libMesh::err << "kp=" << kp << std::endl;
406 
407  // Call embedding matrx with permuted indices
408  return this->_embedding_matrix[i][jp][kp];
409 }
410 
411 
412 
413 
414 
415 
416 // void Tet4::reselect_diagonal (const Diagonal diag)
417 // {
418 // /* Make sure that the element has just been refined. */
419 // libmesh_assert(_children);
420 // libmesh_assert_equal_to (n_children(), 8);
421 // libmesh_assert_equal_to (_children[0]->refinement_flag(), JUST_REFINED);
422 // libmesh_assert_equal_to (_children[1]->refinement_flag(), JUST_REFINED);
423 // libmesh_assert_equal_to (_children[2]->refinement_flag(), JUST_REFINED);
424 // libmesh_assert_equal_to (_children[3]->refinement_flag(), JUST_REFINED);
425 // libmesh_assert_equal_to (_children[4]->refinement_flag(), JUST_REFINED);
426 // libmesh_assert_equal_to (_children[5]->refinement_flag(), JUST_REFINED);
427 // libmesh_assert_equal_to (_children[6]->refinement_flag(), JUST_REFINED);
428 // libmesh_assert_equal_to (_children[7]->refinement_flag(), JUST_REFINED);
429 //
430 // /* Check whether anything has to be changed. */
431 // if (_diagonal_selection!=diag)
432 // {
433 // /* Set new diagonal selection. */
434 // _diagonal_selection = diag;
435 //
436 // /* The first four children do not have to be changed. For the
437 // others, only the nodes have to be changed. Note that we have
438 // to keep track of the nodes ourselves since there is no \p
439 // MeshRefinement object with a valid \p _new_nodes_map
440 // available. */
441 // for (unsigned int c=4; c<this->n_children(); c++)
442 // {
443 // Elem *child = this->child(c);
444 // for (unsigned int nc=0; nc<child->n_nodes(); nc++)
445 // {
446 // /* Unassign the current node. */
447 // child->set_node(nc) = NULL;
448 //
449 // /* We have to find the correct new node now. We know
450 // that it exists somewhere. We make use of the fact
451 // that the embedding matrix for these children consists
452 // of entries 0.0 and 0.5 only. Also, we make use of
453 // the properties of the embedding matrix for the first
454 // (unchanged) four children, which allow us to use a
455 // simple mechanism to find the required node. */
456 //
457 //
458 // unsigned int first_05_in_embedding_matrix = libMesh::invalid_uint;
459 // for (unsigned int n=0; n<this->n_nodes(); n++)
460 // {
461 // if (this->embedding_matrix(c,nc,n) != 0.0)
462 // {
463 // /* It must be 0.5 then. Check whether it's the
464 // first or second time that we get a 0.5
465 // value. */
466 // if (first_05_in_embedding_matrix==libMesh::invalid_uint)
467 // {
468 // /* First time, so just remeber this position. */
469 // first_05_in_embedding_matrix = n;
470 // }
471 // else
472 // {
473 // /* Second time, so we know now which node to
474 // use. */
475 // child->set_node(nc) = this->child(n)->get_node(first_05_in_embedding_matrix);
476 // }
477 //
478 // }
479 // }
480 //
481 // /* Make sure that a node has been found. */
482 // libmesh_assert(child->get_node(nc));
483 // }
484 // }
485 // }
486 // }
487 
488 
489 
490 // void Tet4::reselect_optimal_diagonal (const Diagonal exclude_this)
491 // {
492 // Real diag_01_23 = (this->point(0)+this->point(1)-this->point(2)-this->point(3)).size_sq();
493 // Real diag_02_13 = (this->point(0)-this->point(1)+this->point(2)-this->point(3)).size_sq();
494 // Real diag_03_12 = (this->point(0)-this->point(1)-this->point(2)+this->point(3)).size_sq();
495 //
496 // Diagonal use_this = INVALID_DIAG;
497 // switch (exclude_this)
498 // {
499 // case DIAG_01_23:
500 // use_this = DIAG_02_13;
501 // if (diag_03_12 < diag_02_13)
502 // {
503 // use_this = DIAG_03_12;
504 // }
505 // break;
506 //
507 // case DIAG_02_13:
508 // use_this = DIAG_03_12;
509 // if (diag_01_23 < diag_03_12)
510 // {
511 // use_this = DIAG_01_23;
512 // }
513 // break;
514 //
515 // case DIAG_03_12:
516 // use_this = DIAG_02_13;
517 // if (diag_01_23 < diag_02_13)
518 // {
519 // use_this = DIAG_01_23;
520 // }
521 // break;
522 //
523 // default:
524 // use_this = DIAG_02_13;
525 // if (diag_01_23 < diag_02_13 || diag_03_12 < diag_02_13)
526 // {
527 // if (diag_01_23 < diag_03_12)
528 // {
529 // use_this = DIAG_01_23;
530 // }
531 // else
532 // {
533 // use_this = DIAG_03_12;
534 // }
535 // }
536 // break;
537 // }
538 //
539 // reselect_diagonal (use_this);
540 // }
541 #endif // #ifdef LIBMESH_ENABLE_AMR
542 
543 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo