cell_inf_hex.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 // Local includes
19 #include "libmesh/libmesh_config.h"
20 
21 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
22 
23 
24 // C++ includes
25 #include <algorithm> // for std::min, std::max
26 
27 // Local includes cont'd
28 #include "libmesh/cell_inf_hex.h"
29 #include "libmesh/cell_inf_hex8.h"
30 #include "libmesh/face_quad4.h"
31 #include "libmesh/face_inf_quad4.h"
32 
33 namespace libMesh
34 {
35 
36 
37 
38 
39 // ------------------------------------------------------------
40 // InfHex class member functions
41 dof_id_type InfHex::key (const unsigned int s) const
42 {
43  libmesh_assert_less (s, this->n_sides());
44 
45  switch (s)
46  {
47  case 0: // the face at z = -1
48 
49  return
50  this->compute_key (this->node(0),
51  this->node(1),
52  this->node(2),
53  this->node(3));
54 
55  case 1: // the face at y = -1
56 
57  return
58  this->compute_key (this->node(0),
59  this->node(1),
60  this->node(5),
61  this->node(4));
62 
63  case 2: // the face at x = 1
64 
65  return
66  this->compute_key (this->node(1),
67  this->node(2),
68  this->node(6),
69  this->node(5));
70 
71  case 3: // the face at y = 1
72 
73  return
74  this->compute_key (this->node(2),
75  this->node(3),
76  this->node(7),
77  this->node(6));
78 
79 
80  case 4: // the face at x = -1
81 
82  return
83  this->compute_key (this->node(3),
84  this->node(0),
85  this->node(4),
86  this->node(7));
87  }
88 
89  // We'll never get here.
90  libmesh_error();
91  return 0;
92 }
93 
94 
95 
96 AutoPtr<Elem> InfHex::side (const unsigned int i) const
97 {
98  libmesh_assert_less (i, this->n_sides());
99 
100  /*
101  *Think of a unit cube: (-1,1) x (-1,1)x (-1,1),
102  * with (in general) the normals pointing outwards
103  */
104  switch (i)
105  {
106  case 0: // the face at z = -1
107  // the base, where the infinite element couples to conventional
108  // elements
109  {
110  Elem* face = new Quad4;
111  AutoPtr<Elem> ap_face(face);
112 
113  /*
114  * Oops, here we are, claiming the normal of the face
115  * elements point outwards -- and this is the exception:
116  * For the side built from the base face,
117  * the normal is pointing _into_ the element!
118  * Why is that? - In agreement with build_side(),
119  * which in turn _has_ to build the face in this
120  * way as to enable the cool way \p InfFE re-uses \p FE.
121  */
122  face->set_node(0) = this->get_node(0);
123  face->set_node(1) = this->get_node(1);
124  face->set_node(2) = this->get_node(2);
125  face->set_node(3) = this->get_node(3);
126 
127  return ap_face;
128  }
129 
130  case 1: // the face at y = -1
131  // this face connects to another infinite element
132  {
133  Elem* face = new InfQuad4;
134  AutoPtr<Elem> ap_face(face);
135 
136  face->set_node(0) = this->get_node(0);
137  face->set_node(1) = this->get_node(1);
138  face->set_node(2) = this->get_node(4);
139  face->set_node(3) = this->get_node(5);
140 
141  return ap_face;
142  }
143 
144  case 2: // the face at x = 1
145  // this face connects to another infinite element
146  {
147  Elem* face = new InfQuad4;
148  AutoPtr<Elem> ap_face(face);
149  //AutoPtr<Elem> face(new InfQuad4);
150 
151  face->set_node(0) = this->get_node(1);
152  face->set_node(1) = this->get_node(2);
153  face->set_node(2) = this->get_node(5);
154  face->set_node(3) = this->get_node(6);
155 
156  return ap_face;
157  }
158 
159  case 3: // the face at y = 1
160  // this face connects to another infinite element
161  {
162  Elem* face = new InfQuad4;
163  AutoPtr<Elem> ap_face(face);
164  //AutoPtr<Elem> face(new InfQuad4);
165 
166  face->set_node(0) = this->get_node(2);
167  face->set_node(1) = this->get_node(3);
168  face->set_node(2) = this->get_node(6);
169  face->set_node(3) = this->get_node(7);
170 
171  return ap_face;
172  }
173 
174  case 4: // the face at x = -1
175  // this face connects to another infinite element
176  {
177  Elem* face = new InfQuad4;
178  AutoPtr<Elem> ap_face(face);
179  //AutoPtr<Elem> face(new InfQuad4);
180 
181  face->set_node(0) = this->get_node(3);
182  face->set_node(1) = this->get_node(0);
183  face->set_node(2) = this->get_node(7);
184  face->set_node(3) = this->get_node(4);
185 
186  return ap_face;
187  }
188 
189  default:
190  {
191  libmesh_error();
192  AutoPtr<Elem> ap(NULL); return ap;
193  }
194  }
195 
196  // We'll never get here.
197  libmesh_error();
198  AutoPtr<Elem> ap(NULL); return ap;
199 }
200 
201 
202 
203 bool InfHex::is_child_on_side(const unsigned int c,
204  const unsigned int s) const
205 {
206  libmesh_assert_less (c, this->n_children());
207  libmesh_assert_less (s, this->n_sides());
208 
209  return (s == 0 || c+1 == s || c == s%4);
210 }
211 
212 
213 
214 bool InfHex::is_edge_on_side (const unsigned int e,
215  const unsigned int s) const
216 {
217  libmesh_assert_less (e, this->n_edges());
218  libmesh_assert_less (s, this->n_sides());
219 
220  return (is_node_on_side(InfHex8::edge_nodes_map[e][0],s) &&
222 }
223 
224 
225 
226 
228 {
229  switch (q)
230  {
231 
241  case DIAGONAL:
242  {
243  // Diagonal between node 0 and node 2
244  const Real d02 = this->length(0,2);
245 
246  // Diagonal between node 1 and node 3
247  const Real d13 = this->length(1,3);
248 
249  // Find the biggest and smallest diagonals
250  const Real min = std::min(d02, d13);
251  const Real max = std::max(d02, d13);
252 
253  libmesh_assert_not_equal_to (max, 0.0);
254 
255  return min / max;
256 
257  break;
258  }
259 
267  case TAPER:
268  {
269 
273  const Real d01 = this->length(0,1);
274  const Real d12 = this->length(1,2);
275  const Real d23 = this->length(2,3);
276  const Real d03 = this->length(0,3);
277 
278  std::vector<Real> edge_ratios(2);
279 
280  // Bottom
281  edge_ratios[8] = std::min(d01, d23) / std::max(d01, d23);
282  edge_ratios[9] = std::min(d03, d12) / std::max(d03, d12);
283 
284  return *(std::min_element(edge_ratios.begin(), edge_ratios.end())) ;
285 
286  break;
287  }
288 
289 
297  case STRETCH:
298  {
302  const Real sqrt3 = 1.73205080756888;
303 
307  const Real d02 = this->length(0,2);
308  const Real d13 = this->length(1,3);
309  const Real max_diag = std::max(d02, d13);
310 
311  libmesh_assert_not_equal_to ( max_diag, 0.0 );
312 
316  std::vector<Real> edges(4);
317  edges[0] = this->length(0,1);
318  edges[1] = this->length(1,2);
319  edges[2] = this->length(2,3);
320  edges[3] = this->length(0,3);
321 
322  const Real min_edge = *(std::min_element(edges.begin(), edges.end()));
323  return sqrt3 * min_edge / max_diag ;
324  break;
325  }
326 
327 
332  default:
333  {
334  return Elem::quality(q);
335  }
336  }
337 
338 
339  // Will never get here...
340  libmesh_error();
341  return 0.;
342 }
343 
344 
345 
346 std::pair<Real, Real> InfHex::qual_bounds (const ElemQuality) const
347 {
348  libmesh_not_implemented();
349 
350  std::pair<Real, Real> bounds;
351 
352 /*
353  switch (q)
354  {
355 
356  case ASPECT_RATIO:
357  bounds.first = 1.;
358  bounds.second = 4.;
359  break;
360 
361  case SKEW:
362  bounds.first = 0.;
363  bounds.second = 0.5;
364  break;
365 
366  case SHEAR:
367  case SHAPE:
368  bounds.first = 0.3;
369  bounds.second = 1.;
370  break;
371 
372  case CONDITION:
373  bounds.first = 1.;
374  bounds.second = 8.;
375  break;
376 
377  case JACOBIAN:
378  bounds.first = 0.5;
379  bounds.second = 1.;
380  break;
381 
382  case DISTORTION:
383  bounds.first = 0.6;
384  bounds.second = 1.;
385  break;
386 
387  case TAPER:
388  bounds.first = 0.;
389  bounds.second = 0.4;
390  break;
391 
392  case STRETCH:
393  bounds.first = 0.25;
394  bounds.second = 1.;
395  break;
396 
397  case DIAGONAL:
398  bounds.first = 0.65;
399  bounds.second = 1.;
400  break;
401 
402  case SIZE:
403  bounds.first = 0.5;
404  bounds.second = 1.;
405  break;
406 
407  default:
408  libMesh::out << "Warning: Invalid quality measure chosen." << std::endl;
409  bounds.first = -1;
410  bounds.second = -1;
411  }
412 */
413 
414  return bounds;
415 }
416 
417 
418 
419 
420 
421 const unsigned short int InfHex::_second_order_adjacent_vertices[8][2] =
422 {
423  { 0, 1}, // vertices adjacent to node 8
424  { 1, 2}, // vertices adjacent to node 9
425  { 2, 3}, // vertices adjacent to node 10
426  { 0, 3}, // vertices adjacent to node 11
427 
428  { 4, 5}, // vertices adjacent to node 12
429  { 5, 6}, // vertices adjacent to node 13
430  { 6, 7}, // vertices adjacent to node 14
431  { 4, 7} // vertices adjacent to node 15
432 };
433 
434 
435 
436 const unsigned short int InfHex::_second_order_vertex_child_number[18] =
437 {
438  99,99,99,99,99,99,99,99, // Vertices
439  0,1,2,0, // Edges
440  0,1,2,0,0, // Faces
441  0 // Interior
442 };
443 
444 
445 
446 const unsigned short int InfHex::_second_order_vertex_child_index[18] =
447 {
448  99,99,99,99,99,99,99,99, // Vertices
449  1,2,3,3, // Edges
450  5,6,7,7,2, // Faces
451  6 // Interior
452 };
453 
454 } // namespace libMesh
455 
456 
457 
458 
459 #endif // ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS

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

Hosted By:
SourceForge.net Logo