libMesh
cell_tet10.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_tet10.h"
24 #include "libmesh/edge_edge3.h"
25 #include "libmesh/face_tri6.h"
26 
27 namespace libMesh
28 {
29 
30 
31 
32 // ------------------------------------------------------------
33 // Tet10 class static member initializations
34 const unsigned int Tet10::side_nodes_map[4][6] =
35  {
36  {0, 2, 1, 6, 5, 4}, // Side 0
37  {0, 1, 3, 4, 8, 7}, // Side 1
38  {1, 2, 3, 5, 9, 8}, // Side 2
39  {2, 0, 3, 6, 7, 9} // Side 3
40  };
41 
42 const unsigned int Tet10::edge_nodes_map[6][3] =
43  {
44  {0, 1, 4}, // Edge 0
45  {1, 2, 5}, // Edge 1
46  {0, 2, 6}, // Edge 2
47  {0, 3, 7}, // Edge 3
48  {1, 3, 8}, // Edge 4
49  {2, 3, 9} // Edge 5
50  };
51 
52 
53 
54 // ------------------------------------------------------------
55 // Tet10 class member functions
56 
57 bool Tet10::is_vertex(const unsigned int i) const
58 {
59  if (i < 4)
60  return true;
61  return false;
62 }
63 
64 bool Tet10::is_edge(const unsigned int i) const
65 {
66  if (i < 4)
67  return false;
68  return true;
69 }
70 
71 bool Tet10::is_face(const unsigned int) const
72 {
73  return false;
74 }
75 
76 bool Tet10::is_node_on_side(const unsigned int n,
77  const unsigned int s) const
78 {
79  libmesh_assert_less (s, n_sides());
80  for (unsigned int i = 0; i != 6; ++i)
81  if (side_nodes_map[s][i] == n)
82  return true;
83  return false;
84 }
85 
86 bool Tet10::is_node_on_edge(const unsigned int n,
87  const unsigned int e) const
88 {
89  libmesh_assert_less (e, n_edges());
90  for (unsigned int i = 0; i != 3; ++i)
91  if (edge_nodes_map[e][i] == n)
92  return true;
93  return false;
94 }
95 
96 
97 #ifdef LIBMESH_ENABLE_AMR
98 
99 // This function only works if LIBMESH_ENABLE_AMR...
100 bool Tet10::is_child_on_side(const unsigned int c,
101  const unsigned int s) const
102 {
103  // Table of local IDs for the midege nodes on the side opposite a given node.
104  // See the ASCII art in the header file for this class to confirm this.
105  const unsigned int midedge_nodes_opposite[4][3] =
106  {
107  {5,8,9}, // midedge nodes opposite node 0
108  {6,7,9}, // midedge nodes opposite node 1
109  {4,7,8}, // midedge nodes opposite node 2
110  {4,5,6} // midedge nodes opposite node 3
111  };
112 
113  // Call the base class helper function
114  return Tet::is_child_on_side_helper(c, s, midedge_nodes_opposite);
115 }
116 
117 #else
118 
119 bool Tet10::is_child_on_side(const unsigned int /*c*/,
120  const unsigned int /*s*/) const
121 {
122  libmesh_not_implemented();
123  return false;
124 }
125 
126 #endif //LIBMESH_ENABLE_AMR
127 
128 
129 
130 bool Tet10::has_affine_map() const
131 {
132  // Make sure edges are straight
133  if (!this->point(4).relative_fuzzy_equals
134  ((this->point(0) + this->point(1))/2))
135  return false;
136  if (!this->point(5).relative_fuzzy_equals
137  ((this->point(1) + this->point(2))/2))
138  return false;
139  if (!this->point(6).relative_fuzzy_equals
140  ((this->point(2) + this->point(0))/2))
141  return false;
142  if (!this->point(7).relative_fuzzy_equals
143  ((this->point(3) + this->point(0))/2))
144  return false;
145  if (!this->point(8).relative_fuzzy_equals
146  ((this->point(3) + this->point(1))/2))
147  return false;
148  if (!this->point(9).relative_fuzzy_equals
149  ((this->point(3) + this->point(2))/2))
150  return false;
151  return true;
152 }
153 
154 
155 
156 unsigned int Tet10::which_node_am_i(unsigned int side,
157  unsigned int side_node) const
158 {
159  libmesh_assert_less (side, this->n_sides());
160  libmesh_assert_less (side_node, 6);
161 
162  return Tet10::side_nodes_map[side][side_node];
163 }
164 
165 
166 
167 UniquePtr<Elem> Tet10::build_side_ptr (const unsigned int i,
168  bool proxy)
169 {
170  libmesh_assert_less (i, this->n_sides());
171 
172  if (proxy)
173  return UniquePtr<Elem>(new Side<Tri6,Tet10>(this,i));
174 
175  else
176  {
177  Elem * face = new Tri6;
178  face->subdomain_id() = this->subdomain_id();
179 
180  for (unsigned n=0; n<face->n_nodes(); ++n)
181  face->set_node(n) = this->node_ptr(Tet10::side_nodes_map[i][n]);
182 
183  return UniquePtr<Elem>(face);
184  }
185 
186  libmesh_error_msg("We'll never get here!");
187  return UniquePtr<Elem>();
188 }
189 
190 
191 
192 UniquePtr<Elem> Tet10::build_edge_ptr (const unsigned int i)
193 {
194  libmesh_assert_less (i, this->n_edges());
195 
196  return UniquePtr<Elem>(new SideEdge<Edge3,Tet10>(this,i));
197 }
198 
199 
200 
201 void Tet10::connectivity(const unsigned int sc,
202  const IOPackage iop,
203  std::vector<dof_id_type> & conn) const
204 {
205  libmesh_assert(_nodes);
206  libmesh_assert_less (sc, this->n_sub_elem());
207  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
208 
209  switch (iop)
210  {
211  case TECPLOT:
212  {
213  conn.resize(8);
214  switch (sc)
215  {
216 
217 
218  // Linear sub-tet 0
219  case 0:
220 
221  conn[0] = this->node_id(0)+1;
222  conn[1] = this->node_id(4)+1;
223  conn[2] = this->node_id(6)+1;
224  conn[3] = this->node_id(6)+1;
225  conn[4] = this->node_id(7)+1;
226  conn[5] = this->node_id(7)+1;
227  conn[6] = this->node_id(7)+1;
228  conn[7] = this->node_id(7)+1;
229 
230  return;
231 
232  // Linear sub-tet 1
233  case 1:
234 
235  conn[0] = this->node_id(4)+1;
236  conn[1] = this->node_id(1)+1;
237  conn[2] = this->node_id(5)+1;
238  conn[3] = this->node_id(5)+1;
239  conn[4] = this->node_id(8)+1;
240  conn[5] = this->node_id(8)+1;
241  conn[6] = this->node_id(8)+1;
242  conn[7] = this->node_id(8)+1;
243 
244  return;
245 
246  // Linear sub-tet 2
247  case 2:
248 
249  conn[0] = this->node_id(5)+1;
250  conn[1] = this->node_id(2)+1;
251  conn[2] = this->node_id(6)+1;
252  conn[3] = this->node_id(6)+1;
253  conn[4] = this->node_id(9)+1;
254  conn[5] = this->node_id(9)+1;
255  conn[6] = this->node_id(9)+1;
256  conn[7] = this->node_id(9)+1;
257 
258  return;
259 
260  // Linear sub-tet 3
261  case 3:
262 
263  conn[0] = this->node_id(7)+1;
264  conn[1] = this->node_id(8)+1;
265  conn[2] = this->node_id(9)+1;
266  conn[3] = this->node_id(9)+1;
267  conn[4] = this->node_id(3)+1;
268  conn[5] = this->node_id(3)+1;
269  conn[6] = this->node_id(3)+1;
270  conn[7] = this->node_id(3)+1;
271 
272  return;
273 
274  // Linear sub-tet 4
275  case 4:
276 
277  conn[0] = this->node_id(4)+1;
278  conn[1] = this->node_id(8)+1;
279  conn[2] = this->node_id(6)+1;
280  conn[3] = this->node_id(6)+1;
281  conn[4] = this->node_id(7)+1;
282  conn[5] = this->node_id(7)+1;
283  conn[6] = this->node_id(7)+1;
284  conn[7] = this->node_id(7)+1;
285 
286  return;
287 
288  // Linear sub-tet 5
289  case 5:
290 
291  conn[0] = this->node_id(4)+1;
292  conn[1] = this->node_id(5)+1;
293  conn[2] = this->node_id(6)+1;
294  conn[3] = this->node_id(6)+1;
295  conn[4] = this->node_id(8)+1;
296  conn[5] = this->node_id(8)+1;
297  conn[6] = this->node_id(8)+1;
298  conn[7] = this->node_id(8)+1;
299 
300  return;
301 
302  // Linear sub-tet 6
303  case 6:
304 
305  conn[0] = this->node_id(5)+1;
306  conn[1] = this->node_id(9)+1;
307  conn[2] = this->node_id(6)+1;
308  conn[3] = this->node_id(6)+1;
309  conn[4] = this->node_id(8)+1;
310  conn[5] = this->node_id(8)+1;
311  conn[6] = this->node_id(8)+1;
312  conn[7] = this->node_id(8)+1;
313 
314  return;
315 
316  // Linear sub-tet 7
317  case 7:
318 
319  conn[0] = this->node_id(7)+1;
320  conn[1] = this->node_id(6)+1;
321  conn[2] = this->node_id(9)+1;
322  conn[3] = this->node_id(9)+1;
323  conn[4] = this->node_id(8)+1;
324  conn[5] = this->node_id(8)+1;
325  conn[6] = this->node_id(8)+1;
326  conn[7] = this->node_id(8)+1;
327 
328  return;
329 
330 
331  default:
332  libmesh_error_msg("Invalid sc = " << sc);
333  }
334  }
335 
336  case VTK:
337  {
338  conn.resize(10);
339  conn[0] = this->node_id(0);
340  conn[1] = this->node_id(1);
341  conn[2] = this->node_id(2);
342  conn[3] = this->node_id(3);
343  conn[4] = this->node_id(4);
344  conn[5] = this->node_id(5);
345  conn[6] = this->node_id(6);
346  conn[7] = this->node_id(7);
347  conn[8] = this->node_id(8);
348  conn[9] = this->node_id(9);
349  return;
350  /*
351  conn.resize(4);
352  switch (sc)
353  {
354  // Linear sub-tet 0
355  case 0:
356 
357  conn[0] = this->node_id(0);
358  conn[1] = this->node_id(4);
359  conn[2] = this->node_id(6);
360  conn[3] = this->node_id(7);
361 
362  return;
363 
364  // Linear sub-tet 1
365  case 1:
366 
367  conn[0] = this->node_id(4);
368  conn[1] = this->node_id(1);
369  conn[2] = this->node_id(5);
370  conn[3] = this->node_id(8);
371 
372  return;
373 
374  // Linear sub-tet 2
375  case 2:
376 
377  conn[0] = this->node_id(5);
378  conn[1] = this->node_id(2);
379  conn[2] = this->node_id(6);
380  conn[3] = this->node_id(9);
381 
382  return;
383 
384  // Linear sub-tet 3
385  case 3:
386 
387  conn[0] = this->node_id(7);
388  conn[1] = this->node_id(8);
389  conn[2] = this->node_id(9);
390  conn[3] = this->node_id(3);
391 
392  return;
393 
394  // Linear sub-tet 4
395  case 4:
396 
397  conn[0] = this->node_id(4);
398  conn[1] = this->node_id(8);
399  conn[2] = this->node_id(6);
400  conn[3] = this->node_id(7);
401 
402  return;
403 
404  // Linear sub-tet 5
405  case 5:
406 
407  conn[0] = this->node_id(4);
408  conn[1] = this->node_id(5);
409  conn[2] = this->node_id(6);
410  conn[3] = this->node_id(8);
411 
412  return;
413 
414  // Linear sub-tet 6
415  case 6:
416 
417  conn[0] = this->node_id(5);
418  conn[1] = this->node_id(9);
419  conn[2] = this->node_id(6);
420  conn[3] = this->node_id(8);
421 
422  return;
423 
424  // Linear sub-tet 7
425  case 7:
426 
427  conn[0] = this->node_id(7);
428  conn[1] = this->node_id(6);
429  conn[2] = this->node_id(9);
430  conn[3] = this->node_id(8);
431 
432  return;
433 
434 
435  default:
436 
437  libmesh_error_msg("Invalid sc = " << sc);
438  }
439  */
440  }
441 
442  default:
443  libmesh_error_msg("Unsupported IO package " << iop);
444  }
445 }
446 
447 
448 
449 const unsigned short int Tet10::_second_order_vertex_child_number[10] =
450  {
451  99,99,99,99, // Vertices
452  0,1,0,0,1,2 // Edges
453  };
454 
455 
456 
457 const unsigned short int Tet10::_second_order_vertex_child_index[10] =
458  {
459  99,99,99,99, // Vertices
460  1,2,2,3,3,3 // Edges
461  };
462 
463 
464 
465 std::pair<unsigned short int, unsigned short int>
466 Tet10::second_order_child_vertex (const unsigned int n) const
467 {
468  libmesh_assert_greater_equal (n, this->n_vertices());
469  libmesh_assert_less (n, this->n_nodes());
470  return std::pair<unsigned short int, unsigned short int>
471  (_second_order_vertex_child_number[n],
472  _second_order_vertex_child_index[n]);
473 }
474 
475 
476 
477 unsigned short int Tet10::second_order_adjacent_vertex (const unsigned int n,
478  const unsigned int v) const
479 {
480  libmesh_assert_greater_equal (n, this->n_vertices());
481  libmesh_assert_less (n, this->n_nodes());
482  libmesh_assert_less (v, 2);
483  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
484 }
485 
486 
487 
488 const unsigned short int Tet10::_second_order_adjacent_vertices[6][2] =
489  {
490  {0, 1}, // vertices adjacent to node 4
491  {1, 2}, // vertices adjacent to node 5
492  {0, 2}, // vertices adjacent to node 6
493  {0, 3}, // vertices adjacent to node 7
494  {1, 3}, // vertices adjacent to node 8
495  {2, 3} // vertices adjacent to node 9
496  };
497 
498 
499 
500 
501 
502 #ifdef LIBMESH_ENABLE_AMR
503 
504 const float Tet10::_embedding_matrix[8][10][10] =
505  {
506  // embedding matrix for child 0
507  {
508  // 0 1 2 3 4 5 6 7 8 9
509  { 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
510  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 1
511  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
512  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
513  { 0.375,-0.125, 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 4
514  { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 5
515  { 0.375, 0.,-0.125, 0., 0., 0., 0.75, 0., 0., 0.}, // 6
516  { 0.375, 0., 0.,-0.125, 0., 0., 0., 0.75, 0., 0.}, // 7
517  { 0.,-0.125, 0.,-0.125, 0.5, 0., 0., 0.5, 0.25, 0.}, // 8
518  { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25} // 9
519  },
520 
521  // embedding matrix for child 1
522  {
523  // 0 1 2 3 4 5 6 7 8 9
524  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
525  { 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
526  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 2
527  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
528  {-0.125, 0.375, 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 4
529  { 0., 0.375,-0.125, 0., 0., 0.75, 0., 0., 0., 0.}, // 5
530  {-0.125, 0.,-0.125, 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 6
531  {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 7
532  { 0., 0.375, 0.,-0.125, 0., 0., 0., 0., 0.75, 0.}, // 8
533  { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25} // 9
534  },
535 
536  // embedding matrix for child 2
537  {
538  // 0 1 2 3 4 5 6 7 8 9
539  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
540  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
541  { 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, // 2
542  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 3
543  {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 4
544  { 0.,-0.125, 0.375, 0., 0., 0.75, 0., 0., 0., 0.}, // 5
545  {-0.125, 0., 0.375, 0., 0., 0., 0.75, 0., 0., 0.}, // 6
546  {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 7
547  { 0.,-0.125, 0.,-0.125, 0., 0.5, 0., 0., 0.25, 0.5}, // 8
548  { 0., 0., 0.375,-0.125, 0., 0., 0., 0., 0., 0.75} // 9
549  },
550 
551  // embedding matrix for child 3
552  {
553  // 0 1 2 3 4 5 6 7 8 9
554  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 0
555  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
556  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
557  { 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, // 3
558  {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 4
559  { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5}, // 5
560  {-0.125, 0.,-0.125, 0., 0., 0., 0.25, 0.5, 0., 0.5}, // 6
561  {-0.125, 0., 0., 0.375, 0., 0., 0., 0.75, 0., 0.}, // 7
562  { 0.,-0.125, 0., 0.375, 0., 0., 0., 0., 0.75, 0.}, // 8
563  { 0., 0.,-0.125, 0.375, 0., 0., 0., 0., 0., 0.75} // 9
564  },
565 
566  // embedding matrix for child 4
567  {
568  // 0 1 2 3 4 5 6 7 8 9
569  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
570  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
571  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
572  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
573  {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 4
574  {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 5
575  { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 6
576  { 0.,-0.125, 0.,-0.125, 0.5, 0., 0., 0.5, 0.25, 0.}, // 7
577  {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 8
578  { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25} // 9
579  },
580 
581  // embedding matrix for child 5
582  {
583  // 0 1 2 3 4 5 6 7 8 9
584  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
585  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
586  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
587  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
588  {-0.125, 0.,-0.125, 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 4
589  {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 5
590  { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 6
591  {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 7
592  { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25}, // 8
593  {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25} // 9
594  },
595 
596  // embedding matrix for child 6
597  {
598  // 0 1 2 3 4 5 6 7 8 9
599  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
600  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
601  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
602  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
603  {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 4
604  { 0.,-0.125, 0.,-0.125, 0., 0.5, 0., 0., 0.25, 0.5}, // 5
605  {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 6
606  {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 7
607  { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25}, // 8
608  { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5} // 9
609  },
610 
611  // embedding matrix for child 7
612  {
613  // 0 1 2 3 4 5 6 7 8 9
614  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
615  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
616  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
617  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
618  {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 4
619  { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5}, // 5
620  {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 6
621  { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25}, // 7
622  {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 8
623  {-0.125, 0.,-0.125, 0., 0., 0., 0.25, 0.5, 0., 0.5} // 9
624  }
625  };
626 
627 
628 
629 float Tet10::embedding_matrix (const unsigned int i,
630  const unsigned int j,
631  const unsigned int k) const
632 {
633  // Choose an optimal diagonal, if one has not already been selected
634  this->choose_diagonal();
635 
636  // Permuted j and k indices
637  unsigned int
638  jp=j,
639  kp=k;
640 
641  if ((i>3) && (this->_diagonal_selection!=DIAG_02_13))
642  {
643  // Just the enum value cast to an unsigned int...
644  const unsigned ds = static_cast<unsigned int>(this->_diagonal_selection); // == 1 or 2
645 
646  // Instead of doing a lot of arithmetic, use these
647  // straightforward arrays for the permutations. Note that 3 ->
648  // 3, and the first array consists of "forward" permutations of
649  // the sets {0,1,2}, {4,5,6}, and {7,8,9} while the second array
650  // consists of "reverse" permutations of the same sets.
651  const unsigned int perms[2][10] =
652  {
653  {1, 2, 0, 3, 5, 6, 4, 8, 9, 7},
654  {2, 0, 1, 3, 6, 4, 5, 9, 7, 8}
655  };
656 
657  // Permute j
658  jp = perms[ds-1][j];
659  // if (jp<3)
660  // jp = (jp+ds)%3;
661  // else if (jp>3)
662  // jp = (jp-1+ds)%3 + 1 + 3*((jp-1)/3);
663 
664  // Permute k
665  kp = perms[ds-1][k];
666  // if (kp<3)
667  // kp = (kp+ds)%3;
668  // else if (kp>3)
669  // kp = (kp-1+ds)%3 + 1 + 3*((kp-1)/3);
670  }
671 
672  // Debugging:
673  // libMesh::err << "Selected diagonal " << _diagonal_selection << std::endl;
674  // libMesh::err << "j=" << j << std::endl;
675  // libMesh::err << "k=" << k << std::endl;
676  // libMesh::err << "jp=" << jp << std::endl;
677  // libMesh::err << "kp=" << kp << std::endl;
678 
679  // Call embedding matrix with permuted indices
680  return this->_embedding_matrix[i][jp][kp];
681 }
682 
683 #endif // #ifdef LIBMESH_ENABLE_AMR
684 
685 
686 
687 Real Tet10::volume () const
688 {
689  // Make copies of our points. It makes the subsequent calculations a bit
690  // shorter and avoids dereferencing the same pointer multiple times.
691  Point
692  x0 = point(0), x1 = point(1), x2 = point(2), x3 = point(3), x4 = point(4),
693  x5 = point(5), x6 = point(6), x7 = point(7), x8 = point(8), x9 = point(9);
694 
695  // The constant components of the dx/dxi vector, linear in xi, eta, zeta.
696  // These were copied directly from the output of a Python script.
697  Point dx_dxi[4] =
698  {
699  -3*x0 - x1 + 4*x4, // constant
700  4*x0 - 4*x4 - 4*x7 + 4*x8, // zeta
701  4*x0 - 4*x4 + 4*x5 - 4*x6, // eta
702  4*x0 + 4*x1 - 8*x4 // xi
703  };
704 
705  // The constant components of the dx/deta vector, linear in xi, eta, zeta.
706  // These were copied directly from the output of a Python script.
707  Point dx_deta[4] =
708  {
709  -3*x0 - x2 + 4*x6, // constant
710  4*x0 - 4*x6 - 4*x7 + 4*x9, // zeta
711  4*x0 + 4*x2 - 8*x6, // eta
712  4*x0 - 4*x4 + 4*x5 - 4*x6 // xi
713  };
714 
715  // The constant components of the dx/dzeta vector, linear in xi, eta, zeta.
716  // These were copied directly from the output of a Python script.
717  Point dx_dzeta[4] =
718  {
719  -3*x0 - x3 + 4*x7, // constant
720  4*x0 + 4*x3 - 8*x7, // zeta
721  4*x0 - 4*x6 - 4*x7 + 4*x9, // eta
722  4*x0 - 4*x4 - 4*x7 + 4*x8 // xi
723  };
724 
725  // 2x2x2 conical quadrature rule. Note: there is also a five point
726  // rule for tets with a negative weight which would be cheaper, but
727  // we'll use this one to preclude any possible issues with
728  // cancellation error.
729  const int N = 8;
730  static const Real w[N] =
731  {
732  3.6979856358852914509238091810505e-02L,
733  1.6027040598476613723156741868689e-02L,
734  2.1157006454524061178256145400082e-02L,
735  9.1694299214797439226823542540576e-03L,
736  3.6979856358852914509238091810505e-02L,
737  1.6027040598476613723156741868689e-02L,
738  2.1157006454524061178256145400082e-02L,
739  9.1694299214797439226823542540576e-03L
740  };
741 
742  static const Real xi[N] =
743  {
744  1.2251482265544137786674043037115e-01L,
745  5.4415184401122528879992623629551e-01L,
746  1.2251482265544137786674043037115e-01L,
747  5.4415184401122528879992623629551e-01L,
748  1.2251482265544137786674043037115e-01L,
749  5.4415184401122528879992623629551e-01L,
750  1.2251482265544137786674043037115e-01L,
751  5.4415184401122528879992623629551e-01L
752  };
753 
754  static const Real eta[N] =
755  {
756  1.3605497680284601717109468420738e-01L,
757  7.0679724159396903069267439165167e-02L,
758  5.6593316507280088053551297149570e-01L,
759  2.9399880063162286589079157179842e-01L,
760  1.3605497680284601717109468420738e-01L,
761  7.0679724159396903069267439165167e-02L,
762  5.6593316507280088053551297149570e-01L,
763  2.9399880063162286589079157179842e-01L
764  };
765 
766  static const Real zeta[N] =
767  {
768  1.5668263733681830907933725249176e-01L,
769  8.1395667014670255076709592007207e-02L,
770  6.5838687060044409936029672711329e-02L,
771  3.4202793236766414300604458388142e-02L,
772  5.8474756320489429588282763292971e-01L,
773  3.0377276481470755305409673253211e-01L,
774  2.4571332521171333166171692542182e-01L,
775  1.2764656212038543100867773351792e-01L
776  };
777 
778  Real vol = 0.;
779  for (int q=0; q<N; ++q)
780  {
781  // Compute dx_dxi, dx_deta, dx_dzeta at the current quadrature point.
782  Point
783  dx_dxi_q = dx_dxi[0] + zeta[q]*dx_dxi[1] + eta[q]*dx_dxi[2] + xi[q]*dx_dxi[3],
784  dx_deta_q = dx_deta[0] + zeta[q]*dx_deta[1] + eta[q]*dx_deta[2] + xi[q]*dx_deta[3],
785  dx_dzeta_q = dx_dzeta[0] + zeta[q]*dx_dzeta[1] + eta[q]*dx_dzeta[2] + xi[q]*dx_dzeta[3];
786 
787  // Compute scalar triple product, multiply by weight, and accumulate volume.
788  vol += w[q] * triple_product(dx_dxi_q, dx_deta_q, dx_dzeta_q);
789  }
790 
791  return vol;
792 }
793 
794 
795 
796 } // namespace libMesh
unsigned short int side
Definition: xdr_io.C:49
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)
T triple_product(const TypeVector< T > &a, const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1051
const dof_id_type n_nodes
Definition: tecplot_io.C:67
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.