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