face_tri6.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 // C++ includes
19 
20 // Local includes
21 #include "libmesh/side.h"
22 #include "libmesh/edge_edge3.h"
23 #include "libmesh/face_tri6.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 
31 // ------------------------------------------------------------
32 // Tri6 class static member initializations
33 const unsigned int Tri6::side_nodes_map[3][3] =
34  {
35  {0, 1, 3}, // Side 0
36  {1, 2, 4}, // Side 1
37  {2, 0, 5} // Side 2
38  };
39 
40 
41 #ifdef LIBMESH_ENABLE_AMR
42 
43 const float Tri6::_embedding_matrix[4][6][6] =
44  {
45  // embedding matrix for child 0
46  {
47  // 0 1 2 3 4 5
48  { 1.0, 0.0, 0.0, 0.0, 0.0, 0.0}, // 0
49  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 1
50  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}, // 2
51  {.375, -.125, 0.0, .75, 0.0, 0.0}, // 3
52  { 0.0, -.125, -.125, 0.5, .25, 0.5}, // 4
53  {.375, 0.0, -.125, 0.0, 0.0, .75} // 5
54  },
55 
56  // embedding matrix for child 1
57  {
58  // 0 1 2 3 4 5
59  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 0
60  { 0.0, 1.0, 0.0, 0.0, 0.0, 0.0}, // 1
61  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 2
62  {-.125, .375, 0.0, .75, 0.0, 0.0}, // 3
63  { 0.0, .375, -.125, 0.0, .75, 0.0}, // 4
64  {-.125, 0.0, -.125, 0.5, 0.5, .25} // 5
65  },
66 
67  // embedding matrix for child 2
68  {
69  // 0 1 2 3 4 5
70  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}, // 0
71  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 1
72  { 0.0, 0.0, 1.0, 0.0, 0.0, 0.0}, // 2
73  {-.125, -.125, 0.0, .25, 0.5, 0.5}, // 3
74  { 0.0, -.125, .375, 0.0, .75, 0.0}, // 4
75  {-.125, 0.0, .375, 0.0, 0.0, .75} // 5
76  },
77 
78  // embedding matrix for child 3
79  {
80  // 0 1 2 3 4 5
81  { 0.0, 0.0, 0.0, 1.0, 0.0, 0.0}, // 0
82  { 0.0, 0.0, 0.0, 0.0, 1.0, 0.0}, // 1
83  { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0}, // 2
84  {-.125, 0.0, -.125, 0.5, 0.5, .25}, // 3
85  {-.125, -.125, 0.0, .25, 0.5, 0.5}, // 4
86  { 0.0, -.125, -.125, 0.5, .25, 0.5} // 5
87  }
88  };
89 
90 #endif
91 
92 
93 
94 // ------------------------------------------------------------
95 // Tri6 class member functions
96 
97 bool Tri6::is_vertex(const unsigned int i) const
98 {
99  if (i < 3)
100  return true;
101  return false;
102 }
103 
104 bool Tri6::is_edge(const unsigned int i) const
105 {
106  if (i < 3)
107  return false;
108  return true;
109 }
110 
111 bool Tri6::is_face(const unsigned int) const
112 {
113  return false;
114 }
115 
116 bool Tri6::is_node_on_side(const unsigned int n,
117  const unsigned int s) const
118 {
119  libmesh_assert_less (s, n_sides());
120  for (unsigned int i = 0; i != 3; ++i)
121  if (side_nodes_map[s][i] == n)
122  return true;
123  return false;
124 }
125 
126 
127 
129 {
130  // Make sure edges are straight
131  if (!this->point(3).relative_fuzzy_equals
132  ((this->point(0) + this->point(1))/2.))
133  return false;
134  if (!this->point(4).relative_fuzzy_equals
135  ((this->point(1) + this->point(2))/2.))
136  return false;
137  if (!this->point(5).relative_fuzzy_equals
138  ((this->point(2) + this->point(0))/2.))
139  return false;
140 
141  return true;
142 }
143 
144 
145 
146 dof_id_type Tri6::key (const unsigned int s) const
147 {
148  libmesh_assert_less (s, this->n_sides());
149 
150  switch (s)
151  {
152  case 0:
153 
154  return
155  this->compute_key (this->node_id(3));
156 
157  case 1:
158 
159  return
160  this->compute_key (this->node_id(4));
161 
162  case 2:
163 
164  return
165  this->compute_key (this->node_id(5));
166 
167  default:
168  libmesh_error_msg("Invalid side s = " << s);
169  }
170 
171  libmesh_error_msg("We'll never get here!");
172  return 0;
173 }
174 
175 
176 
177 UniquePtr<Elem> Tri6::build_side_ptr (const unsigned int i,
178  bool proxy)
179 {
180  libmesh_assert_less (i, this->n_sides());
181 
182  if (proxy)
183  return UniquePtr<Elem>(new Side<Edge3,Tri6>(this,i));
184 
185  else
186  {
187  Elem * edge = new Edge3;
188  edge->subdomain_id() = this->subdomain_id();
189 
190  // Set the nodes
191  for (unsigned n=0; n<edge->n_nodes(); ++n)
192  edge->set_node(n) = this->node_ptr(Tri6::side_nodes_map[i][n]);
193 
194  return UniquePtr<Elem>(edge);
195  }
196 
197  libmesh_error_msg("We'll never get here!");
198  return UniquePtr<Elem>();
199 }
200 
201 
202 void Tri6::connectivity(const unsigned int sf,
203  const IOPackage iop,
204  std::vector<dof_id_type> & conn) const
205 {
206  libmesh_assert_less (sf, 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(4);
214  switch(sf)
215  {
216  case 0:
217  // linear sub-triangle 0
218  conn[0] = this->node_id(0)+1;
219  conn[1] = this->node_id(3)+1;
220  conn[2] = this->node_id(5)+1;
221  conn[3] = this->node_id(5)+1;
222 
223  return;
224 
225  case 1:
226  // linear sub-triangle 1
227  conn[0] = this->node_id(3)+1;
228  conn[1] = this->node_id(1)+1;
229  conn[2] = this->node_id(4)+1;
230  conn[3] = this->node_id(4)+1;
231 
232  return;
233 
234  case 2:
235  // linear sub-triangle 2
236  conn[0] = this->node_id(5)+1;
237  conn[1] = this->node_id(4)+1;
238  conn[2] = this->node_id(2)+1;
239  conn[3] = this->node_id(2)+1;
240 
241  return;
242 
243  case 3:
244  // linear sub-triangle 3
245  conn[0] = this->node_id(3)+1;
246  conn[1] = this->node_id(4)+1;
247  conn[2] = this->node_id(5)+1;
248  conn[3] = this->node_id(5)+1;
249 
250  return;
251 
252  default:
253  libmesh_error_msg("Invalid sf = " << sf);
254  }
255  }
256 
257  case VTK:
258  {
259  // VTK_QUADRATIC_TRIANGLE has same numbering as libmesh TRI6
260  conn.resize(6);
261  conn[0] = this->node_id(0);
262  conn[1] = this->node_id(1);
263  conn[2] = this->node_id(2);
264  conn[3] = this->node_id(3);
265  conn[4] = this->node_id(4);
266  conn[5] = this->node_id(5);
267  return;
268 
269  // Used to write out linear sub-triangles for VTK...
270  /*
271  conn.resize(3);
272  switch(sf)
273  {
274  case 0:
275  // linear sub-triangle 0
276  conn[0] = this->node_id(0);
277  conn[1] = this->node_id(3);
278  conn[2] = this->node_id(5);
279 
280  return;
281 
282  case 1:
283  // linear sub-triangle 1
284  conn[0] = this->node_id(3);
285  conn[1] = this->node_id(1);
286  conn[2] = this->node_id(4);
287 
288  return;
289 
290  case 2:
291  // linear sub-triangle 2
292  conn[0] = this->node_id(5);
293  conn[1] = this->node_id(4);
294  conn[2] = this->node_id(2);
295 
296  return;
297 
298  case 3:
299  // linear sub-triangle 3
300  conn[0] = this->node_id(3);
301  conn[1] = this->node_id(4);
302  conn[2] = this->node_id(5);
303 
304  return;
305 
306  default:
307  libmesh_error_msg("Invalid sf = " << sf);
308  }
309  */
310  }
311 
312  default:
313  libmesh_error_msg("Unsupported IO package " << iop);
314  }
315 }
316 
317 
318 
320 {
321  // This might have curved edges, or might be a curved surface in
322  // 3-space, in which case the full bounding box can be larger than
323  // the bounding box of just the nodes.
324  //
325  //
326  // FIXME - I haven't yet proven the formula below to be correct for
327  // quadratics in 2D - RHS
328  Point pmin, pmax;
329 
330  for (unsigned d=0; d<LIBMESH_DIM; ++d)
331  {
332  Real center = this->point(0)(d);
333  for (unsigned int p=1; p != 6; ++p)
334  center += this->point(p)(d);
335  center /= 6;
336 
337  Real hd = std::abs(center - this->point(0)(d));
338  for (unsigned int p=1; p != 6; ++p)
339  hd = std::max(hd, std::abs(center - this->point(p)(d)));
340 
341  pmin(d) = center - hd;
342  pmax(d) = center + hd;
343  }
344 
345  return BoundingBox(pmin, pmax);
346 }
347 
348 
349 
351 {
352  // Make copies of our points. It makes the subsequent calculations a bit
353  // shorter and avoids dereferencing the same pointer multiple times.
354  Point
355  x0 = point(0), x1 = point(1), x2 = point(2),
356  x3 = point(3), x4 = point(4), x5 = point(5);
357 
358  // Construct constant data vectors.
359  // \vec{x}_{\xi} = \vec{a1}*xi + \vec{b1}*eta + \vec{c1}
360  // \vec{x}_{\eta} = \vec{a2}*xi + \vec{b2}*eta + \vec{c2}
361  Point
362  a1 = 4*x0 + 4*x1 - 8*x3,
363  b1 = 4*x0 - 4*x3 + 4*x4 - 4*x5, /*=a2*/
364  c1 = -3*x0 - 1*x1 + 4*x3,
365  b2 = 4*x0 + 4*x2 - 8*x5,
366  c2 = -3*x0 - 1*x2 + 4*x5;
367 
368  // If a1 == b1 == a2 == b2 == 0, this is a TRI6 with straight sides,
369  // and we can use the TRI3 formula to compute the volume.
370  if (a1.relative_fuzzy_equals(Point(0,0,0)) &&
371  b1.relative_fuzzy_equals(Point(0,0,0)) &&
372  b2.relative_fuzzy_equals(Point(0,0,0)))
373  return 0.5 * cross_norm(c1, c2);
374 
375  // 7-point rule, exact for quintics.
376  const unsigned int N = 7;
377 
378  // Parameters of the quadrature rule
379  const static Real
380  w1 = Real(31)/480 + std::sqrt(15.0L)/2400,
381  w2 = Real(31)/480 - std::sqrt(15.0L)/2400,
382  q1 = Real(2)/7 + std::sqrt(15.0L)/21,
383  q2 = Real(2)/7 - std::sqrt(15.0L)/21;
384 
385  const static Real xi[N] = {Real(1)/3, q1, q1, 1-2*q1, q2, q2, 1-2*q2};
386  const static Real eta[N] = {Real(1)/3, q1, 1-2*q1, q1, q2, 1-2*q2, q2};
387  const static Real wts[N] = {Real(9)/80, w1, w1, w1, w2, w2, w2};
388 
389  // Approximate the area with quadrature
390  Real vol=0.;
391  for (unsigned int q=0; q<N; ++q)
392  vol += wts[q] * cross_norm(xi[q]*a1 + eta[q]*b1 + c1,
393  xi[q]*b1 + eta[q]*b2 + c2);
394 
395  return vol;
396 }
397 
398 
399 
400 unsigned short int Tri6::second_order_adjacent_vertex (const unsigned int n,
401  const unsigned int v) const
402 {
403  libmesh_assert_greater_equal (n, this->n_vertices());
404  libmesh_assert_less (n, this->n_nodes());
405  libmesh_assert_less (v, 2);
406  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
407 }
408 
409 
410 
411 const unsigned short int Tri6::_second_order_adjacent_vertices[3][2] =
412  {
413  {0, 1}, // vertices adjacent to node 3
414  {1, 2}, // vertices adjacent to node 4
415  {0, 2} // vertices adjacent to node 5
416  };
417 
418 
419 
420 std::pair<unsigned short int, unsigned short int>
421 Tri6::second_order_child_vertex (const unsigned int n) const
422 {
423  libmesh_assert_greater_equal (n, this->n_vertices());
424  libmesh_assert_less (n, this->n_nodes());
425  return std::pair<unsigned short int, unsigned short int>
428 }
429 
430 
431 
432 const unsigned short int Tri6::_second_order_vertex_child_number[6] =
433  {
434  99,99,99, // Vertices
435  0,1,0 // Edges
436  };
437 
438 
439 
440 const unsigned short int Tri6::_second_order_vertex_child_index[6] =
441  {
442  99,99,99, // Vertices
443  1,2,2 // Edges
444  };
445 
446 } // namespace libMesh
double abs(double a)
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1723
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1068
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy) libmesh_override
Definition: face_tri6.C:177
virtual unsigned int n_sub_elem() const libmesh_override
Definition: face_tri6.h:73
The base class for all geometric element types.
Definition: elem.h:86
virtual bool is_edge(const unsigned int i) const libmesh_override
Definition: face_tri6.C:104
long double max(long double a, double b)
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const libmesh_override
Definition: face_tri6.C:400
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const libmesh_override
Definition: face_tri6.C:202
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
virtual bool has_affine_map() const libmesh_override
Definition: face_tri6.C:128
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1658
virtual unsigned int n_nodes() const libmesh_override
Definition: face_tri6.h:68
virtual bool is_face(const unsigned int i) const libmesh_override
Definition: face_tri6.C:111
virtual unsigned int n_vertices() const libmesh_override
Definition: face_tri.h:91
subdomain_id_type subdomain_id() const
Definition: elem.h:1733
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const libmesh_override
Definition: face_tri6.C:116
A 1D geometric element with 3 nodes.
Definition: edge_edge3.h:43
virtual BoundingBox loose_bounding_box() const libmesh_override
Definition: face_tri6.C:319
Proxy class for efficiently representing an Elem&#39;s side.
Definition: side.h:48
static const unsigned short int _second_order_vertex_child_number[6]
Definition: face_tri6.h:219
virtual unsigned int n_sides() const libmesh_override
Definition: face_tri.h:86
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual dof_id_type key() const libmesh_override
Definition: face_tri.C:62
const Point & point(const unsigned int i) const
Definition: elem.h:1595
static const unsigned short int _second_order_adjacent_vertices[3][2]
Definition: face_tri6.h:214
virtual std::pair< unsigned short int, unsigned short int > second_order_child_vertex(const unsigned int n) const libmesh_override
Definition: face_tri6.C:421
static const unsigned int side_nodes_map[3][3]
Definition: face_tri6.h:165
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1617
static const unsigned short int _second_order_vertex_child_index[6]
Definition: face_tri6.h:224
static dof_id_type compute_key(dof_id_type n0)
Definition: elem.h:2336
virtual Real volume() const libmesh_override
Definition: face_tri6.C:350
static const float _embedding_matrix[4][6][6]
Definition: face_tri6.h:202
A geometric point in (x,y,z) space.
Definition: point.h:38
bool relative_fuzzy_equals(const TypeVector< T > &rhs, Real tol=TOLERANCE) const
Definition: type_vector.h:963
virtual bool is_vertex(const unsigned int i) const libmesh_override
Definition: face_tri6.C:97
uint8_t dof_id_type
Definition: id_types.h:64