face_quad4.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_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 
128  bool proxy)
129 {
130  libmesh_assert_less (i, this->n_sides());
131 
132  if (proxy)
133  return UniquePtr<Elem>(new Side<Edge2,Quad4>(this,i));
134 
135  else
136  {
137  Elem * edge = new Edge2;
138  edge->subdomain_id() = this->subdomain_id();
139 
140  // Set the nodes
141  for (unsigned n=0; n<edge->n_nodes(); ++n)
142  edge->set_node(n) = this->node_ptr(Quad4::side_nodes_map[i][n]);
143 
144  return UniquePtr<Elem>(edge);
145  }
146 
147  libmesh_error_msg("We'll never get here!");
148  return UniquePtr<Elem>();
149 }
150 
151 
152 
153 
154 
155 void Quad4::connectivity(const unsigned int libmesh_dbg_var(sf),
156  const IOPackage iop,
157  std::vector<dof_id_type> & conn) const
158 {
159  libmesh_assert_less (sf, this->n_sub_elem());
160  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
161 
162  // Create storage.
163  conn.resize(4);
164 
165  switch (iop)
166  {
167  case TECPLOT:
168  {
169  conn[0] = this->node_id(0)+1;
170  conn[1] = this->node_id(1)+1;
171  conn[2] = this->node_id(2)+1;
172  conn[3] = this->node_id(3)+1;
173  return;
174  }
175 
176  case VTK:
177  {
178  conn[0] = this->node_id(0);
179  conn[1] = this->node_id(1);
180  conn[2] = this->node_id(2);
181  conn[3] = this->node_id(3);
182  return;
183  }
184 
185  default:
186  libmesh_error_msg("Unsupported IO package " << iop);
187  }
188 }
189 
190 
191 
193 {
194  // Make copies of our points. It makes the subsequent calculations a bit
195  // shorter and avoids dereferencing the same pointer multiple times.
196  Point
197  x0 = point(0), x1 = point(1),
198  x2 = point(2), x3 = point(3);
199 
200  // Construct constant data vectors.
201  // \vec{x}_{\xi} = \vec{a1}*eta + \vec{b1}
202  // \vec{x}_{\eta} = \vec{a2}*xi + \vec{b2}
203  // This is copy-pasted directly from the output of a Python script.
204  Point
205  a1 = x0/4 - x1/4 + x2/4 - x3/4,
206  b1 = -x0/4 + x1/4 + x2/4 - x3/4,
207  a2 = a1,
208  b2 = -x0/4 - x1/4 + x2/4 + x3/4;
209 
210  // Check for quick return for parallelogram QUAD4.
211  if (a1.relative_fuzzy_equals(Point(0,0,0)))
212  return 4. * b1.cross(b2).norm();
213 
214  // Otherwise, use 2x2 quadrature to approximate the surface area.
215 
216  // 4-point rule, exact for bi-cubics. The weights for this rule are
217  // all equal to 1.
218  const Real q[2] = {-std::sqrt(3.)/3, std::sqrt(3.)/3.};
219 
220  Real vol=0.;
221  for (unsigned int i=0; i<2; ++i)
222  for (unsigned int j=0; j<2; ++j)
223  vol += cross_norm(q[j]*a1 + b1,
224  q[i]*a2 + b2);
225 
226  return vol;
227 }
228 
229 } // namespace libMesh
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:1723
virtual bool is_vertex(const unsigned int i) const libmesh_override
Definition: face_quad4.C:92
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const libmesh_override
Definition: face_quad4.C:155
virtual bool has_affine_map() const libmesh_override
Definition: face_quad4.C:119
T cross_norm(const TypeVector< T > &b, const TypeVector< T > &c)
Definition: type_vector.h:1068
static const unsigned int side_nodes_map[4][2]
Definition: face_quad4.h:122
The base class for all geometric element types.
Definition: elem.h:86
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy) libmesh_override
Definition: face_quad4.C:127
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &) const
Definition: type_vector.h:859
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1658
virtual bool is_face(const unsigned int i) const libmesh_override
Definition: face_quad4.C:102
subdomain_id_type subdomain_id() const
Definition: elem.h:1733
PetscErrorCode Vec Mat libmesh_dbg_var(j)
static const float _embedding_matrix[4][4][4]
Definition: face_quad4.h:154
virtual bool is_node_on_side(const unsigned int n, const unsigned int s) const libmesh_override
Definition: face_quad4.C:107
Proxy class for efficiently representing an Elem&#39;s side.
Definition: side.h:48
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1595
virtual bool is_edge(const unsigned int i) const libmesh_override
Definition: face_quad4.C:97
virtual unsigned int n_sides() const libmesh_override
Definition: face_quad.h:85
A 1D geometric element with 2 nodes.
Definition: edge_edge2.h:43
virtual unsigned int n_sub_elem() const libmesh_override
Definition: face_quad4.h:68
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1617
virtual Real volume() const libmesh_override
Definition: face_quad4.C:192
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