face_quad4.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 std::unique_ptr<Elem> Quad4::build_side_ptr (const unsigned int i,
128  bool proxy)
129 {
130  libmesh_assert_less (i, this->n_sides());
131 
132  if (proxy)
133  return libmesh_make_unique<Side<Edge2,Quad4>>(this,i);
134 
135  else
136  {
137  std::unique_ptr<Elem> edge = libmesh_make_unique<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 edge;
145  }
146 }
147 
148 
149 
150 
151 
152 void Quad4::connectivity(const unsigned int libmesh_dbg_var(sf),
153  const IOPackage iop,
154  std::vector<dof_id_type> & conn) const
155 {
156  libmesh_assert_less (sf, this->n_sub_elem());
157  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
158 
159  // Create storage.
160  conn.resize(4);
161 
162  switch (iop)
163  {
164  case TECPLOT:
165  {
166  conn[0] = this->node_id(0)+1;
167  conn[1] = this->node_id(1)+1;
168  conn[2] = this->node_id(2)+1;
169  conn[3] = this->node_id(3)+1;
170  return;
171  }
172 
173  case VTK:
174  {
175  conn[0] = this->node_id(0);
176  conn[1] = this->node_id(1);
177  conn[2] = this->node_id(2);
178  conn[3] = this->node_id(3);
179  return;
180  }
181 
182  default:
183  libmesh_error_msg("Unsupported IO package " << iop);
184  }
185 }
186 
187 
188 
190 {
191  // Make copies of our points. It makes the subsequent calculations a bit
192  // shorter and avoids dereferencing the same pointer multiple times.
193  Point
194  x0 = point(0), x1 = point(1),
195  x2 = point(2), x3 = point(3);
196 
197  // Construct constant data vectors.
198  // \vec{x}_{\xi} = \vec{a1}*eta + \vec{b1}
199  // \vec{x}_{\eta} = \vec{a2}*xi + \vec{b2}
200  // This is copy-pasted directly from the output of a Python script.
201  Point
202  a1 = x0/4 - x1/4 + x2/4 - x3/4,
203  b1 = -x0/4 + x1/4 + x2/4 - x3/4,
204  a2 = a1,
205  b2 = -x0/4 - x1/4 + x2/4 + x3/4;
206 
207  // Check for quick return for parallelogram QUAD4.
208  if (a1.relative_fuzzy_equals(Point(0,0,0)))
209  return 4. * b1.cross(b2).norm();
210 
211  // Otherwise, use 2x2 quadrature to approximate the surface area.
212 
213  // 4-point rule, exact for bi-cubics. The weights for this rule are
214  // all equal to 1.
215  const Real q[2] = {-std::sqrt(3.)/3, std::sqrt(3.)/3.};
216 
217  Real vol=0.;
218  for (unsigned int i=0; i<2; ++i)
219  for (unsigned int j=0; j<2; ++j)
220  vol += cross_norm(q[j]*a1 + b1,
221  q[i]*a2 + b2);
222 
223  return vol;
224 }
225 
226 } // namespace libMesh
virtual bool is_vertex(const unsigned int i) const libmesh_override
Definition: face_quad4.C:92
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy) libmesh_override
Definition: face_quad4.C:127
virtual void connectivity(const unsigned int sf, const IOPackage iop, std::vector< dof_id_type > &conn) const libmesh_override
Definition: face_quad4.C:152
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:1092
static const unsigned int side_nodes_map[4][2]
Definition: face_quad4.h:122
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1875
virtual bool is_face(const unsigned int i) const libmesh_override
Definition: face_quad4.C:102
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition: type_vector.h:879
subdomain_id_type subdomain_id() const
Definition: elem.h:1952
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
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1810
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
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:1832
virtual Real volume() const libmesh_override
Definition: face_quad4.C:189
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:987