face_quad4.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_edge2.h"
23 #include "libmesh/face_quad4.h"
24 
25 namespace libMesh
26 {
27 
28 
29 
30 
31 // ------------------------------------------------------------
32 // Quad class static member initialization
33 const unsigned int Quad4::side_nodes_map[4][2] =
34 {
35  {0, 1}, // Side 0
36  {1, 2}, // Side 1
37  {2, 3}, // Side 2
38  {3, 0} // Side 3
39 };
40 
41 
42 #ifdef LIBMESH_ENABLE_AMR
43 
44 const float Quad4::_embedding_matrix[4][4][4] =
45 {
46  // embedding matrix for child 0
47  {
48  // 0 1 2 3
49  {1.0, 0.0, 0.0, 0.0}, // 0
50  {0.5, 0.5, 0.0, 0.0}, // 1
51  {.25, .25, .25, .25}, // 2
52  {0.5, 0.0, 0.0, 0.5} // 3
53  },
54 
55  // embedding matrix for child 1
56  {
57  // 0 1 2 3
58  {0.5, 0.5, 0.0, 0.0}, // 0
59  {0.0, 1.0, 0.0, 0.0}, // 1
60  {0.0, 0.5, 0.5, 0.0}, // 2
61  {.25, .25, .25, .25} // 3
62  },
63 
64  // embedding matrix for child 2
65  {
66  // 0 1 2 3
67  {0.5, 0.0, 0.0, 0.5}, // 0
68  {.25, .25, .25, .25}, // 1
69  {0.0, 0.0, 0.5, 0.5}, // 2
70  {0.0, 0.0, 0.0, 1.0} // 3
71  },
72 
73  // embedding matrix for child 3
74  {
75  // 0 1 2 3
76  {.25, .25, .25, .25}, // 0
77  {0.0, 0.5, 0.5, 0.0}, // 1
78  {0.0, 0.0, 1.0, 0.0}, // 2
79  {0.0, 0.0, 0.5, 0.5} // 3
80  }
81 };
82 
83 #endif
84 
85 
86 
87 
88 
89 // ------------------------------------------------------------
90 // Quad4 class member functions
91 
92 bool Quad4::is_vertex(const unsigned int) const
93 {
94  return true;
95 }
96 
97 bool Quad4::is_edge(const unsigned int) const
98 {
99  return false;
100 }
101 
102 bool Quad4::is_face(const unsigned int) const
103 {
104  return false;
105 }
106 
107 bool Quad4::is_node_on_side(const unsigned int n,
108  const unsigned int s) const
109 {
110  libmesh_assert_less (s, n_sides());
111  for (unsigned int i = 0; i != 2; ++i)
112  if (side_nodes_map[s][i] == n)
113  return true;
114  return false;
115 }
116 
117 
118 
120 {
121  Point v = this->point(3) - this->point(0);
122  return (v.relative_fuzzy_equals(this->point(2) - this->point(1)));
123 }
124 
125 
126 
127 AutoPtr<Elem> Quad4::build_side (const unsigned int i,
128  bool proxy) const
129 {
130  libmesh_assert_less (i, this->n_sides());
131 
132  if (proxy)
133  {
134  AutoPtr<Elem> ap(new Side<Edge2,Quad4>(this,i));
135  return ap;
136  }
137 
138  else
139  {
140  AutoPtr<Elem> edge(new Edge2);
141  edge->subdomain_id() = this->subdomain_id();
142 
143  switch (i)
144  {
145  case 0:
146  {
147  edge->set_node(0) = this->get_node(0);
148  edge->set_node(1) = this->get_node(1);
149 
150  return edge;
151  }
152  case 1:
153  {
154  edge->set_node(0) = this->get_node(1);
155  edge->set_node(1) = this->get_node(2);
156 
157  return edge;
158  }
159  case 2:
160  {
161  edge->set_node(0) = this->get_node(2);
162  edge->set_node(1) = this->get_node(3);
163 
164  return edge;
165  }
166  case 3:
167  {
168  edge->set_node(0) = this->get_node(3);
169  edge->set_node(1) = this->get_node(0);
170 
171  return edge;
172  }
173  default:
174  {
175  libmesh_error();
176  }
177  }
178  }
179 
180  // We will never get here...
181  AutoPtr<Elem> ap(NULL); return ap;
182 }
183 
184 
185 
186 
187 
188 void Quad4::connectivity(const unsigned int libmesh_dbg_var(sf),
189  const IOPackage iop,
190  std::vector<dof_id_type>& conn) const
191 {
192  libmesh_assert_less (sf, this->n_sub_elem());
193  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
194 
195  // Create storage.
196  conn.resize(4);
197 
198  switch (iop)
199  {
200  case TECPLOT:
201  {
202  conn[0] = this->node(0)+1;
203  conn[1] = this->node(1)+1;
204  conn[2] = this->node(2)+1;
205  conn[3] = this->node(3)+1;
206  return;
207  }
208 
209  case VTK:
210  {
211  conn[0] = this->node(0);
212  conn[1] = this->node(1);
213  conn[2] = this->node(2);
214  conn[3] = this->node(3);
215  return;
216  }
217 
218  default:
219  libmesh_error();
220  }
221 
222  libmesh_error();
223 }
224 
225 
226 
228 {
229  // The A,B,C,D naming scheme here corresponds exactly to the
230  // libmesh counter-clockwise numbering scheme.
231 
232  // 3 2 D C
233  // QUAD4: o-----------o o-----------o
234  // | | | |
235  // | | | |
236  // | | | |
237  // | | | |
238  // | | | |
239  // o-----------o o-----------o
240  // 0 1 A B
241 
242  // Vector pointing from A to C
243  Point AC ( this->point(2) - this->point(0) );
244 
245  // Vector pointing from A to B
246  Point AB ( this->point(1) - this->point(0) );
247 
248  // Vector pointing from A to D
249  Point AD ( this->point(3) - this->point(0) );
250 
251  // The diagonal vector minus the side vectors
252  Point AC_AB_AD (AC - AB - AD);
253 
254  // Check for quick return for planar QUAD4. This will
255  // be the most common case, occuring for all purely 2D meshes.
256  if (AC_AB_AD == Point(0.,0.,0.))
257  return AB.cross(AD).size();
258 
259  else
260  {
261  // Use 2x2 quadrature to approximate the surface area. (The
262  // true integral is too difficult to compute analytically.) The
263  // accuracy here is exactly the same as would be obtained via a
264  // call to Elem::volume(), however it is a bit more optimized to
265  // do it this way. The technique used is to integrate the magnitude
266  // of the normal vector over the whole area. See for example,
267  //
268  // Y. Zhang, C. Bajaj, G. Xu. Surface Smoothing and Quality
269  // Improvement of Quadrilateral/Hexahedral Meshes with Geometric
270  // Flow. The special issue of the Journal Communications in
271  // Numerical Methods in Engineering (CNME), submitted as an
272  // invited paper, 2006.
273  // http://www.ices.utexas.edu/~jessica/paper/quadhexgf/quadhex_geomflow_CNM.pdf
274 
275  // 4-point rule
276  const Real q[2] = {0.5 - std::sqrt(3.) / 6.,
277  0.5 + std::sqrt(3.) / 6.};
278 
279  Real vol=0.;
280  for (unsigned int i=0; i<2; ++i)
281  for (unsigned int j=0; j<2; ++j)
282  vol += (AB + q[i]*AC_AB_AD).cross(AD + q[j]*AC_AB_AD).size();
283 
284  return 0.25*vol;
285  }
286 }
287 
288 } // namespace libMesh

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

Hosted By:
SourceForge.net Logo