libMesh
edge_edge4.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 
19 
20 // Local includes
21 #include "libmesh/edge_edge4.h"
22 
23 namespace libMesh
24 {
25 
26 #ifdef LIBMESH_ENABLE_AMR
27 
28 const float Edge4::_embedding_matrix[2][4][4] =
29  {
30  // embedding matrix for child 0
31 
32  {
33  // 0 1 2 3 // Shape function index
34  {1.0, 0.0, 0.0, 0.0}, // left, xi = -1
35  {-0.0625, -0.0625, 0.5625, 0.5625}, // right, xi = 0
36  {0.3125, 0.0625, 0.9375, -0.3125}, // middle left, xi = -2/3
37  {0.0, 0.0, 1.0, 0.0} // middle right, xi = -1/3
38  },
39 
40  // embedding matrix for child 1
41  {
42  // 0 1 2 3 // Shape function index
43  {-0.0625, -0.0625, 0.5625, 0.5625}, // left, xi = 0
44  {0.0, 1.0, 0.0, 0.0}, // right, xi = 1
45  {0.0, 0.0, 0.0, 1.0}, // middle left, xi = 1/3
46  {0.0625, 0.3125, -0.3125, 0.9375} // middle right, xi = 2/3
47  }
48  };
49 
50 #endif
51 
52 bool Edge4::is_vertex(const unsigned int i) const
53 {
54  return (i==0) || (i==1);
55 }
56 
57 bool Edge4::is_edge(const unsigned int i) const
58 {
59  return (i==2) || (i==3);
60 }
61 
62 bool Edge4::is_face(const unsigned int ) const
63 {
64  return false;
65 }
66 
67 bool Edge4::is_node_on_side(const unsigned int n,
68  const unsigned int s) const
69 {
70  libmesh_assert_less (s, 2);
71  libmesh_assert_less (n, 4);
72  return (s == n);
73 }
74 
75 bool Edge4::is_node_on_edge(const unsigned int,
76  const unsigned int libmesh_dbg_var(e)) const
77 {
78  libmesh_assert_equal_to (e, 0);
79  return true;
80 }
81 
82 
83 
84 bool Edge4::has_affine_map() const
85 {
86  if (!this->point(2).relative_fuzzy_equals
87  ((this->point(0)*2. + this->point(1))/3.))
88  return false;
89  if (!this->point(3).relative_fuzzy_equals
90  ((this->point(0) + this->point(1)*2.)/3.))
91  return false;
92  return true;
93 }
94 
95 
96 
97 void Edge4::connectivity(const unsigned int sc,
98  const IOPackage iop,
99  std::vector<dof_id_type> & conn) const
100 {
101  libmesh_assert_less_equal (sc, 2);
102  libmesh_assert_less (sc, this->n_sub_elem());
103  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
104 
105  // Create storage
106  conn.resize(2);
107 
108  switch(iop)
109  {
110  case TECPLOT:
111  {
112  switch (sc)
113  {
114  case 0:
115  conn[0] = this->node_id(0)+1;
116  conn[1] = this->node_id(2)+1;
117  return;
118 
119  case 1:
120  conn[0] = this->node_id(2)+1;
121  conn[1] = this->node_id(3)+1;
122  return;
123 
124  case 2:
125  conn[0] = this->node_id(3)+1;
126  conn[1] = this->node_id(1)+1;
127  return;
128 
129  default:
130  libmesh_error_msg("Invalid sc = " << sc);
131  }
132 
133  }
134 
135  case VTK:
136  {
137 
138  switch (sc)
139  {
140  case 0:
141  conn[0] = this->node_id(0);
142  conn[1] = this->node_id(2);
143  return;
144 
145  case 1:
146  conn[0] = this->node_id(2);
147  conn[1] = this->node_id(3);
148  return;
149 
150  case 2:
151  conn[0] = this->node_id(3);
152  conn[1] = this->node_id(1);
153  return;
154 
155  default:
156  libmesh_error_msg("Invalid sc = " << sc);
157  }
158  }
159 
160  default:
161  libmesh_error_msg("Unsupported IO package " << iop);
162  }
163 }
164 
165 
166 
167 BoundingBox Edge4::loose_bounding_box () const
168 {
169  // This might be a curved line through 2-space or 3-space, in which
170  // case the full bounding box can be larger than the bounding box of
171  // just the nodes.
172  //
173  // FIXME - I haven't yet proven the formula below to be correct for
174  // cubics - RHS
175  Point pmin, pmax;
176 
177  for (unsigned d=0; d<LIBMESH_DIM; ++d)
178  {
179  Real center = (this->point(2)(d) + this->point(3)(d))/2;
180  Real hd = std::max(std::abs(center - this->point(0)(d)),
181  std::abs(center - this->point(1)(d)));
182 
183  pmin(d) = center - hd;
184  pmax(d) = center + hd;
185  }
186 
187  return BoundingBox(pmin, pmax);
188 }
189 
190 
191 
192 dof_id_type Edge4::key () const
193 {
194  return this->compute_key(this->node_id(0),
195  this->node_id(1),
196  this->node_id(2),
197  this->node_id(3));
198 }
199 
200 
201 
202 Real Edge4::volume () const
203 {
204  // Make copies of our points. It makes the subsequent calculations a bit
205  // shorter and avoids dereferencing the same pointer multiple times.
206  Point
207  x0 = point(0),
208  x1 = point(1),
209  x2 = point(2),
210  x3 = point(3);
211 
212  // Construct constant data vectors.
213  // \vec{x}_{\xi} = \vec{a1}*xi**2 + \vec{b1}*xi + \vec{c1}
214  // This is copy-pasted directly from the output of a Python script.
215  Point
216  a1 = -27*x0/16 + 27*x1/16 + 81*x2/16 - 81*x3/16,
217  b1 = 9*x0/8 + 9*x1/8 - 9*x2/8 - 9*x3/8,
218  c1 = x0/16 - x1/16 - 27*x2/16 + 27*x3/16;
219 
220  // 4 point quadrature, 7th-order accurate
221  const unsigned int N = 4;
222  const Real q[N] = {-std::sqrt(525 + 70*std::sqrt(30.)) / 35,
223  -std::sqrt(525 - 70*std::sqrt(30.)) / 35,
224  std::sqrt(525 - 70*std::sqrt(30.)) / 35,
225  std::sqrt(525 + 70*std::sqrt(30.)) / 35};
226  const Real w[N] = {(18 - std::sqrt(30.)) / 36,
227  (18 + std::sqrt(30.)) / 36,
228  (18 + std::sqrt(30.)) / 36,
229  (18 - std::sqrt(30.)) / 36};
230 
231  Real vol=0.;
232  for (unsigned int i=0; i<N; ++i)
233  vol += w[i] * (q[i]*q[i]*a1 + q[i]*b1 + c1).norm();
234 
235  return vol;
236 }
237 
238 } // namespace libMesh
double abs(double a)
The libMesh namespace provides an interface to certain functionality in the library.
long double max(long double a, double b)
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IOPackage
libMesh interfaces with several different software packages for the purposes of creating, reading, and writing mesh files.
uint8_t dof_id_type
Definition: id_types.h:64