face_quad9.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 // C++ includes
19 
20 // Local includes
21 #include "libmesh/side.h"
22 #include "libmesh/edge_edge3.h"
23 #include "libmesh/face_quad9.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 
31 // ------------------------------------------------------------
32 // Quad9 class static member initializations
33 const unsigned int Quad9::side_nodes_map[4][3] =
34 {
35  {0, 1, 4}, // Side 0
36  {1, 2, 5}, // Side 1
37  {2, 3, 6}, // Side 2
38  {3, 0, 7} // Side 3
39 };
40 
41 
42 #ifdef LIBMESH_ENABLE_AMR
43 
44 const float Quad9::_embedding_matrix[4][9][9] =
45 {
46  // embedding matrix for child 0
47  {
48  // 0 1 2 3 4 5 6 7 8
49  { 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 0
50  { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1
51  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 2
52  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 3
53  { 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 4
54  { 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.00000, 0.750000 }, // 5
55  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.750000 }, // 6
56  { 0.375000, 0.00000, 0.00000, -0.125000, 0.00000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 7
57  { 0.140625, -0.0468750, 0.0156250, -0.0468750, 0.281250, -0.0937500, -0.0937500, 0.281250, 0.562500 } // 8
58  },
59 
60  // embedding matrix for child 1
61  {
62  // 0 1 2 3 4 5 6 7 8
63  { 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 0
64  { 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 1
65  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 2
66  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 3
67  { -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 4
68  { 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 5
69  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.750000 }, // 6
70  { 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.00000, 0.750000 }, // 7
71  { -0.0468750, 0.140625, -0.0468750, 0.0156250, 0.281250, 0.281250, -0.0937500, -0.0937500, 0.562500 } // 8
72  },
73 
74  // embedding matrix for child 2
75  {
76  // 0 1 2 3 4 5 6 7 8
77  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000 }, // 0
78  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 1
79  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 2
80  { 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 3
81  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.750000 }, // 4
82  { 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.00000, 0.750000 }, // 5
83  { 0.00000, 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 6
84  { -0.125000, 0.00000, 0.00000, 0.375000, 0.00000, 0.00000, 0.00000, 0.750000, 0.00000 }, // 7
85  { -0.0468750, 0.0156250, -0.0468750, 0.140625, -0.0937500, -0.0937500, 0.281250, 0.281250, 0.562500 } // 8
86  },
87 
88  // embedding matrix for child 3
89  {
90  // 0 1 2 3 4 5 6 7 8
91  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000 }, // 0
92  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000 }, // 1
93  { 0.00000, 0.00000, 1.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000 }, // 2
94  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 1.00000, 0.00000, 0.00000 }, // 3
95  { 0.00000, 0.00000, 0.00000, 0.00000, 0.00000, 0.375000, 0.00000, -0.125000, 0.750000 }, // 4
96  { 0.00000, -0.125000, 0.375000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000, 0.00000 }, // 5
97  { 0.00000, 0.00000, 0.375000, -0.125000, 0.00000, 0.00000, 0.750000, 0.00000, 0.00000 }, // 6
98  { 0.00000, 0.00000, 0.00000, 0.00000, -0.125000, 0.00000, 0.375000, 0.00000, 0.750000 }, // 7
99  { 0.0156250, -0.0468750, 0.140625, -0.0468750, -0.0937500, 0.281250, 0.281250, -0.0937500, 0.562500 } // 8
100  }
101 };
102 
103 #endif
104 
105 
106 
107 // ------------------------------------------------------------
108 // Quad9 class member functions
109 
110 bool Quad9::is_vertex(const unsigned int i) const
111 {
112  if (i < 4)
113  return true;
114  return false;
115 }
116 
117 bool Quad9::is_edge(const unsigned int i) const
118 {
119  if (i < 4)
120  return false;
121  if (i > 7)
122  return false;
123  return true;
124 }
125 
126 bool Quad9::is_face(const unsigned int i) const
127 {
128  if (i > 7)
129  return true;
130  return false;
131 }
132 
133 bool Quad9::is_node_on_side(const unsigned int n,
134  const unsigned int s) const
135 {
136  libmesh_assert_less (s, n_sides());
137  for (unsigned int i = 0; i != 3; ++i)
138  if (side_nodes_map[s][i] == n)
139  return true;
140  return false;
141 }
142 
143 
144 
146 {
147  // make sure corners form a parallelogram
148  Point v = this->point(1) - this->point(0);
149  if (!v.relative_fuzzy_equals(this->point(2) - this->point(3)))
150  return false;
151  // make sure "horizontal" sides are straight
152  v /= 2;
153  if (!v.relative_fuzzy_equals(this->point(4) - this->point(0)) ||
154  !v.relative_fuzzy_equals(this->point(6) - this->point(3)))
155  return false;
156  // make sure "vertical" sides are straight
157  // and the center node is centered
158  v = (this->point(3) - this->point(0))/2;
159  if (!v.relative_fuzzy_equals(this->point(7) - this->point(0)) ||
160  !v.relative_fuzzy_equals(this->point(5) - this->point(1)) ||
161  !v.relative_fuzzy_equals(this->point(8) - this->point(4)))
162  return false;
163  return true;
164 }
165 
166 
167 
168 dof_id_type Quad9::key (const unsigned int s) const
169 {
170  libmesh_assert_less (s, this->n_sides());
171 
172  switch (s)
173  {
174  case 0:
175 
176  return
177  this->compute_key (this->node(4));
178 
179  case 1:
180 
181  return
182  this->compute_key (this->node(5));
183 
184  case 2:
185 
186  return
187  this->compute_key (this->node(6));
188 
189  case 3:
190 
191  return
192  this->compute_key (this->node(7));
193  }
194 
195 
196  // We will never get here... Look at the code above.
197  libmesh_error();
198  return 0;
199 }
200 
201 
202 
203 AutoPtr<Elem> Quad9::build_side (const unsigned int i,
204  bool proxy) const
205 {
206  libmesh_assert_less (i, this->n_sides());
207 
208  if (proxy)
209  {
210  AutoPtr<Elem> ap(new Side<Edge3,Quad9>(this,i));
211  return ap;
212  }
213 
214  else
215  {
216  AutoPtr<Elem> edge(new Edge3);
217  edge->subdomain_id() = this->subdomain_id();
218 
219  switch (i)
220  {
221  case 0:
222  {
223  edge->set_node(0) = this->get_node(0);
224  edge->set_node(1) = this->get_node(1);
225  edge->set_node(2) = this->get_node(4);
226 
227  return edge;
228  }
229  case 1:
230  {
231  edge->set_node(0) = this->get_node(1);
232  edge->set_node(1) = this->get_node(2);
233  edge->set_node(2) = this->get_node(5);
234 
235  return edge;
236  }
237  case 2:
238  {
239  edge->set_node(0) = this->get_node(2);
240  edge->set_node(1) = this->get_node(3);
241  edge->set_node(2) = this->get_node(6);
242 
243  return edge;
244  }
245  case 3:
246  {
247  edge->set_node(0) = this->get_node(3);
248  edge->set_node(1) = this->get_node(0);
249  edge->set_node(2) = this->get_node(7);
250 
251  return edge;
252  }
253  default:
254  {
255  libmesh_error();
256  }
257  }
258  }
259 
260  // We will never get here...
261  AutoPtr<Elem> ap(NULL); return ap;
262 }
263 
264 
265 
266 
267 
268 
269 
270 void Quad9::connectivity(const unsigned int sf,
271  const IOPackage iop,
272  std::vector<dof_id_type>& conn) const
273 {
274  libmesh_assert_less (sf, this->n_sub_elem());
275  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
276 
277  conn.resize(4);
278 
279  switch (iop)
280  {
281  case TECPLOT:
282  {
283  switch(sf)
284  {
285  case 0:
286  // linear sub-quad 0
287  conn[0] = this->node(0)+1;
288  conn[1] = this->node(4)+1;
289  conn[2] = this->node(8)+1;
290  conn[3] = this->node(7)+1;
291  return;
292 
293  case 1:
294  // linear sub-quad 1
295  conn[0] = this->node(4)+1;
296  conn[1] = this->node(1)+1;
297  conn[2] = this->node(5)+1;
298  conn[3] = this->node(8)+1;
299  return;
300 
301  case 2:
302  // linear sub-quad 2
303  conn[0] = this->node(7)+1;
304  conn[1] = this->node(8)+1;
305  conn[2] = this->node(6)+1;
306  conn[3] = this->node(3)+1;
307  return;
308 
309  case 3:
310  // linear sub-quad 3
311  conn[0] = this->node(8)+1;
312  conn[1] = this->node(5)+1;
313  conn[2] = this->node(2)+1;
314  conn[3] = this->node(6)+1;
315  return;
316 
317  default:
318  libmesh_error();
319  }
320  }
321 
322  case VTK:
323  {
324  conn.resize(9);
325  conn[0] = this->node(0);
326  conn[1] = this->node(1);
327  conn[2] = this->node(2);
328  conn[3] = this->node(3);
329  conn[4] = this->node(4);
330  conn[5] = this->node(5);
331  conn[6] = this->node(6);
332  conn[7] = this->node(7);
333  conn[8] = this->node(8);
334  return;
335 
336  /*
337  switch(sf)
338  {
339  case 0:
340  // linear sub-quad 0
341  conn[0] = this->node(0);
342  conn[1] = this->node(4);
343  conn[2] = this->node(8);
344  conn[3] = this->node(7);
345 
346  return;
347 
348  case 1:
349  // linear sub-quad 1
350  conn[0] = this->node(4);
351  conn[1] = this->node(1);
352  conn[2] = this->node(5);
353  conn[3] = this->node(8);
354 
355  return;
356 
357  case 2:
358  // linear sub-quad 2
359  conn[0] = this->node(7);
360  conn[1] = this->node(8);
361  conn[2] = this->node(6);
362  conn[3] = this->node(3);
363 
364  return;
365 
366  case 3:
367  // linear sub-quad 3
368  conn[0] = this->node(8);
369  conn[1] = this->node(5);
370  conn[2] = this->node(2);
371  conn[3] = this->node(6);
372 
373  return;
374 
375  default:
376  libmesh_error();
377  }*/
378  }
379 
380  default:
381  {
382  libmesh_error();
383  }
384  }
385 
386  libmesh_error();
387 }
388 
389 
390 
391 
392 unsigned int Quad9::n_second_order_adjacent_vertices (const unsigned int n) const
393 {
394  switch (n)
395  {
396  case 4:
397  case 5:
398  case 6:
399  case 7:
400  return 2;
401 
402  case 8:
403  return 4;
404 
405  default:
406  libmesh_error();
407  }
408 
409  libmesh_error();
410  return libMesh::invalid_uint;
411 }
412 
413 
414 
415 unsigned short int Quad9::second_order_adjacent_vertex (const unsigned int n,
416  const unsigned int v) const
417 {
418  libmesh_assert_greater_equal (n, this->n_vertices());
419  libmesh_assert_less (n, this->n_nodes());
420 
421  switch (n)
422  {
423  case 8:
424  {
425  libmesh_assert_less (v, 4);
426  return static_cast<unsigned short int>(v);
427  }
428 
429  default:
430  {
431  libmesh_assert_less (v, 2);
432  // use the matrix that we inherited from \p Quad
433  return _second_order_adjacent_vertices[n-this->n_vertices()][v];
434  }
435  }
436 }
437 
438 
439 
440 std::pair<unsigned short int, unsigned short int>
441 Quad9::second_order_child_vertex (const unsigned int n) const
442 {
443  libmesh_assert_greater_equal (n, this->n_vertices());
444  libmesh_assert_less (n, this->n_nodes());
445  /*
446  * the _second_order_vertex_child_* vectors are
447  * stored in face_quad.C, since they are identical
448  * for Quad8 and Quad9 (for the first 4 higher-order nodes)
449  */
450  return std::pair<unsigned short int, unsigned short int>
453 }
454 
455 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo