libMesh
fe_boundary.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 // C++ includes
21 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 #include <cmath> // for std::sqrt
23 
24 
25 // Local includes
26 #include "libmesh/libmesh_common.h"
27 #include "libmesh/fe.h"
28 #include "libmesh/quadrature.h"
29 #include "libmesh/elem.h"
30 #include "libmesh/libmesh_logging.h"
31 #include "libmesh/tensor_value.h" // May be necessary if destructors
32 // get instantiated here
33 
34 namespace libMesh
35 {
36 
37 //-------------------------------------------------------
38 // Full specializations for useless methods in 0D, 1D
39 #define REINIT_ERROR(_dim, _type, _func) \
40  template <> \
41  void FE<_dim,_type>::_func(const Elem *, \
42  const unsigned int, \
43  const Real, \
44  const std::vector<Point> * const, \
45  const std::vector<Real> * const)
46 
47 #define SIDEMAP_ERROR(_dim, _type, _func) \
48  template <> \
49  void FE<_dim,_type>::_func(const Elem *, \
50  const Elem *, \
51  const unsigned int, \
52  const std::vector<Point> &, \
53  std::vector<Point> &)
54 
55 #define FACE_EDGE_SHAPE_ERROR(_dim, _func) \
56  template <> \
57  void FEMap::_func<_dim>(const std::vector<Point> &, \
58  const Elem *) \
59  { \
60  libmesh_error_msg("ERROR: This method makes no sense for low-D elements!"); \
61  }
62 
63 
64 // 0D error instantiations
65 REINIT_ERROR(0, CLOUGH, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D CLOUGH elements!"); }
66 REINIT_ERROR(0, CLOUGH, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D CLOUGH elements!"); }
67 SIDEMAP_ERROR(0, CLOUGH, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D CLOUGH elements!"); }
68 
69 REINIT_ERROR(0, HERMITE, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D HERMITE elements!"); }
70 REINIT_ERROR(0, HERMITE, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D HERMITE elements!"); }
71 SIDEMAP_ERROR(0, HERMITE, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D HERMITE elements!"); }
72 
73 REINIT_ERROR(0, HIERARCHIC, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D HIERARCHIC elements!"); }
74 REINIT_ERROR(0, HIERARCHIC, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D HIERARCHIC elements!"); }
75 SIDEMAP_ERROR(0, HIERARCHIC, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D HIERARCHIC elements!"); }
76 
77 REINIT_ERROR(0, L2_HIERARCHIC, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D L2_HIERARCHIC elements!"); }
78 REINIT_ERROR(0, L2_HIERARCHIC, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D L2_HIERARCHIC elements!"); }
79 SIDEMAP_ERROR(0, L2_HIERARCHIC, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D L2_HIERARCHIC elements!"); }
80 
81 REINIT_ERROR(0, LAGRANGE, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D LAGRANGE elements!"); }
82 REINIT_ERROR(0, LAGRANGE, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D LAGRANGE elements!"); }
83 SIDEMAP_ERROR(0, LAGRANGE, side_map) { libmesh_error_msg("ERROR: Cannot side_map LAGRANGE elements!"); }
84 
85 REINIT_ERROR(0, LAGRANGE_VEC, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D LAGRANGE_VEC elements!"); }
86 REINIT_ERROR(0, LAGRANGE_VEC, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D LAGRANGE_VEC elements!"); }
87 SIDEMAP_ERROR(0, LAGRANGE_VEC, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D LAGRANGE_VEC elements!"); }
88 
89 REINIT_ERROR(0, L2_LAGRANGE, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D L2_LAGRANGE elements!"); }
90 REINIT_ERROR(0, L2_LAGRANGE, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D L2_LAGRANGE elements!"); }
91 SIDEMAP_ERROR(0, L2_LAGRANGE, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D L2_LAGRANGE elements!"); }
92 
93 REINIT_ERROR(0, MONOMIAL, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D MONOMIAL elements!"); }
94 REINIT_ERROR(0, MONOMIAL, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D MONOMIAL elements!"); }
95 SIDEMAP_ERROR(0, MONOMIAL, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D MONOMIAL elements!"); }
96 
97 REINIT_ERROR(0, SCALAR, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D SCALAR elements!"); }
98 REINIT_ERROR(0, SCALAR, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D SCALAR elements!"); }
99 SIDEMAP_ERROR(0, SCALAR, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D SCALAR elements!"); }
100 
101 REINIT_ERROR(0, XYZ, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D XYZ elements!"); }
102 REINIT_ERROR(0, XYZ, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D XYZ elements!"); }
103 SIDEMAP_ERROR(0, XYZ, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D XYZ elements!"); }
104 
105 REINIT_ERROR(0, NEDELEC_ONE, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D NEDELEC_ONE elements!"); }
106 REINIT_ERROR(0, NEDELEC_ONE, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D NEDELEC_ONE elements!"); }
107 SIDEMAP_ERROR(0, NEDELEC_ONE, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D NEDELEC_ONE elements!"); }
108 
109 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
110 REINIT_ERROR(0, BERNSTEIN, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D BERNSTEIN elements!"); }
111 REINIT_ERROR(0, BERNSTEIN, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D BERNSTEIN elements!"); }
112 SIDEMAP_ERROR(0, BERNSTEIN, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D BERNSTEIN elements!"); }
113 
114 REINIT_ERROR(0, SZABAB, reinit) { libmesh_error_msg("ERROR: Cannot reinit 0D SZABAB elements!"); }
115 REINIT_ERROR(0, SZABAB, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 0D SZABAB elements!"); }
116 SIDEMAP_ERROR(0, SZABAB, side_map) { libmesh_error_msg("ERROR: Cannot side_map 0D SZABAB elements!"); }
117 #endif
118 
119 // 1D error instantiations
120 REINIT_ERROR(1, CLOUGH, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D CLOUGH elements!"); }
121 REINIT_ERROR(1, HERMITE, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D HERMITE elements!"); }
122 REINIT_ERROR(1, HIERARCHIC, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D HIERARCHIC elements!"); }
123 REINIT_ERROR(1, L2_HIERARCHIC, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D L2_HIERARCHIC elements!"); }
124 REINIT_ERROR(1, LAGRANGE, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D LAGRANGE elements!"); }
125 REINIT_ERROR(1, LAGRANGE_VEC, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D LAGRANGE_VEC elements!"); }
126 REINIT_ERROR(1, L2_LAGRANGE, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D L2_LAGRANGE elements!"); }
127 REINIT_ERROR(1, XYZ, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D XYZ elements!"); }
128 REINIT_ERROR(1, MONOMIAL, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D MONOMIAL elements!"); }
129 REINIT_ERROR(1, SCALAR, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D SCALAR elements!"); }
130 REINIT_ERROR(1, NEDELEC_ONE, reinit) { libmesh_error_msg("ERROR: Cannot reinit 1D NEDELEC_ONE elements!"); }
131 REINIT_ERROR(1, NEDELEC_ONE, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D NEDELEC_ONE elements!"); }
132 SIDEMAP_ERROR(1, NEDELEC_ONE, side_map) { libmesh_error_msg("ERROR: Cannot side_map 1D NEDELEC_ONE elements!"); }
133 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
134 REINIT_ERROR(1, BERNSTEIN, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D BERNSTEIN elements!"); }
135 REINIT_ERROR(1, SZABAB, edge_reinit) { libmesh_error_msg("ERROR: Cannot edge_reinit 1D SZABAB elements!"); }
136 #endif
137 
138 
139 //-------------------------------------------------------
140 // Methods for 2D, 3D
141 template <unsigned int Dim, FEFamily T>
142 void FE<Dim,T>::reinit(const Elem * elem,
143  const unsigned int s,
144  const Real /* tolerance */,
145  const std::vector<Point> * const pts,
146  const std::vector<Real> * const weights)
147 {
148  libmesh_assert(elem);
149  libmesh_assert (this->qrule != libmesh_nullptr || pts != libmesh_nullptr);
150  // We now do this for 1D elements!
151  // libmesh_assert_not_equal_to (Dim, 1);
152 
153  // We're (possibly re-) calculating now! FIXME - we currently
154  // expect to be able to use side_map and JxW later, but we could
155  // optimize further here.
156  this->_fe_map->calculations_started = false;
157  this->_fe_map->get_JxW();
158  this->_fe_map->get_xyz();
159  this->determine_calculations();
160 
161  // Build the side of interest
162  const UniquePtr<const Elem> side(elem->build_side_ptr(s));
163 
164  // Find the max p_level to select
165  // the right quadrature rule for side integration
166  unsigned int side_p_level = elem->p_level();
167  if (elem->neighbor_ptr(s) != libmesh_nullptr)
168  side_p_level = std::max(side_p_level, elem->neighbor_ptr(s)->p_level());
169 
170  // Initialize the shape functions at the user-specified
171  // points
172  if (pts != libmesh_nullptr)
173  {
174  // The shape functions do not correspond to the qrule
175  this->shapes_on_quadrature = false;
176 
177  // Initialize the face shape functions
178  this->_fe_map->template init_face_shape_functions<Dim>(*pts, side.get());
179 
180  // Compute the Jacobian*Weight on the face for integration
181  if (weights != libmesh_nullptr)
182  {
183  this->_fe_map->compute_face_map (Dim, *weights, side.get());
184  }
185  else
186  {
187  std::vector<Real> dummy_weights (pts->size(), 1.);
188  this->_fe_map->compute_face_map (Dim, dummy_weights, side.get());
189  }
190  }
191  // If there are no user specified points, we use the
192  // quadrature rule
193  else
194  {
195  // initialize quadrature rule
196  this->qrule->init(side->type(), side_p_level);
197 
198  if (this->qrule->shapes_need_reinit())
199  this->shapes_on_quadrature = false;
200 
201  // FIXME - could this break if the same FE object was used
202  // for both volume and face integrals? - RHS
203  // We might not need to reinitialize the shape functions
204  if ((this->get_type() != elem->type()) ||
205  (side->type() != last_side) ||
206  (this->get_p_level() != side_p_level) ||
207  this->shapes_need_reinit() ||
208  !this->shapes_on_quadrature)
209  {
210  // Set the element type and p_level
211  this->elem_type = elem->type();
212 
213  // Set the last_side
214  last_side = side->type();
215 
216  // Set the last p level
217  this->_p_level = side_p_level;
218 
219  // Initialize the face shape functions
220  this->_fe_map->template init_face_shape_functions<Dim>(this->qrule->get_points(), side.get());
221  }
222 
223  // Compute the Jacobian*Weight on the face for integration
224  this->_fe_map->compute_face_map (Dim, this->qrule->get_weights(), side.get());
225 
226  // The shape functions correspond to the qrule
227  this->shapes_on_quadrature = true;
228  }
229 
230  // make a copy of the Jacobian for integration
231  const std::vector<Real> JxW_int(this->_fe_map->get_JxW());
232 
233  // make a copy of shape on quadrature info
234  bool shapes_on_quadrature_side = this->shapes_on_quadrature;
235 
236  // Find where the integration points are located on the
237  // full element.
238  const std::vector<Point> * ref_qp;
239  if (pts != libmesh_nullptr)
240  ref_qp = pts;
241  else
242  ref_qp = &this->qrule->get_points();
243 
244  std::vector<Point> qp;
245  this->side_map(elem, side.get(), s, *ref_qp, qp);
246 
247  // compute the shape function and derivative values
248  // at the points qp
249  this->reinit (elem, &qp);
250 
251  this->shapes_on_quadrature = shapes_on_quadrature_side;
252 
253  // copy back old data
254  this->_fe_map->get_JxW() = JxW_int;
255 }
256 
257 
258 
259 template <unsigned int Dim, FEFamily T>
260 void FE<Dim,T>::edge_reinit(const Elem * elem,
261  const unsigned int e,
262  const Real tolerance,
263  const std::vector<Point> * const pts,
264  const std::vector<Real> * const weights)
265 {
266  libmesh_assert(elem);
267  libmesh_assert (this->qrule != libmesh_nullptr || pts != libmesh_nullptr);
268  // We don't do this for 1D elements!
269  libmesh_assert_not_equal_to (Dim, 1);
270 
271  // We're (possibly re-) calculating now! Time to determine what.
272  // FIXME - we currently just assume that we're using JxW and calling
273  // edge_map later.
274  this->_fe_map->calculations_started = false;
275  this->_fe_map->get_JxW();
276  this->_fe_map->get_xyz();
277  this->determine_calculations();
278 
279  // Build the side of interest
280  const UniquePtr<const Elem> edge(elem->build_edge_ptr(e));
281 
282  // Initialize the shape functions at the user-specified
283  // points
284  if (pts != libmesh_nullptr)
285  {
286  // The shape functions do not correspond to the qrule
287  this->shapes_on_quadrature = false;
288 
289  // Initialize the edge shape functions
290  this->_fe_map->template init_edge_shape_functions<Dim> (*pts, edge.get());
291 
292  // Compute the Jacobian*Weight on the face for integration
293  if (weights != libmesh_nullptr)
294  {
295  this->_fe_map->compute_edge_map (Dim, *weights, edge.get());
296  }
297  else
298  {
299  std::vector<Real> dummy_weights (pts->size(), 1.);
300  this->_fe_map->compute_edge_map (Dim, dummy_weights, edge.get());
301  }
302  }
303  // If there are no user specified points, we use the
304  // quadrature rule
305  else
306  {
307  // initialize quadrature rule
308  this->qrule->init(edge->type(), elem->p_level());
309 
310  if (this->qrule->shapes_need_reinit())
311  this->shapes_on_quadrature = false;
312 
313  // We might not need to reinitialize the shape functions
314  if ((this->get_type() != elem->type()) ||
315  (edge->type() != static_cast<int>(last_edge)) || // Comparison between enum and unsigned, cast the unsigned to int
316  this->shapes_need_reinit() ||
317  !this->shapes_on_quadrature)
318  {
319  // Set the element type
320  this->elem_type = elem->type();
321 
322  // Set the last_edge
323  last_edge = edge->type();
324 
325  // Initialize the edge shape functions
326  this->_fe_map->template init_edge_shape_functions<Dim> (this->qrule->get_points(), edge.get());
327  }
328 
329  // Compute the Jacobian*Weight on the face for integration
330  this->_fe_map->compute_edge_map (Dim, this->qrule->get_weights(), edge.get());
331 
332  // The shape functions correspond to the qrule
333  this->shapes_on_quadrature = true;
334  }
335 
336  // make a copy of the Jacobian for integration
337  const std::vector<Real> JxW_int(this->_fe_map->get_JxW());
338 
339  // Find where the integration points are located on the
340  // full element.
341  std::vector<Point> qp;
342  this->inverse_map (elem, this->_fe_map->get_xyz(), qp, tolerance);
343 
344  // compute the shape function and derivative values
345  // at the points qp
346  this->reinit (elem, &qp);
347 
348  // copy back old data
349  this->_fe_map->get_JxW() = JxW_int;
350 }
351 
352 template <unsigned int Dim, FEFamily T>
353 void FE<Dim,T>::side_map (const Elem * elem,
354  const Elem * side,
355  const unsigned int s,
356  const std::vector<Point> & reference_side_points,
357  std::vector<Point> & reference_points)
358 {
359  // We're calculating mappings - we need at least first order info
360  this->calculate_phi = true;
361  this->determine_calculations();
362 
363  unsigned int side_p_level = elem->p_level();
364  if (elem->neighbor_ptr(s) != libmesh_nullptr)
365  side_p_level = std::max(side_p_level, elem->neighbor_ptr(s)->p_level());
366 
367  if (side->type() != last_side ||
368  side_p_level != this->_p_level ||
369  !this->shapes_on_quadrature)
370  {
371  // Set the element type
372  this->elem_type = elem->type();
373  this->_p_level = side_p_level;
374 
375  // Set the last_side
376  last_side = side->type();
377 
378  // Initialize the face shape functions
379  this->_fe_map->template init_face_shape_functions<Dim>(reference_side_points, side);
380  }
381 
382  const unsigned int n_points =
383  cast_int<unsigned int>(reference_side_points.size());
384  reference_points.resize(n_points);
385  for (unsigned int i = 0; i < n_points; i++)
386  reference_points[i].zero();
387 
388  std::vector<unsigned int> elem_nodes_map;
389  elem_nodes_map.resize(side->n_nodes());
390  for (unsigned int j = 0; j < side->n_nodes(); j++)
391  for (unsigned int i = 0; i < elem->n_nodes(); i++)
392  if (side->node_id(j) == elem->node_id(i))
393  elem_nodes_map[j] = i;
394  std::vector<Point> refspace_nodes;
395  this->get_refspace_nodes(elem->type(), refspace_nodes);
396 
397  const std::vector<std::vector<Real>> & psi_map = this->_fe_map->get_psi();
398 
399  for (std::size_t i=0; i<psi_map.size(); i++) // sum over the nodes
400  {
401  const Point & side_node = refspace_nodes[elem_nodes_map[i]];
402  for (unsigned int p=0; p<n_points; p++)
403  reference_points[p].add_scaled (side_node, psi_map[i][p]);
404  }
405 }
406 
407 template<unsigned int Dim>
408 void FEMap::init_face_shape_functions(const std::vector<Point> & qp,
409  const Elem * side)
410 {
411  // Start logging the shape function initialization
412  LOG_SCOPE("init_face_shape_functions()", "FEMap");
413 
414  libmesh_assert(side);
415 
416  // We're calculating now!
417  this->determine_calculations();
418 
419  // The element type and order to use in
420  // the map
421  const Order mapping_order (side->default_order());
422  const ElemType mapping_elem_type (side->type());
423 
424  // The number of quadrature points.
425  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
426 
427  const unsigned int n_mapping_shape_functions =
428  FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
429  mapping_order);
430 
431  // resize the vectors to hold current data
432  // Psi are the shape functions used for the FE mapping
433  if (calculate_xyz)
434  this->psi_map.resize (n_mapping_shape_functions);
435 
436  if (Dim > 1)
437  {
438  if (calculate_dxyz)
439  this->dpsidxi_map.resize (n_mapping_shape_functions);
440  if (calculate_d2xyz)
441  this->d2psidxi2_map.resize (n_mapping_shape_functions);
442  }
443 
444  if (Dim == 3)
445  {
446  if (calculate_dxyz)
447  this->dpsideta_map.resize (n_mapping_shape_functions);
448  if (calculate_d2xyz)
449  {
450  this->d2psidxideta_map.resize (n_mapping_shape_functions);
451  this->d2psideta2_map.resize (n_mapping_shape_functions);
452  }
453  }
454 
455  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
456  {
457  // Allocate space to store the values of the shape functions
458  // and their first and second derivatives at the quadrature points.
459  if (calculate_xyz)
460  this->psi_map[i].resize (n_qp);
461  if (Dim > 1)
462  {
463  if (calculate_dxyz)
464  this->dpsidxi_map[i].resize (n_qp);
465  if (calculate_d2xyz)
466  this->d2psidxi2_map[i].resize (n_qp);
467  }
468  if (Dim == 3)
469  {
470  if (calculate_dxyz)
471  this->dpsideta_map[i].resize (n_qp);
472  if (calculate_d2xyz)
473  {
474  this->d2psidxideta_map[i].resize (n_qp);
475  this->d2psideta2_map[i].resize (n_qp);
476  }
477  }
478 
479  // Compute the value of shape function i, and its first and
480  // second derivatives at quadrature point p
481  // (Lagrange shape functions are used for the mapping)
482  for (unsigned int p=0; p<n_qp; p++)
483  {
484  if (calculate_xyz)
485  this->psi_map[i][p] = FE<Dim-1,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
486  if (Dim > 1)
487  {
488  if (calculate_dxyz)
489  this->dpsidxi_map[i][p] = FE<Dim-1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
490  if (calculate_d2xyz)
491  this->d2psidxi2_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 0, qp[p]);
492  }
493  // libMesh::out << "this->d2psidxi2_map["<<i<<"][p]=" << d2psidxi2_map[i][p] << std::endl;
494 
495  // If we are in 3D, then our sides are 2D faces.
496  // For the second derivatives, we must also compute the cross
497  // derivative d^2() / dxi deta
498  if (Dim == 3)
499  {
500  if (calculate_dxyz)
501  this->dpsideta_map[i][p] = FE<Dim-1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
502  if (calculate_d2xyz)
503  {
504  this->d2psidxideta_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 1, qp[p]);
505  this->d2psideta2_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 2, qp[p]);
506  }
507  }
508  }
509  }
510 }
511 
512 template<unsigned int Dim>
513 void FEMap::init_edge_shape_functions(const std::vector<Point> & qp,
514  const Elem * edge)
515 {
516  // Start logging the shape function initialization
517  LOG_SCOPE("init_edge_shape_functions()", "FEMap");
518 
519  libmesh_assert(edge);
520 
521  // We're calculating now!
522  this->determine_calculations();
523 
524  // The element type and order to use in
525  // the map
526  const Order mapping_order (edge->default_order());
527  const ElemType mapping_elem_type (edge->type());
528 
529  // The number of quadrature points.
530  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
531 
532  const unsigned int n_mapping_shape_functions =
533  FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
534  mapping_order);
535 
536  // resize the vectors to hold current data
537  // Psi are the shape functions used for the FE mapping
538  if (calculate_xyz)
539  this->psi_map.resize (n_mapping_shape_functions);
540  if (calculate_dxyz)
541  this->dpsidxi_map.resize (n_mapping_shape_functions);
542  if (calculate_d2xyz)
543  this->d2psidxi2_map.resize (n_mapping_shape_functions);
544 
545  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
546  {
547  // Allocate space to store the values of the shape functions
548  // and their first and second derivatives at the quadrature points.
549  if (calculate_xyz)
550  this->psi_map[i].resize (n_qp);
551  if (calculate_dxyz)
552  this->dpsidxi_map[i].resize (n_qp);
553  if (calculate_d2xyz)
554  this->d2psidxi2_map[i].resize (n_qp);
555 
556  // Compute the value of shape function i, and its first and
557  // second derivatives at quadrature point p
558  // (Lagrange shape functions are used for the mapping)
559  for (unsigned int p=0; p<n_qp; p++)
560  {
561  if (calculate_xyz)
562  this->psi_map[i][p] = FE<1,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
563  if (calculate_dxyz)
564  this->dpsidxi_map[i][p] = FE<1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
565  if (calculate_d2xyz)
566  this->d2psidxi2_map[i][p] = FE<1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 0, qp[p]);
567  }
568  }
569 }
570 
571 
572 
573 void FEMap::compute_face_map(int dim, const std::vector<Real> & qw,
574  const Elem * side)
575 {
576  libmesh_assert(side);
577 
578  // We're calculating now!
579  this->determine_calculations();
580 
581  LOG_SCOPE("compute_face_map()", "FEMap");
582 
583  // The number of quadrature points.
584  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
585 
586  switch (dim)
587  {
588  case 1:
589  {
590  // A 1D finite element, currently assumed to be in 1D space
591  // This means the boundary is a "0D finite element", a
592  // NODEELEM.
593 
594  // Resize the vectors to hold data at the quadrature points
595  {
596  if (calculate_xyz)
597  this->xyz.resize(n_qp);
598  if (calculate_dxyz)
599  normals.resize(n_qp);
600 
601  if (calculate_dxyz)
602  this->JxW.resize(n_qp);
603  }
604 
605  // If we have no quadrature points, there's nothing else to do
606  if (!n_qp)
607  break;
608 
609  // We need to look back at the full edge to figure out the normal
610  // vector
611  const Elem * elem = side->parent();
612  libmesh_assert (elem);
613  if (calculate_dxyz)
614  {
615  if (side->node_id(0) == elem->node_id(0))
616  normals[0] = Point(-1.);
617  else
618  {
619  libmesh_assert_equal_to (side->node_id(0),
620  elem->node_id(1));
621  normals[0] = Point(1.);
622  }
623  }
624 
625  // Calculate x at the point
626  if (calculate_xyz)
627  libmesh_assert_equal_to (this->psi_map.size(), 1);
628  // In the unlikely event we have multiple quadrature
629  // points, they'll be in the same place
630  for (unsigned int p=0; p<n_qp; p++)
631  {
632  if (calculate_xyz)
633  {
634  this->xyz[p].zero();
635  this->xyz[p].add_scaled (side->point(0), this->psi_map[0][p]);
636  }
637  if (calculate_dxyz)
638  {
639  normals[p] = normals[0];
640  this->JxW[p] = 1.0*qw[p];
641  }
642  }
643 
644  // done computing the map
645  break;
646  }
647 
648  case 2:
649  {
650  // A 2D finite element living in either 2D or 3D space.
651  // This means the boundary is a 1D finite element, i.e.
652  // and EDGE2 or EDGE3.
653  // Resize the vectors to hold data at the quadrature points
654  {
655  if (calculate_xyz)
656  this->xyz.resize(n_qp);
657  if (calculate_dxyz)
658  {
659  this->dxyzdxi_map.resize(n_qp);
660  this->tangents.resize(n_qp);
661  this->normals.resize(n_qp);
662 
663  this->JxW.resize(n_qp);
664  }
665 
666  if (calculate_d2xyz)
667  {
668  this->d2xyzdxi2_map.resize(n_qp);
669  this->curvatures.resize(n_qp);
670  }
671  }
672 
673  // Clear the entities that will be summed
674  // Compute the tangent & normal at the quadrature point
675  for (unsigned int p=0; p<n_qp; p++)
676  {
677  if (calculate_xyz)
678  this->xyz[p].zero();
679  if (calculate_dxyz)
680  {
681  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
682  this->dxyzdxi_map[p].zero();
683  }
684  if (calculate_d2xyz)
685  this->d2xyzdxi2_map[p].zero();
686  }
687 
688  const unsigned int n_mapping_shape_functions =
690  side->default_order());
691 
692  // compute x, dxdxi at the quadrature points
693  for (unsigned int i=0; i<n_mapping_shape_functions; i++) // sum over the nodes
694  {
695  const Point & side_point = side->point(i);
696 
697  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
698  {
699  if (calculate_xyz)
700  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
701  if (calculate_dxyz)
702  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
703  if (calculate_d2xyz)
704  this->d2xyzdxi2_map[p].add_scaled(side_point, this->d2psidxi2_map[i][p]);
705  }
706  }
707 
708  // Compute the tangent & normal at the quadrature point
709  if (calculate_dxyz)
710  {
711  for (unsigned int p=0; p<n_qp; p++)
712  {
713  // The first tangent comes from just the edge's Jacobian
714  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
715 
716 #if LIBMESH_DIM == 2
717  // For a 2D element living in 2D, the normal is given directly
718  // from the entries in the edge Jacobian.
719  this->normals[p] = (Point(this->dxyzdxi_map[p](1), -this->dxyzdxi_map[p](0), 0.)).unit();
720 
721 #elif LIBMESH_DIM == 3
722  // For a 2D element living in 3D, there is a second tangent.
723  // For the second tangent, we need to refer to the full
724  // element's (not just the edge's) Jacobian.
725  const Elem * elem = side->parent();
726  libmesh_assert(elem);
727 
728  // Inverse map xyz[p] to a reference point on the parent...
729  Point reference_point = FE<2,LAGRANGE>::inverse_map(elem, this->xyz[p]);
730 
731  // Get dxyz/dxi and dxyz/deta from the parent map.
732  Point dx_dxi = FE<2,LAGRANGE>::map_xi (elem, reference_point);
733  Point dx_deta = FE<2,LAGRANGE>::map_eta(elem, reference_point);
734 
735  // The second tangent vector is formed by crossing these vectors.
736  tangents[p][1] = dx_dxi.cross(dx_deta).unit();
737 
738  // Finally, the normal in this case is given by crossing these
739  // two tangents.
740  normals[p] = tangents[p][0].cross(tangents[p][1]).unit();
741 #endif
742 
743 
744  // The curvature is computed via the familiar Frenet formula:
745  // curvature = [d^2(x) / d (xi)^2] dot [normal]
746  // For a reference, see:
747  // F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill, p. 310
748  //
749  // Note: The sign convention here is different from the
750  // 3D case. Concave-upward curves (smiles) have a positive
751  // curvature. Concave-downward curves (frowns) have a
752  // negative curvature. Be sure to take that into account!
753  if (calculate_d2xyz)
754  {
755  const Real numerator = this->d2xyzdxi2_map[p] * this->normals[p];
756  const Real denominator = this->dxyzdxi_map[p].norm_sq();
757  libmesh_assert_not_equal_to (denominator, 0);
758  curvatures[p] = numerator / denominator;
759  }
760  }
761 
762  // compute the jacobian at the quadrature points
763  for (unsigned int p=0; p<n_qp; p++)
764  {
765  const Real the_jac = this->dxyzdxi_map[p].norm();
766 
767  libmesh_assert_greater (the_jac, 0.);
768 
769  this->JxW[p] = the_jac*qw[p];
770  }
771  }
772 
773  // done computing the map
774  break;
775  }
776 
777 
778 
779  case 3:
780  {
781  // A 3D finite element living in 3D space.
782  // Resize the vectors to hold data at the quadrature points
783  {
784  if (calculate_xyz)
785  this->xyz.resize(n_qp);
786  if (calculate_dxyz)
787  {
788  this->dxyzdxi_map.resize(n_qp);
789  this->dxyzdeta_map.resize(n_qp);
790  this->tangents.resize(n_qp);
791  this->normals.resize(n_qp);
792  this->JxW.resize(n_qp);
793  }
794  if (calculate_d2xyz)
795  {
796  this->d2xyzdxi2_map.resize(n_qp);
797  this->d2xyzdxideta_map.resize(n_qp);
798  this->d2xyzdeta2_map.resize(n_qp);
799  this->curvatures.resize(n_qp);
800  }
801  }
802 
803  // Clear the entities that will be summed
804  for (unsigned int p=0; p<n_qp; p++)
805  {
806  if (calculate_xyz)
807  this->xyz[p].zero();
808  if (calculate_dxyz)
809  {
810  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
811  this->dxyzdxi_map[p].zero();
812  this->dxyzdeta_map[p].zero();
813  }
814  if (calculate_d2xyz)
815  {
816  this->d2xyzdxi2_map[p].zero();
817  this->d2xyzdxideta_map[p].zero();
818  this->d2xyzdeta2_map[p].zero();
819  }
820  }
821 
822  const unsigned int n_mapping_shape_functions =
824  side->default_order());
825 
826  // compute x, dxdxi at the quadrature points
827  for (unsigned int i=0; i<n_mapping_shape_functions; i++) // sum over the nodes
828  {
829  const Point & side_point = side->point(i);
830 
831  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
832  {
833  if (calculate_xyz)
834  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
835  if (calculate_dxyz)
836  {
837  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
838  this->dxyzdeta_map[p].add_scaled(side_point, this->dpsideta_map[i][p]);
839  }
840  if (calculate_d2xyz)
841  {
842  this->d2xyzdxi2_map[p].add_scaled (side_point, this->d2psidxi2_map[i][p]);
843  this->d2xyzdxideta_map[p].add_scaled(side_point, this->d2psidxideta_map[i][p]);
844  this->d2xyzdeta2_map[p].add_scaled (side_point, this->d2psideta2_map[i][p]);
845  }
846  }
847  }
848 
849  // Compute the tangents, normal, and curvature at the quadrature point
850  if (calculate_dxyz)
851  {
852  for (unsigned int p=0; p<n_qp; p++)
853  {
854  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
855  this->normals[p] = n.unit();
856  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
857  this->tangents[p][1] = n.cross(this->dxyzdxi_map[p]).unit();
858 
859  if (calculate_d2xyz)
860  {
861  // Compute curvature using the typical nomenclature
862  // of the first and second fundamental forms.
863  // For reference, see:
864  // 1) http://mathworld.wolfram.com/MeanCurvature.html
865  // (note -- they are using inward normal)
866  // 2) F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill
867  const Real L = -this->d2xyzdxi2_map[p] * this->normals[p];
868  const Real M = -this->d2xyzdxideta_map[p] * this->normals[p];
869  const Real N = -this->d2xyzdeta2_map[p] * this->normals[p];
870  const Real E = this->dxyzdxi_map[p].norm_sq();
871  const Real F = this->dxyzdxi_map[p] * this->dxyzdeta_map[p];
872  const Real G = this->dxyzdeta_map[p].norm_sq();
873 
874  const Real numerator = E*N -2.*F*M + G*L;
875  const Real denominator = E*G - F*F;
876  libmesh_assert_not_equal_to (denominator, 0.);
877  curvatures[p] = 0.5*numerator/denominator;
878  }
879  }
880 
881  // compute the jacobian at the quadrature points, see
882  // http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/thtoc.htm
883  for (unsigned int p=0; p<n_qp; p++)
884  {
885  const Real g11 = (dxdxi_map(p)*dxdxi_map(p) +
886  dydxi_map(p)*dydxi_map(p) +
887  dzdxi_map(p)*dzdxi_map(p));
888 
889  const Real g12 = (dxdxi_map(p)*dxdeta_map(p) +
890  dydxi_map(p)*dydeta_map(p) +
891  dzdxi_map(p)*dzdeta_map(p));
892 
893  const Real g21 = g12;
894 
895  const Real g22 = (dxdeta_map(p)*dxdeta_map(p) +
896  dydeta_map(p)*dydeta_map(p) +
897  dzdeta_map(p)*dzdeta_map(p));
898 
899 
900  const Real the_jac = std::sqrt(g11*g22 - g12*g21);
901 
902  libmesh_assert_greater (the_jac, 0.);
903 
904  this->JxW[p] = the_jac*qw[p];
905  }
906  }
907 
908  // done computing the map
909  break;
910  }
911 
912 
913  default:
914  libmesh_error_msg("Invalid dimension dim = " << dim);
915  }
916 }
917 
918 
919 
920 
922  const std::vector<Real> & qw,
923  const Elem * edge)
924 {
925  libmesh_assert(edge);
926 
927  if (dim == 2)
928  {
929  // A 2D finite element living in either 2D or 3D space.
930  // The edges here are the sides of the element, so the
931  // (misnamed) compute_face_map function does what we want
932  this->compute_face_map(dim, qw, edge);
933  return;
934  }
935 
936  libmesh_assert_equal_to (dim, 3); // 1D is unnecessary and currently unsupported
937 
938  LOG_SCOPE("compute_edge_map()", "FEMap");
939 
940  // We're calculating now!
941  this->determine_calculations();
942 
943  // The number of quadrature points.
944  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
945 
946  // Resize the vectors to hold data at the quadrature points
947  if (calculate_xyz)
948  this->xyz.resize(n_qp);
949  if (calculate_dxyz)
950  {
951  this->dxyzdxi_map.resize(n_qp);
952  this->dxyzdeta_map.resize(n_qp);
953  this->tangents.resize(n_qp);
954  this->normals.resize(n_qp);
955  this->JxW.resize(n_qp);
956  }
957  if (calculate_d2xyz)
958  {
959  this->d2xyzdxi2_map.resize(n_qp);
960  this->d2xyzdxideta_map.resize(n_qp);
961  this->d2xyzdeta2_map.resize(n_qp);
962  this->curvatures.resize(n_qp);
963  }
964 
965  // Clear the entities that will be summed
966  for (unsigned int p=0; p<n_qp; p++)
967  {
968  if (calculate_xyz)
969  this->xyz[p].zero();
970  if (calculate_dxyz)
971  {
972  this->tangents[p].resize(1);
973  this->dxyzdxi_map[p].zero();
974  this->dxyzdeta_map[p].zero();
975  }
976  if (calculate_d2xyz)
977  {
978  this->d2xyzdxi2_map[p].zero();
979  this->d2xyzdxideta_map[p].zero();
980  this->d2xyzdeta2_map[p].zero();
981  }
982  }
983 
984  // compute x, dxdxi at the quadrature points
985  for (std::size_t i=0; i<this->psi_map.size(); i++) // sum over the nodes
986  {
987  const Point & edge_point = edge->point(i);
988 
989  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
990  {
991  if (calculate_xyz)
992  this->xyz[p].add_scaled (edge_point, this->psi_map[i][p]);
993  if (calculate_dxyz)
994  this->dxyzdxi_map[p].add_scaled (edge_point, this->dpsidxi_map[i][p]);
995  if (calculate_d2xyz)
996  this->d2xyzdxi2_map[p].add_scaled (edge_point, this->d2psidxi2_map[i][p]);
997  }
998  }
999 
1000  // Compute the tangents at the quadrature point
1001  // FIXME: normals (plural!) and curvatures are uncalculated
1002  if (calculate_dxyz)
1003  for (unsigned int p=0; p<n_qp; p++)
1004  {
1005  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
1006  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
1007 
1008  // compute the jacobian at the quadrature points
1009  const Real the_jac = std::sqrt(this->dxdxi_map(p)*this->dxdxi_map(p) +
1010  this->dydxi_map(p)*this->dydxi_map(p) +
1011  this->dzdxi_map(p)*this->dzdxi_map(p));
1012 
1013  libmesh_assert_greater (the_jac, 0.);
1014 
1015  this->JxW[p] = the_jac*qw[p];
1016  }
1017 }
1018 
1019 
1020 // Explicit FEMap Instantiations
1021 FACE_EDGE_SHAPE_ERROR(0,init_face_shape_functions)
1022 template void FEMap::init_face_shape_functions<1>(const std::vector<Point> &, const Elem *);
1023 template void FEMap::init_face_shape_functions<2>(const std::vector<Point> &, const Elem *);
1024 template void FEMap::init_face_shape_functions<3>(const std::vector<Point> &, const Elem *);
1025 
1026 FACE_EDGE_SHAPE_ERROR(0,init_edge_shape_functions)
1027 template void FEMap::init_edge_shape_functions<1>(const std::vector<Point> &, const Elem *);
1028 template void FEMap::init_edge_shape_functions<2>(const std::vector<Point> &, const Elem *);
1029 template void FEMap::init_edge_shape_functions<3>(const std::vector<Point> &, const Elem *);
1030 
1031 //--------------------------------------------------------------
1032 // Explicit FE instantiations
1033 template void FE<1,LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1034 template void FE<1,LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1035 template void FE<1,LAGRANGE_VEC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1036 template void FE<1,LAGRANGE_VEC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1037 template void FE<1,L2_LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1038 template void FE<1,L2_LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1039 template void FE<1,HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1040 template void FE<1,HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1041 template void FE<1,L2_HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1042 template void FE<1,L2_HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1043 template void FE<1,CLOUGH>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1044 template void FE<1,CLOUGH>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1045 template void FE<1,HERMITE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1046 template void FE<1,HERMITE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1047 template void FE<1,MONOMIAL>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1048 template void FE<1,MONOMIAL>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1049 template void FE<1,SCALAR>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1050 template void FE<1,SCALAR>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1051 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1052 template void FE<1,BERNSTEIN>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1053 template void FE<1,BERNSTEIN>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1054 template void FE<1,SZABAB>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1055 template void FE<1,SZABAB>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1056 #endif
1057 template void FE<1,XYZ>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1058 template void FE<1,XYZ>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1059 
1060 template void FE<2,LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1061 template void FE<2,LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1062 template void FE<2,LAGRANGE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1063 template void FE<2,LAGRANGE_VEC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1064 template void FE<2,LAGRANGE_VEC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1065 template void FE<2,LAGRANGE_VEC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1066 template void FE<2,L2_LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1067 template void FE<2,L2_LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1068 template void FE<2,L2_LAGRANGE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1069 template void FE<2,HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1070 template void FE<2,HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1071 template void FE<2,HIERARCHIC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1072 template void FE<2,L2_HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1073 template void FE<2,L2_HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1074 template void FE<2,L2_HIERARCHIC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1075 template void FE<2,CLOUGH>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1076 template void FE<2,CLOUGH>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1077 template void FE<2,CLOUGH>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1078 template void FE<2,HERMITE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1079 template void FE<2,HERMITE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1080 template void FE<2,HERMITE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1081 template void FE<2,MONOMIAL>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1082 template void FE<2,MONOMIAL>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1083 template void FE<2,MONOMIAL>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1084 template void FE<2,SCALAR>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1085 template void FE<2,SCALAR>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1086 template void FE<2,SCALAR>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1087 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1088 template void FE<2,BERNSTEIN>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1089 template void FE<2,BERNSTEIN>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1090 template void FE<2,BERNSTEIN>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1091 template void FE<2,SZABAB>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1092 template void FE<2,SZABAB>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1093 template void FE<2,SZABAB>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1094 #endif
1095 template void FE<2,SUBDIVISION>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1096 template void FE<2,XYZ>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1097 template void FE<2,XYZ>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1098 template void FE<2,XYZ>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1099 template void FE<2,NEDELEC_ONE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1100 template void FE<2,NEDELEC_ONE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1101 template void FE<2,NEDELEC_ONE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1102 
1103 // Intel 9.1 complained it needed this in devel mode.
1104 //template void FE<2,XYZ>::init_face_shape_functions(const std::vector<Point> &, const Elem *);
1105 
1106 template void FE<3,LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1107 template void FE<3,LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1108 template void FE<3,LAGRANGE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1109 template void FE<3,LAGRANGE_VEC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1110 template void FE<3,LAGRANGE_VEC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1111 template void FE<3,LAGRANGE_VEC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1112 template void FE<3,L2_LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1113 template void FE<3,L2_LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1114 template void FE<3,L2_LAGRANGE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1115 template void FE<3,HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1116 template void FE<3,HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1117 template void FE<3,HIERARCHIC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1118 template void FE<3,L2_HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1119 template void FE<3,L2_HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1120 template void FE<3,L2_HIERARCHIC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1121 template void FE<3,CLOUGH>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1122 template void FE<3,CLOUGH>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1123 template void FE<3,CLOUGH>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1124 template void FE<3,HERMITE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1125 template void FE<3,HERMITE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1126 template void FE<3,HERMITE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1127 template void FE<3,MONOMIAL>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1128 template void FE<3,MONOMIAL>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1129 template void FE<3,MONOMIAL>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1130 template void FE<3,SCALAR>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1131 template void FE<3,SCALAR>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1132 template void FE<3,SCALAR>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1133 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1134 template void FE<3,BERNSTEIN>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1135 template void FE<3,BERNSTEIN>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1136 template void FE<3,BERNSTEIN>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1137 template void FE<3,SZABAB>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1138 template void FE<3,SZABAB>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1139 template void FE<3,SZABAB>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1140 #endif
1141 template void FE<3,XYZ>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1142 template void FE<3,XYZ>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1143 template void FE<3,XYZ>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1144 template void FE<3,NEDELEC_ONE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1145 template void FE<3,NEDELEC_ONE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1146 template void FE<3,NEDELEC_ONE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1147 
1148 // Intel 9.1 complained it needed this in devel mode.
1149 //template void FE<3,XYZ>::init_face_shape_functions(const std::vector<Point> &, const Elem *);
1150 
1151 } // namespace libMesh
void init_edge_shape_functions(const std::vector< Point > &qp, const Elem *edge)
Same as before, but for an edge.
Definition: fe_boundary.C:513
virtual void edge_reinit(const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
Reinitializes all the physical element-dependent data based on the edge.
Definition: fe_boundary.C:260
virtual void side_map(const Elem *elem, const Elem *side, const unsigned int s, const std::vector< Point > &reference_side_points, std::vector< Point > &reference_points) libmesh_override
Computes the reference space quadrature points on the side of an element based on the side quadrature...
Definition: fe_boundary.C:353
std::vector< std::vector< Real > > dpsideta_map
Map for the derivative of the side function, d(psi)/d(eta).
Definition: fe_map.h:820
std::vector< std::vector< Real > > psi_map
Map for the side shape functions, psi.
Definition: fe_map.h:808
virtual UniquePtr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
std::vector< std::vector< Real > > d2psideta2_map
Map for the second derivatives (in eta) of the side shape functions.
Definition: fe_map.h:841
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
virtual ElemType type() const =0
unsigned int p_level() const
Definition: elem.h:2422
bool calculate_dxyz
Should we calculate mapping gradients?
Definition: fe_map.h:884
static Point map_xi(const Elem *elem, const Point &reference_point)
Definition: fe_map.C:2002
void compute_edge_map(int dim, const std::vector< Real > &qw, const Elem *side)
Same as before, but for an edge.
Definition: fe_boundary.C:921
SIDEMAP_ERROR(0, CLOUGH, side_map)
Definition: fe_boundary.C:67
unsigned int dim
REINIT_ERROR(0, CLOUGH, reinit)
Definition: fe_boundary.C:65
ElemType
Defines an enum for geometric element types.
const Elem * parent() const
Definition: elem.h:2346
Real dxdxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:546
unsigned short int side
Definition: xdr_io.C:49
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
const class libmesh_nullptr_t libmesh_nullptr
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
The libMesh namespace provides an interface to certain functionality in the library.
Real dydxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:554
const Number zero
.
Definition: libmesh.h:178
long double max(long double a, double b)
std::vector< std::vector< Real > > d2psidxideta_map
Map for the second (cross) derivatives in xi, eta of the side shape functions.
Definition: fe_map.h:834
std::vector< RealGradient > d2xyzdxideta_map
Vector of mixed second partial derivatives in xi-eta: d^2(x)/d(xi)d(eta) d^2(y)/d(xi)d(eta) d^2(z)/d(...
Definition: fe_map.h:645
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
virtual unsigned int n_nodes() const =0
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:1967
std::vector< RealGradient > dxyzdxi_map
Vector of partial derivatives: d(x)/d(xi), d(y)/d(xi), d(z)/d(xi)
Definition: fe_map.h:621
std::vector< RealGradient > d2xyzdeta2_map
Vector of second partial derivatives in eta: d^2(x)/d(eta)^2.
Definition: fe_map.h:651
std::vector< RealGradient > d2xyzdxi2_map
Vector of second partial derivatives in xi: d^2(x)/d(xi)^2, d^2(y)/d(xi)^2, d^2(z)/d(xi)^2.
Definition: fe_map.h:639
bool calculate_d2xyz
Should we calculate mapping hessians?
Definition: fe_map.h:889
Real dzdxi_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:562
bool calculate_xyz
Should we calculate physical point locations?
Definition: fe_map.h:879
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition: type_vector.h:879
virtual Order default_order() const =0
Real dxdeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:570
Real dydeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:578
std::vector< Real > curvatures
The mean curvature (= one half the sum of the principal curvatures) on the boundary at the quadrature...
Definition: fe_map.h:858
std::vector< std::vector< Real > > d2psidxi2_map
Map for the second derivatives (in xi) of the side shape functions.
Definition: fe_map.h:827
TypeVector< T > unit() const
Definition: type_vector.C:37
std::vector< Real > JxW
Jacobian*Weight values at quadrature points.
Definition: fe_map.h:868
std::vector< Point > normals
Normal vectors on boundary at quadrature points.
Definition: fe_map.h:851
std::vector< Point > xyz
The spatial locations of the quadrature points.
Definition: fe_map.h:615
static Point map_eta(const Elem *elem, const Point &reference_point)
Definition: fe_map.C:2030
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=libmesh_nullptr, const std::vector< Real > *const weights=libmesh_nullptr) libmesh_override
This is at the core of this class.
Definition: fe.C:115
std::vector< std::vector< Real > > dpsidxi_map
Map for the derivative of the side functions, d(psi)/d(xi).
Definition: fe_map.h:814
virtual unsigned int n_shape_functions() const libmesh_override
Definition: fe.C:36
std::vector< std::vector< Point > > tangents
Tangent vectors on boundary at quadrature points.
Definition: fe_map.h:846
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
std::vector< RealGradient > dxyzdeta_map
Vector of partial derivatives: d(x)/d(eta), d(y)/d(eta), d(z)/d(eta)
Definition: fe_map.h:627
virtual UniquePtr< Elem > build_edge_ptr(const unsigned int i)=0
void determine_calculations()
Determine which values are to be calculated.
Definition: fe_map.h:524
virtual void compute_face_map(int dim, const std::vector< Real > &qw, const Elem *side)
Same as compute_map, but for a side.
Definition: fe_boundary.C:573
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1831
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
void init_face_shape_functions(const std::vector< Point > &qp, const Elem *side)
Initializes the reference to physical element map for a side.
Definition: fe_boundary.C:408
Real dzdeta_map(const unsigned int p) const
Used in FEMap::compute_map(), which should be be usable in derived classes, and therefore protected...
Definition: fe_map.h:586
static Point inverse_map(const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_map.C:1566