cell_tet10.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_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}, // Side 0
45  {1, 2, 5}, // Side 1
46  {0, 2, 6}, // Side 2
47  {0, 3, 7}, // Side 3
48  {1, 3, 8}, // Side 4
49  {2, 3, 9} // Side 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 
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 AutoPtr<Elem> Tet10::build_side (const unsigned int i,
157  bool proxy) const
158 {
159  libmesh_assert_less (i, this->n_sides());
160 
161  if (proxy)
162  {
163  AutoPtr<Elem> ap(new Side<Tri6,Tet10>(this,i));
164  return ap;
165  }
166 
167  else
168  {
169  AutoPtr<Elem> face(new Tri6);
170  face->subdomain_id() = this->subdomain_id();
171 
172  switch (i)
173  {
174  case 0:
175  {
176  face->set_node(0) = this->get_node(0);
177  face->set_node(1) = this->get_node(2);
178  face->set_node(2) = this->get_node(1);
179  face->set_node(3) = this->get_node(6);
180  face->set_node(4) = this->get_node(5);
181  face->set_node(5) = this->get_node(4);
182 
183  return face;
184  }
185  case 1:
186  {
187  face->set_node(0) = this->get_node(0);
188  face->set_node(1) = this->get_node(1);
189  face->set_node(2) = this->get_node(3);
190  face->set_node(3) = this->get_node(4);
191  face->set_node(4) = this->get_node(8);
192  face->set_node(5) = this->get_node(7);
193 
194  return face;
195  }
196  case 2:
197  {
198  face->set_node(0) = this->get_node(1);
199  face->set_node(1) = this->get_node(2);
200  face->set_node(2) = this->get_node(3);
201  face->set_node(3) = this->get_node(5);
202  face->set_node(4) = this->get_node(9);
203  face->set_node(5) = this->get_node(8);
204 
205  return face;
206  }
207  case 3:
208  {
209  face->set_node(0) = this->get_node(2);
210  face->set_node(1) = this->get_node(0);
211  face->set_node(2) = this->get_node(3);
212  face->set_node(3) = this->get_node(6);
213  face->set_node(4) = this->get_node(7);
214  face->set_node(5) = this->get_node(9);
215 
216  return face;
217  }
218  default:
219  {
220  libmesh_error();
221  }
222  }
223  }
224 
225 
226  // We'll never get here.
227  libmesh_error();
228  AutoPtr<Elem> ap(NULL); return ap;
229 }
230 
231 
232 
233 AutoPtr<Elem> Tet10::build_edge (const unsigned int i) const
234 {
235  libmesh_assert_less (i, this->n_edges());
236 
237  return AutoPtr<Elem>(new SideEdge<Edge3,Tet10>(this,i));
238 }
239 
240 
241 
242 void Tet10::connectivity(const unsigned int sc,
243  const IOPackage iop,
244  std::vector<dof_id_type>& conn) const
245 {
247  libmesh_assert_less (sc, this->n_sub_elem());
248  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
249 
250  switch (iop)
251  {
252  case TECPLOT:
253  {
254  conn.resize(8);
255  switch (sc)
256  {
257 
258 
259  // Linear sub-tet 0
260  case 0:
261 
262  conn[0] = this->node(0)+1;
263  conn[1] = this->node(4)+1;
264  conn[2] = this->node(6)+1;
265  conn[3] = this->node(6)+1;
266  conn[4] = this->node(7)+1;
267  conn[5] = this->node(7)+1;
268  conn[6] = this->node(7)+1;
269  conn[7] = this->node(7)+1;
270 
271  return;
272 
273  // Linear sub-tet 1
274  case 1:
275 
276  conn[0] = this->node(4)+1;
277  conn[1] = this->node(1)+1;
278  conn[2] = this->node(5)+1;
279  conn[3] = this->node(5)+1;
280  conn[4] = this->node(8)+1;
281  conn[5] = this->node(8)+1;
282  conn[6] = this->node(8)+1;
283  conn[7] = this->node(8)+1;
284 
285  return;
286 
287  // Linear sub-tet 2
288  case 2:
289 
290  conn[0] = this->node(5)+1;
291  conn[1] = this->node(2)+1;
292  conn[2] = this->node(6)+1;
293  conn[3] = this->node(6)+1;
294  conn[4] = this->node(9)+1;
295  conn[5] = this->node(9)+1;
296  conn[6] = this->node(9)+1;
297  conn[7] = this->node(9)+1;
298 
299  return;
300 
301  // Linear sub-tet 3
302  case 3:
303 
304  conn[0] = this->node(7)+1;
305  conn[1] = this->node(8)+1;
306  conn[2] = this->node(9)+1;
307  conn[3] = this->node(9)+1;
308  conn[4] = this->node(3)+1;
309  conn[5] = this->node(3)+1;
310  conn[6] = this->node(3)+1;
311  conn[7] = this->node(3)+1;
312 
313  return;
314 
315  // Linear sub-tet 4
316  case 4:
317 
318  conn[0] = this->node(4)+1;
319  conn[1] = this->node(8)+1;
320  conn[2] = this->node(6)+1;
321  conn[3] = this->node(6)+1;
322  conn[4] = this->node(7)+1;
323  conn[5] = this->node(7)+1;
324  conn[6] = this->node(7)+1;
325  conn[7] = this->node(7)+1;
326 
327  return;
328 
329  // Linear sub-tet 5
330  case 5:
331 
332  conn[0] = this->node(4)+1;
333  conn[1] = this->node(5)+1;
334  conn[2] = this->node(6)+1;
335  conn[3] = this->node(6)+1;
336  conn[4] = this->node(8)+1;
337  conn[5] = this->node(8)+1;
338  conn[6] = this->node(8)+1;
339  conn[7] = this->node(8)+1;
340 
341  return;
342 
343  // Linear sub-tet 6
344  case 6:
345 
346  conn[0] = this->node(5)+1;
347  conn[1] = this->node(9)+1;
348  conn[2] = this->node(6)+1;
349  conn[3] = this->node(6)+1;
350  conn[4] = this->node(8)+1;
351  conn[5] = this->node(8)+1;
352  conn[6] = this->node(8)+1;
353  conn[7] = this->node(8)+1;
354 
355  return;
356 
357  // Linear sub-tet 7
358  case 7:
359 
360  conn[0] = this->node(7)+1;
361  conn[1] = this->node(6)+1;
362  conn[2] = this->node(9)+1;
363  conn[3] = this->node(9)+1;
364  conn[4] = this->node(8)+1;
365  conn[5] = this->node(8)+1;
366  conn[6] = this->node(8)+1;
367  conn[7] = this->node(8)+1;
368 
369  return;
370 
371 
372  default:
373 
374  libmesh_error();
375  }
376  }
377 
378  case VTK:
379  {
380  conn.resize(10);
381  conn[0] = this->node(0);
382  conn[1] = this->node(1);
383  conn[2] = this->node(2);
384  conn[3] = this->node(3);
385  conn[4] = this->node(4);
386  conn[5] = this->node(5);
387  conn[6] = this->node(6);
388  conn[7] = this->node(7);
389  conn[8] = this->node(8);
390  conn[9] = this->node(9);
391  return;
392  /*
393  conn.resize(4);
394  switch (sc)
395  {
396  // Linear sub-tet 0
397  case 0:
398 
399  conn[0] = this->node(0);
400  conn[1] = this->node(4);
401  conn[2] = this->node(6);
402  conn[3] = this->node(7);
403 
404  return;
405 
406  // Linear sub-tet 1
407  case 1:
408 
409  conn[0] = this->node(4);
410  conn[1] = this->node(1);
411  conn[2] = this->node(5);
412  conn[3] = this->node(8);
413 
414  return;
415 
416  // Linear sub-tet 2
417  case 2:
418 
419  conn[0] = this->node(5);
420  conn[1] = this->node(2);
421  conn[2] = this->node(6);
422  conn[3] = this->node(9);
423 
424  return;
425 
426  // Linear sub-tet 3
427  case 3:
428 
429  conn[0] = this->node(7);
430  conn[1] = this->node(8);
431  conn[2] = this->node(9);
432  conn[3] = this->node(3);
433 
434  return;
435 
436  // Linear sub-tet 4
437  case 4:
438 
439  conn[0] = this->node(4);
440  conn[1] = this->node(8);
441  conn[2] = this->node(6);
442  conn[3] = this->node(7);
443 
444  return;
445 
446  // Linear sub-tet 5
447  case 5:
448 
449  conn[0] = this->node(4);
450  conn[1] = this->node(5);
451  conn[2] = this->node(6);
452  conn[3] = this->node(8);
453 
454  return;
455 
456  // Linear sub-tet 6
457  case 6:
458 
459  conn[0] = this->node(5);
460  conn[1] = this->node(9);
461  conn[2] = this->node(6);
462  conn[3] = this->node(8);
463 
464  return;
465 
466  // Linear sub-tet 7
467  case 7:
468 
469  conn[0] = this->node(7);
470  conn[1] = this->node(6);
471  conn[2] = this->node(9);
472  conn[3] = this->node(8);
473 
474  return;
475 
476 
477  default:
478 
479  libmesh_error();
480  }
481  */
482  }
483 
484  default:
485  libmesh_error();
486  }
487 
488  libmesh_error();
489 }
490 
491 
492 
493 const unsigned short int Tet10::_second_order_vertex_child_number[10] =
494 {
495  99,99,99,99, // Vertices
496  0,1,0,0,1,2 // Edges
497 };
498 
499 
500 
501 const unsigned short int Tet10::_second_order_vertex_child_index[10] =
502 {
503  99,99,99,99, // Vertices
504  1,2,2,3,3,3 // Edges
505 };
506 
507 
508 
509 std::pair<unsigned short int, unsigned short int>
510 Tet10::second_order_child_vertex (const unsigned int n) const
511 {
512  libmesh_assert_greater_equal (n, this->n_vertices());
513  libmesh_assert_less (n, this->n_nodes());
514  return std::pair<unsigned short int, unsigned short int>
517 }
518 
519 
520 
521 unsigned short int Tet10::second_order_adjacent_vertex (const unsigned int n,
522  const unsigned int v) const
523 {
524  libmesh_assert_greater_equal (n, this->n_vertices());
525  libmesh_assert_less (n, this->n_nodes());
526  libmesh_assert_less (v, 2);
527  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
528 }
529 
530 
531 
532 const unsigned short int Tet10::_second_order_adjacent_vertices[6][2] =
533 {
534  {0, 1}, // vertices adjacent to node 4
535  {1, 2}, // vertices adjacent to node 5
536  {0, 2}, // vertices adjacent to node 6
537  {0, 3}, // vertices adjacent to node 7
538  {1, 3}, // vertices adjacent to node 8
539  {2, 3} // vertices adjacent to node 9
540 };
541 
542 
543 
544 
545 
546 #ifdef LIBMESH_ENABLE_AMR
547 
548 const float Tet10::_embedding_matrix[8][10][10] =
549 {
550  // embedding matrix for child 0
551  {
552  // 0 1 2 3 4 5 6 7 8 9
553  { 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.}, // 0
554  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 1
555  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
556  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
557  { 0.375,-0.125, 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 4
558  { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 5
559  { 0.375, 0.,-0.125, 0., 0., 0., 0.75, 0., 0., 0.}, // 6
560  { 0.375, 0., 0.,-0.125, 0., 0., 0., 0.75, 0., 0.}, // 7
561  { 0.,-0.125, 0.,-0.125, 0.5, 0., 0., 0.5, 0.25, 0.}, // 8
562  { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25} // 9
563  },
564 
565  // embedding matrix for child 1
566  {
567  // 0 1 2 3 4 5 6 7 8 9
568  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
569  { 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.}, // 1
570  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 2
571  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
572  {-0.125, 0.375, 0., 0., 0.75, 0., 0., 0., 0., 0.}, // 4
573  { 0., 0.375,-0.125, 0., 0., 0.75, 0., 0., 0., 0.}, // 5
574  {-0.125, 0.,-0.125, 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 6
575  {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 7
576  { 0., 0.375, 0.,-0.125, 0., 0., 0., 0., 0.75, 0.}, // 8
577  { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25} // 9
578  },
579 
580  // embedding matrix for child 2
581  {
582  // 0 1 2 3 4 5 6 7 8 9
583  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
584  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
585  { 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.}, // 2
586  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 3
587  {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 4
588  { 0.,-0.125, 0.375, 0., 0., 0.75, 0., 0., 0., 0.}, // 5
589  {-0.125, 0., 0.375, 0., 0., 0., 0.75, 0., 0., 0.}, // 6
590  {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 7
591  { 0.,-0.125, 0.,-0.125, 0., 0.5, 0., 0., 0.25, 0.5}, // 8
592  { 0., 0., 0.375,-0.125, 0., 0., 0., 0., 0., 0.75} // 9
593  },
594 
595  // embedding matrix for child 3
596  {
597  // 0 1 2 3 4 5 6 7 8 9
598  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 0
599  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
600  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
601  { 0., 0., 0., 1., 0., 0., 0., 0., 0., 0.}, // 3
602  {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 4
603  { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5}, // 5
604  {-0.125, 0.,-0.125, 0., 0., 0., 0.25, 0.5, 0., 0.5}, // 6
605  {-0.125, 0., 0., 0.375, 0., 0., 0., 0.75, 0., 0.}, // 7
606  { 0.,-0.125, 0., 0.375, 0., 0., 0., 0., 0.75, 0.}, // 8
607  { 0., 0.,-0.125, 0.375, 0., 0., 0., 0., 0., 0.75} // 9
608  },
609 
610  // embedding matrix for child 4
611  {
612  // 0 1 2 3 4 5 6 7 8 9
613  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
614  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
615  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
616  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
617  {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 4
618  {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 5
619  { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 6
620  { 0.,-0.125, 0.,-0.125, 0.5, 0., 0., 0.5, 0.25, 0.}, // 7
621  {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 8
622  { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25} // 9
623  },
624 
625  // embedding matrix for child 5
626  {
627  // 0 1 2 3 4 5 6 7 8 9
628  { 0., 0., 0., 0., 1., 0., 0., 0., 0., 0.}, // 0
629  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
630  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 2
631  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
632  {-0.125, 0.,-0.125, 0., 0.5, 0.5, 0.25, 0., 0., 0.}, // 4
633  {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 5
634  { 0.,-0.125,-0.125, 0., 0.5, 0.25, 0.5, 0., 0., 0.}, // 6
635  {-0.125, 0., 0.,-0.125, 0.5, 0., 0., 0.25, 0.5, 0.}, // 7
636  { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25}, // 8
637  {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25} // 9
638  },
639 
640  // embedding matrix for child 6
641  {
642  // 0 1 2 3 4 5 6 7 8 9
643  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
644  { 0., 0., 0., 0., 0., 1., 0., 0., 0., 0.}, // 1
645  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
646  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 3
647  {-0.125,-0.125, 0., 0., 0.25, 0.5, 0.5, 0., 0., 0.}, // 4
648  { 0.,-0.125, 0.,-0.125, 0., 0.5, 0., 0., 0.25, 0.5}, // 5
649  {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 6
650  {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 7
651  { 0., 0.,-0.125,-0.125, 0., 0.5, 0., 0., 0.5, 0.25}, // 8
652  { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5} // 9
653  },
654 
655  // embedding matrix for child 7
656  {
657  // 0 1 2 3 4 5 6 7 8 9
658  { 0., 0., 0., 0., 0., 0., 1., 0., 0., 0.}, // 0
659  { 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.}, // 1
660  { 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.}, // 2
661  { 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.}, // 3
662  {-0.125,-0.125,-0.125,-0.125, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25}, // 4
663  { 0.,-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5}, // 5
664  {-0.125, 0., 0.,-0.125, 0., 0., 0.5, 0.25, 0., 0.5}, // 6
665  { 0., 0.,-0.125,-0.125, 0., 0., 0.5, 0.5, 0., 0.25}, // 7
666  {-0.125,-0.125, 0., 0., 0.25, 0., 0., 0.5, 0.5, 0.}, // 8
667  {-0.125, 0.,-0.125, 0., 0., 0., 0.25, 0.5, 0., 0.5} // 9
668  }
669 };
670 
671 
672 
673 
674 
675 
676 float Tet10::embedding_matrix (const unsigned int i,
677  const unsigned int j,
678  const unsigned int k) const
679 {
680  // Choose an optimal diagonal, if one has not already been selected
681  this->choose_diagonal();
682 
683  // Permuted j and k indices
684  unsigned int
685  jp=j,
686  kp=k;
687 
688  if ((i>3) && (this->_diagonal_selection!=DIAG_02_13))
689  {
690  // Just the enum value cast to an unsigned int...
691  const unsigned ds = static_cast<unsigned int>(this->_diagonal_selection); // == 1 or 2
692 
693  // Instead of doing a lot of arithmetic, use these
694  // straightforward arrays for the permutations. Note that 3 ->
695  // 3, and the first array consists of "forward" permutations of
696  // the sets {0,1,2}, {4,5,6}, and {7,8,9} while the second array
697  // consists of "reverse" permutations of the same sets.
698  const unsigned int perms[2][10] =
699  {
700  {1, 2, 0, 3, 5, 6, 4, 8, 9, 7},
701  {2, 0, 1, 3, 6, 4, 5, 9, 7, 8}
702  };
703 
704  // Permute j
705  jp = perms[ds-1][j];
706  // if (jp<3)
707  // jp = (jp+ds)%3;
708  // else if (jp>3)
709  // jp = (jp-1+ds)%3 + 1 + 3*((jp-1)/3);
710 
711  // Permute k
712  kp = perms[ds-1][k];
713  // if (kp<3)
714  // kp = (kp+ds)%3;
715  // else if (kp>3)
716  // kp = (kp-1+ds)%3 + 1 + 3*((kp-1)/3);
717  }
718 
719  // Debugging:
720  // libMesh::err << "Selected diagonal " << _diagonal_selection << std::endl;
721  // libMesh::err << "j=" << j << std::endl;
722  // libMesh::err << "k=" << k << std::endl;
723  // libMesh::err << "jp=" << jp << std::endl;
724  // libMesh::err << "kp=" << kp << std::endl;
725 
726  // Call embedding matrx with permuted indices
727  return this->_embedding_matrix[i][jp][kp];
728 }
729 
730 #endif // #ifdef LIBMESH_ENABLE_AMR
731 
732 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo