fe_boundary.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 
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 != nullptr || pts != 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 std::unique_ptr<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) != 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 != 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 != 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 != 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 != nullptr || pts != 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 std::unique_ptr<const Elem> edge(elem->build_edge_ptr(e));
281 
282  // Initialize the shape functions at the user-specified
283  // points
284  if (pts != 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 != 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) != 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  // sum over the nodes
400  for (auto i : index_range(psi_map))
401  {
402  const Point & side_node = refspace_nodes[elem_nodes_map[i]];
403  for (unsigned int p=0; p<n_points; p++)
404  reference_points[p].add_scaled (side_node, psi_map[i][p]);
405  }
406 }
407 
408 template<unsigned int Dim>
409 void FEMap::init_face_shape_functions(const std::vector<Point> & qp,
410  const Elem * side)
411 {
412  // Start logging the shape function initialization
413  LOG_SCOPE("init_face_shape_functions()", "FEMap");
414 
415  libmesh_assert(side);
416 
417  // We're calculating now!
418  this->determine_calculations();
419 
420  // The element type and order to use in
421  // the map
422  const Order mapping_order (side->default_order());
423  const ElemType mapping_elem_type (side->type());
424 
425  // The number of quadrature points.
426  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
427 
428  const unsigned int n_mapping_shape_functions =
429  FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
430  mapping_order);
431 
432  // resize the vectors to hold current data
433  // Psi are the shape functions used for the FE mapping
434  if (calculate_xyz)
435  this->psi_map.resize (n_mapping_shape_functions);
436 
437  if (Dim > 1)
438  {
439  if (calculate_dxyz)
440  this->dpsidxi_map.resize (n_mapping_shape_functions);
441  if (calculate_d2xyz)
442  this->d2psidxi2_map.resize (n_mapping_shape_functions);
443  }
444 
445  if (Dim == 3)
446  {
447  if (calculate_dxyz)
448  this->dpsideta_map.resize (n_mapping_shape_functions);
449  if (calculate_d2xyz)
450  {
451  this->d2psidxideta_map.resize (n_mapping_shape_functions);
452  this->d2psideta2_map.resize (n_mapping_shape_functions);
453  }
454  }
455 
456  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
457  {
458  // Allocate space to store the values of the shape functions
459  // and their first and second derivatives at the quadrature points.
460  if (calculate_xyz)
461  this->psi_map[i].resize (n_qp);
462  if (Dim > 1)
463  {
464  if (calculate_dxyz)
465  this->dpsidxi_map[i].resize (n_qp);
466  if (calculate_d2xyz)
467  this->d2psidxi2_map[i].resize (n_qp);
468  }
469  if (Dim == 3)
470  {
471  if (calculate_dxyz)
472  this->dpsideta_map[i].resize (n_qp);
473  if (calculate_d2xyz)
474  {
475  this->d2psidxideta_map[i].resize (n_qp);
476  this->d2psideta2_map[i].resize (n_qp);
477  }
478  }
479 
480  // Compute the value of shape function i, and its first and
481  // second derivatives at quadrature point p
482  // (Lagrange shape functions are used for the mapping)
483  for (unsigned int p=0; p<n_qp; p++)
484  {
485  if (calculate_xyz)
486  this->psi_map[i][p] = FE<Dim-1,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
487  if (Dim > 1)
488  {
489  if (calculate_dxyz)
490  this->dpsidxi_map[i][p] = FE<Dim-1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
491  if (calculate_d2xyz)
492  this->d2psidxi2_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 0, qp[p]);
493  }
494  // libMesh::out << "this->d2psidxi2_map["<<i<<"][p]=" << d2psidxi2_map[i][p] << std::endl;
495 
496  // If we are in 3D, then our sides are 2D faces.
497  // For the second derivatives, we must also compute the cross
498  // derivative d^2() / dxi deta
499  if (Dim == 3)
500  {
501  if (calculate_dxyz)
502  this->dpsideta_map[i][p] = FE<Dim-1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
503  if (calculate_d2xyz)
504  {
505  this->d2psidxideta_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 1, qp[p]);
506  this->d2psideta2_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 2, qp[p]);
507  }
508  }
509  }
510  }
511 }
512 
513 template<unsigned int Dim>
514 void FEMap::init_edge_shape_functions(const std::vector<Point> & qp,
515  const Elem * edge)
516 {
517  // Start logging the shape function initialization
518  LOG_SCOPE("init_edge_shape_functions()", "FEMap");
519 
520  libmesh_assert(edge);
521 
522  // We're calculating now!
523  this->determine_calculations();
524 
525  // The element type and order to use in
526  // the map
527  const Order mapping_order (edge->default_order());
528  const ElemType mapping_elem_type (edge->type());
529 
530  // The number of quadrature points.
531  const unsigned int n_qp = cast_int<unsigned int>(qp.size());
532 
533  const unsigned int n_mapping_shape_functions =
534  FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
535  mapping_order);
536 
537  // resize the vectors to hold current data
538  // Psi are the shape functions used for the FE mapping
539  if (calculate_xyz)
540  this->psi_map.resize (n_mapping_shape_functions);
541  if (calculate_dxyz)
542  this->dpsidxi_map.resize (n_mapping_shape_functions);
543  if (calculate_d2xyz)
544  this->d2psidxi2_map.resize (n_mapping_shape_functions);
545 
546  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
547  {
548  // Allocate space to store the values of the shape functions
549  // and their first and second derivatives at the quadrature points.
550  if (calculate_xyz)
551  this->psi_map[i].resize (n_qp);
552  if (calculate_dxyz)
553  this->dpsidxi_map[i].resize (n_qp);
554  if (calculate_d2xyz)
555  this->d2psidxi2_map[i].resize (n_qp);
556 
557  // Compute the value of shape function i, and its first and
558  // second derivatives at quadrature point p
559  // (Lagrange shape functions are used for the mapping)
560  for (unsigned int p=0; p<n_qp; p++)
561  {
562  if (calculate_xyz)
563  this->psi_map[i][p] = FE<1,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
564  if (calculate_dxyz)
565  this->dpsidxi_map[i][p] = FE<1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
566  if (calculate_d2xyz)
567  this->d2psidxi2_map[i][p] = FE<1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 0, qp[p]);
568  }
569  }
570 }
571 
572 
573 
574 void FEMap::compute_face_map(int dim, const std::vector<Real> & qw,
575  const Elem * side)
576 {
577  libmesh_assert(side);
578 
579  // We're calculating now!
580  this->determine_calculations();
581 
582  LOG_SCOPE("compute_face_map()", "FEMap");
583 
584  // The number of quadrature points.
585  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
586 
587  switch (dim)
588  {
589  case 1:
590  {
591  // A 1D finite element, currently assumed to be in 1D space
592  // This means the boundary is a "0D finite element", a
593  // NODEELEM.
594 
595  // Resize the vectors to hold data at the quadrature points
596  {
597  if (calculate_xyz)
598  this->xyz.resize(n_qp);
599  if (calculate_dxyz)
600  normals.resize(n_qp);
601 
602  if (calculate_dxyz)
603  this->JxW.resize(n_qp);
604  }
605 
606  // If we have no quadrature points, there's nothing else to do
607  if (!n_qp)
608  break;
609 
610  // We need to look back at the full edge to figure out the normal
611  // vector
612  const Elem * elem = side->parent();
613  libmesh_assert (elem);
614  if (calculate_dxyz)
615  {
616  if (side->node_id(0) == elem->node_id(0))
617  normals[0] = Point(-1.);
618  else
619  {
620  libmesh_assert_equal_to (side->node_id(0),
621  elem->node_id(1));
622  normals[0] = Point(1.);
623  }
624  }
625 
626  // Calculate x at the point
627  if (calculate_xyz)
628  libmesh_assert_equal_to (this->psi_map.size(), 1);
629  // In the unlikely event we have multiple quadrature
630  // points, they'll be in the same place
631  for (unsigned int p=0; p<n_qp; p++)
632  {
633  if (calculate_xyz)
634  {
635  this->xyz[p].zero();
636  this->xyz[p].add_scaled (side->point(0), this->psi_map[0][p]);
637  }
638  if (calculate_dxyz)
639  {
640  normals[p] = normals[0];
641  this->JxW[p] = 1.0*qw[p];
642  }
643  }
644 
645  // done computing the map
646  break;
647  }
648 
649  case 2:
650  {
651  // A 2D finite element living in either 2D or 3D space.
652  // This means the boundary is a 1D finite element, i.e.
653  // and EDGE2 or EDGE3.
654  // Resize the vectors to hold data at the quadrature points
655  {
656  if (calculate_xyz)
657  this->xyz.resize(n_qp);
658  if (calculate_dxyz)
659  {
660  this->dxyzdxi_map.resize(n_qp);
661  this->tangents.resize(n_qp);
662  this->normals.resize(n_qp);
663 
664  this->JxW.resize(n_qp);
665  }
666 
667  if (calculate_d2xyz)
668  {
669  this->d2xyzdxi2_map.resize(n_qp);
670  this->curvatures.resize(n_qp);
671  }
672  }
673 
674  // Clear the entities that will be summed
675  // Compute the tangent & normal at the quadrature point
676  for (unsigned int p=0; p<n_qp; p++)
677  {
678  if (calculate_xyz)
679  this->xyz[p].zero();
680  if (calculate_dxyz)
681  {
682  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
683  this->dxyzdxi_map[p].zero();
684  }
685  if (calculate_d2xyz)
686  this->d2xyzdxi2_map[p].zero();
687  }
688 
689  const unsigned int n_mapping_shape_functions =
690  FE<2,LAGRANGE>::n_shape_functions (side->type(),
691  side->default_order());
692 
693  // compute x, dxdxi at the quadrature points
694  for (unsigned int i=0; i<n_mapping_shape_functions; i++) // sum over the nodes
695  {
696  const Point & side_point = side->point(i);
697 
698  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
699  {
700  if (calculate_xyz)
701  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
702  if (calculate_dxyz)
703  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
704  if (calculate_d2xyz)
705  this->d2xyzdxi2_map[p].add_scaled(side_point, this->d2psidxi2_map[i][p]);
706  }
707  }
708 
709  // Compute the tangent & normal at the quadrature point
710  if (calculate_dxyz)
711  {
712  for (unsigned int p=0; p<n_qp; p++)
713  {
714  // The first tangent comes from just the edge's Jacobian
715  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
716 
717 #if LIBMESH_DIM == 2
718  // For a 2D element living in 2D, the normal is given directly
719  // from the entries in the edge Jacobian.
720  this->normals[p] = (Point(this->dxyzdxi_map[p](1), -this->dxyzdxi_map[p](0), 0.)).unit();
721 
722 #elif LIBMESH_DIM == 3
723  // For a 2D element living in 3D, there is a second tangent.
724  // For the second tangent, we need to refer to the full
725  // element's (not just the edge's) Jacobian.
726  const Elem * elem = side->parent();
727  libmesh_assert(elem);
728 
729  // Inverse map xyz[p] to a reference point on the parent...
730  Point reference_point = FE<2,LAGRANGE>::inverse_map(elem, this->xyz[p]);
731 
732  // Get dxyz/dxi and dxyz/deta from the parent map.
733  Point dx_dxi = FE<2,LAGRANGE>::map_xi (elem, reference_point);
734  Point dx_deta = FE<2,LAGRANGE>::map_eta(elem, reference_point);
735 
736  // The second tangent vector is formed by crossing these vectors.
737  tangents[p][1] = dx_dxi.cross(dx_deta).unit();
738 
739  // Finally, the normal in this case is given by crossing these
740  // two tangents.
741  normals[p] = tangents[p][0].cross(tangents[p][1]).unit();
742 #endif
743 
744 
745  // The curvature is computed via the familiar Frenet formula:
746  // curvature = [d^2(x) / d (xi)^2] dot [normal]
747  // For a reference, see:
748  // F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill, p. 310
749  //
750  // Note: The sign convention here is different from the
751  // 3D case. Concave-upward curves (smiles) have a positive
752  // curvature. Concave-downward curves (frowns) have a
753  // negative curvature. Be sure to take that into account!
754  if (calculate_d2xyz)
755  {
756  const Real numerator = this->d2xyzdxi2_map[p] * this->normals[p];
757  const Real denominator = this->dxyzdxi_map[p].norm_sq();
758  libmesh_assert_not_equal_to (denominator, 0);
759  curvatures[p] = numerator / denominator;
760  }
761  }
762 
763  // compute the jacobian at the quadrature points
764  for (unsigned int p=0; p<n_qp; p++)
765  {
766  const Real the_jac = this->dxyzdxi_map[p].norm();
767 
768  libmesh_assert_greater (the_jac, 0.);
769 
770  this->JxW[p] = the_jac*qw[p];
771  }
772  }
773 
774  // done computing the map
775  break;
776  }
777 
778 
779 
780  case 3:
781  {
782  // A 3D finite element living in 3D space.
783  // Resize the vectors to hold data at the quadrature points
784  {
785  if (calculate_xyz)
786  this->xyz.resize(n_qp);
787  if (calculate_dxyz)
788  {
789  this->dxyzdxi_map.resize(n_qp);
790  this->dxyzdeta_map.resize(n_qp);
791  this->tangents.resize(n_qp);
792  this->normals.resize(n_qp);
793  this->JxW.resize(n_qp);
794  }
795  if (calculate_d2xyz)
796  {
797  this->d2xyzdxi2_map.resize(n_qp);
798  this->d2xyzdxideta_map.resize(n_qp);
799  this->d2xyzdeta2_map.resize(n_qp);
800  this->curvatures.resize(n_qp);
801  }
802  }
803 
804  // Clear the entities that will be summed
805  for (unsigned int p=0; p<n_qp; p++)
806  {
807  if (calculate_xyz)
808  this->xyz[p].zero();
809  if (calculate_dxyz)
810  {
811  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
812  this->dxyzdxi_map[p].zero();
813  this->dxyzdeta_map[p].zero();
814  }
815  if (calculate_d2xyz)
816  {
817  this->d2xyzdxi2_map[p].zero();
818  this->d2xyzdxideta_map[p].zero();
819  this->d2xyzdeta2_map[p].zero();
820  }
821  }
822 
823  const unsigned int n_mapping_shape_functions =
824  FE<3,LAGRANGE>::n_shape_functions (side->type(),
825  side->default_order());
826 
827  // compute x, dxdxi at the quadrature points
828  for (unsigned int i=0; i<n_mapping_shape_functions; i++) // sum over the nodes
829  {
830  const Point & side_point = side->point(i);
831 
832  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
833  {
834  if (calculate_xyz)
835  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
836  if (calculate_dxyz)
837  {
838  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
839  this->dxyzdeta_map[p].add_scaled(side_point, this->dpsideta_map[i][p]);
840  }
841  if (calculate_d2xyz)
842  {
843  this->d2xyzdxi2_map[p].add_scaled (side_point, this->d2psidxi2_map[i][p]);
844  this->d2xyzdxideta_map[p].add_scaled(side_point, this->d2psidxideta_map[i][p]);
845  this->d2xyzdeta2_map[p].add_scaled (side_point, this->d2psideta2_map[i][p]);
846  }
847  }
848  }
849 
850  // Compute the tangents, normal, and curvature at the quadrature point
851  if (calculate_dxyz)
852  {
853  for (unsigned int p=0; p<n_qp; p++)
854  {
855  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
856  this->normals[p] = n.unit();
857  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
858  this->tangents[p][1] = n.cross(this->dxyzdxi_map[p]).unit();
859 
860  if (calculate_d2xyz)
861  {
862  // Compute curvature using the typical nomenclature
863  // of the first and second fundamental forms.
864  // For reference, see:
865  // 1) http://mathworld.wolfram.com/MeanCurvature.html
866  // (note -- they are using inward normal)
867  // 2) F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill
868  const Real L = -this->d2xyzdxi2_map[p] * this->normals[p];
869  const Real M = -this->d2xyzdxideta_map[p] * this->normals[p];
870  const Real N = -this->d2xyzdeta2_map[p] * this->normals[p];
871  const Real E = this->dxyzdxi_map[p].norm_sq();
872  const Real F = this->dxyzdxi_map[p] * this->dxyzdeta_map[p];
873  const Real G = this->dxyzdeta_map[p].norm_sq();
874 
875  const Real numerator = E*N -2.*F*M + G*L;
876  const Real denominator = E*G - F*F;
877  libmesh_assert_not_equal_to (denominator, 0.);
878  curvatures[p] = 0.5*numerator/denominator;
879  }
880  }
881 
882  // compute the jacobian at the quadrature points, see
883  // http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/thtoc.htm
884  for (unsigned int p=0; p<n_qp; p++)
885  {
886  const Real g11 = (dxdxi_map(p)*dxdxi_map(p) +
887  dydxi_map(p)*dydxi_map(p) +
888  dzdxi_map(p)*dzdxi_map(p));
889 
890  const Real g12 = (dxdxi_map(p)*dxdeta_map(p) +
891  dydxi_map(p)*dydeta_map(p) +
892  dzdxi_map(p)*dzdeta_map(p));
893 
894  const Real g21 = g12;
895 
896  const Real g22 = (dxdeta_map(p)*dxdeta_map(p) +
897  dydeta_map(p)*dydeta_map(p) +
898  dzdeta_map(p)*dzdeta_map(p));
899 
900 
901  const Real the_jac = std::sqrt(g11*g22 - g12*g21);
902 
903  libmesh_assert_greater (the_jac, 0.);
904 
905  this->JxW[p] = the_jac*qw[p];
906  }
907  }
908 
909  // done computing the map
910  break;
911  }
912 
913 
914  default:
915  libmesh_error_msg("Invalid dimension dim = " << dim);
916  }
917 }
918 
919 
920 
921 
922 void FEMap::compute_edge_map(int dim,
923  const std::vector<Real> & qw,
924  const Elem * edge)
925 {
926  libmesh_assert(edge);
927 
928  if (dim == 2)
929  {
930  // A 2D finite element living in either 2D or 3D space.
931  // The edges here are the sides of the element, so the
932  // (misnamed) compute_face_map function does what we want
933  this->compute_face_map(dim, qw, edge);
934  return;
935  }
936 
937  libmesh_assert_equal_to (dim, 3); // 1D is unnecessary and currently unsupported
938 
939  LOG_SCOPE("compute_edge_map()", "FEMap");
940 
941  // We're calculating now!
942  this->determine_calculations();
943 
944  // The number of quadrature points.
945  const unsigned int n_qp = cast_int<unsigned int>(qw.size());
946 
947  // Resize the vectors to hold data at the quadrature points
948  if (calculate_xyz)
949  this->xyz.resize(n_qp);
950  if (calculate_dxyz)
951  {
952  this->dxyzdxi_map.resize(n_qp);
953  this->dxyzdeta_map.resize(n_qp);
954  this->tangents.resize(n_qp);
955  this->normals.resize(n_qp);
956  this->JxW.resize(n_qp);
957  }
958  if (calculate_d2xyz)
959  {
960  this->d2xyzdxi2_map.resize(n_qp);
961  this->d2xyzdxideta_map.resize(n_qp);
962  this->d2xyzdeta2_map.resize(n_qp);
963  this->curvatures.resize(n_qp);
964  }
965 
966  // Clear the entities that will be summed
967  for (unsigned int p=0; p<n_qp; p++)
968  {
969  if (calculate_xyz)
970  this->xyz[p].zero();
971  if (calculate_dxyz)
972  {
973  this->tangents[p].resize(1);
974  this->dxyzdxi_map[p].zero();
975  this->dxyzdeta_map[p].zero();
976  }
977  if (calculate_d2xyz)
978  {
979  this->d2xyzdxi2_map[p].zero();
980  this->d2xyzdxideta_map[p].zero();
981  this->d2xyzdeta2_map[p].zero();
982  }
983  }
984 
985  // compute x, dxdxi at the quadrature points
986  for (unsigned int i=0,
987  psi_map_size=cast_int<unsigned int>(psi_map.size());
988  i != psi_map_size; i++) // sum over the nodes
989  {
990  const Point & edge_point = edge->point(i);
991 
992  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
993  {
994  if (calculate_xyz)
995  this->xyz[p].add_scaled (edge_point, this->psi_map[i][p]);
996  if (calculate_dxyz)
997  this->dxyzdxi_map[p].add_scaled (edge_point, this->dpsidxi_map[i][p]);
998  if (calculate_d2xyz)
999  this->d2xyzdxi2_map[p].add_scaled (edge_point, this->d2psidxi2_map[i][p]);
1000  }
1001  }
1002 
1003  // Compute the tangents at the quadrature point
1004  // FIXME: normals (plural!) and curvatures are uncalculated
1005  if (calculate_dxyz)
1006  for (unsigned int p=0; p<n_qp; p++)
1007  {
1008  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
1009  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
1010 
1011  // compute the jacobian at the quadrature points
1012  const Real the_jac = std::sqrt(this->dxdxi_map(p)*this->dxdxi_map(p) +
1013  this->dydxi_map(p)*this->dydxi_map(p) +
1014  this->dzdxi_map(p)*this->dzdxi_map(p));
1015 
1016  libmesh_assert_greater (the_jac, 0.);
1017 
1018  this->JxW[p] = the_jac*qw[p];
1019  }
1020 }
1021 
1022 
1023 // Explicit FEMap Instantiations
1024 FACE_EDGE_SHAPE_ERROR(0,init_face_shape_functions)
1025 template void FEMap::init_face_shape_functions<1>(const std::vector<Point> &, const Elem *);
1026 template void FEMap::init_face_shape_functions<2>(const std::vector<Point> &, const Elem *);
1027 template void FEMap::init_face_shape_functions<3>(const std::vector<Point> &, const Elem *);
1028 
1029 FACE_EDGE_SHAPE_ERROR(0,init_edge_shape_functions)
1030 template void FEMap::init_edge_shape_functions<1>(const std::vector<Point> &, const Elem *);
1031 template void FEMap::init_edge_shape_functions<2>(const std::vector<Point> &, const Elem *);
1032 template void FEMap::init_edge_shape_functions<3>(const std::vector<Point> &, const Elem *);
1033 
1034 //--------------------------------------------------------------
1035 // Explicit FE instantiations
1036 template void FE<1,LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1037 template void FE<1,LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1038 template void FE<1,LAGRANGE_VEC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1039 template void FE<1,LAGRANGE_VEC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1040 template void FE<1,L2_LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1041 template void FE<1,L2_LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1042 template void FE<1,HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1043 template void FE<1,HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1044 template void FE<1,L2_HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1045 template void FE<1,L2_HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1046 template void FE<1,CLOUGH>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1047 template void FE<1,CLOUGH>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1048 template void FE<1,HERMITE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1049 template void FE<1,HERMITE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1050 template void FE<1,MONOMIAL>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1051 template void FE<1,MONOMIAL>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1052 template void FE<1,SCALAR>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1053 template void FE<1,SCALAR>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1054 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1055 template void FE<1,BERNSTEIN>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1056 template void FE<1,BERNSTEIN>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1057 template void FE<1,SZABAB>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1058 template void FE<1,SZABAB>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1059 #endif
1060 template void FE<1,XYZ>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1061 template void FE<1,XYZ>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1062 
1063 template void FE<2,LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1064 template void FE<2,LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1065 template void FE<2,LAGRANGE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1066 template void FE<2,LAGRANGE_VEC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1067 template void FE<2,LAGRANGE_VEC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1068 template void FE<2,LAGRANGE_VEC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1069 template void FE<2,L2_LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1070 template void FE<2,L2_LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1071 template void FE<2,L2_LAGRANGE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1072 template void FE<2,HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1073 template void FE<2,HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1074 template void FE<2,HIERARCHIC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1075 template void FE<2,L2_HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1076 template void FE<2,L2_HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1077 template void FE<2,L2_HIERARCHIC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1078 template void FE<2,CLOUGH>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1079 template void FE<2,CLOUGH>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1080 template void FE<2,CLOUGH>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1081 template void FE<2,HERMITE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1082 template void FE<2,HERMITE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1083 template void FE<2,HERMITE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1084 template void FE<2,MONOMIAL>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1085 template void FE<2,MONOMIAL>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1086 template void FE<2,MONOMIAL>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1087 template void FE<2,SCALAR>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1088 template void FE<2,SCALAR>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1089 template void FE<2,SCALAR>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1090 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1091 template void FE<2,BERNSTEIN>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1092 template void FE<2,BERNSTEIN>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1093 template void FE<2,BERNSTEIN>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1094 template void FE<2,SZABAB>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1095 template void FE<2,SZABAB>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1096 template void FE<2,SZABAB>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1097 #endif
1098 template void FE<2,SUBDIVISION>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1099 template void FE<2,XYZ>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1100 template void FE<2,XYZ>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1101 template void FE<2,XYZ>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1102 template void FE<2,NEDELEC_ONE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1103 template void FE<2,NEDELEC_ONE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1104 template void FE<2,NEDELEC_ONE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1105 
1106 // Intel 9.1 complained it needed this in devel mode.
1107 //template void FE<2,XYZ>::init_face_shape_functions(const std::vector<Point> &, const Elem *);
1108 
1109 template void FE<3,LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1110 template void FE<3,LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1111 template void FE<3,LAGRANGE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1112 template void FE<3,LAGRANGE_VEC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1113 template void FE<3,LAGRANGE_VEC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1114 template void FE<3,LAGRANGE_VEC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1115 template void FE<3,L2_LAGRANGE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1116 template void FE<3,L2_LAGRANGE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1117 template void FE<3,L2_LAGRANGE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1118 template void FE<3,HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1119 template void FE<3,HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1120 template void FE<3,HIERARCHIC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1121 template void FE<3,L2_HIERARCHIC>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1122 template void FE<3,L2_HIERARCHIC>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1123 template void FE<3,L2_HIERARCHIC>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1124 template void FE<3,CLOUGH>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1125 template void FE<3,CLOUGH>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1126 template void FE<3,CLOUGH>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1127 template void FE<3,HERMITE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1128 template void FE<3,HERMITE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1129 template void FE<3,HERMITE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1130 template void FE<3,MONOMIAL>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1131 template void FE<3,MONOMIAL>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1132 template void FE<3,MONOMIAL>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1133 template void FE<3,SCALAR>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1134 template void FE<3,SCALAR>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1135 template void FE<3,SCALAR>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1136 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1137 template void FE<3,BERNSTEIN>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1138 template void FE<3,BERNSTEIN>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1139 template void FE<3,BERNSTEIN>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1140 template void FE<3,SZABAB>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1141 template void FE<3,SZABAB>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1142 template void FE<3,SZABAB>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1143 #endif
1144 template void FE<3,XYZ>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1145 template void FE<3,XYZ>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1146 template void FE<3,XYZ>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1147 template void FE<3,NEDELEC_ONE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1148 template void FE<3,NEDELEC_ONE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1149 template void FE<3,NEDELEC_ONE>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1150 
1151 // Intel 9.1 complained it needed this in devel mode.
1152 //template void FE<3,XYZ>::init_face_shape_functions(const std::vector<Point> &, const Elem *);
1153 
1154 } // namespace libMesh
libMesh::FEMap::init_edge_shape_functions
void init_edge_shape_functions(const std::vector< Point > &qp, const Elem *edge)
Definition: fe_boundary.C:514
libMesh::FEMap::dzdeta_map
Real dzdeta_map(const unsigned int p) const
Definition: fe_map.h:589
libMesh::FEMap::dpsideta_map
std::vector< std::vector< Real > > dpsideta_map
Definition: fe_map.h:823
libMesh::FEMap::psi_map
std::vector< std::vector< Real > > psi_map
Definition: fe_map.h:811
libMesh::FEMap::d2psideta2_map
std::vector< std::vector< Real > > d2psideta2_map
Definition: fe_map.h:844
libMesh::FE::side_map
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) override
Definition: fe_boundary.C:353
libMesh::FE::shape
static OutputShape shape(const ElemType t, const Order o, const unsigned int i, const Point &p)
libMesh::FEMap::calculate_dxyz
bool calculate_dxyz
Definition: fe_map.h:887
libMesh::FE::map_xi
static Point map_xi(const Elem *elem, const Point &reference_point)
Definition: fe_map.C:2012
libMesh::FEMap::compute_edge_map
void compute_edge_map(int dim, const std::vector< Real > &qw, const Elem *side)
Definition: fe_boundary.C:922
libMesh::SIDEMAP_ERROR
SIDEMAP_ERROR(0, CLOUGH, side_map)
Definition: fe_boundary.C:67
libMesh::REINIT_ERROR
REINIT_ERROR(0, CLOUGH, reinit)
Definition: fe_boundary.C:65
libMesh::FEMap::dxdeta_map
Real dxdeta_map(const unsigned int p) const
Definition: fe_map.h:573
side
unsigned short int side
Definition: xdr_io.C:50
libMesh::FEMap::dzdxi_map
Real dzdxi_map(const unsigned int p) const
Definition: fe_map.h:565
libMesh::Elem
The base class for all geometric element types.
Definition: elem.h:100
libMesh::index_range
IntRange< std::size_t > index_range(const std::vector< T > &vec)
Definition: int_range.h:104
libMesh::FE::shape_deriv
static OutputShape shape_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::Elem::p_level
unsigned int p_level() const
Definition: elem.h:2555
libMesh::zero
const Number zero
Definition: libmesh.h:239
std::max
long double max(long double a, double b)
Definition: libmesh_augment_std_namespace.h:27
libMesh::FEMap::d2psidxideta_map
std::vector< std::vector< Real > > d2psidxideta_map
Definition: fe_map.h:837
tensor_value.h
libMesh::FEMap::d2xyzdxideta_map
std::vector< RealGradient > d2xyzdxideta_map
Definition: fe_map.h:648
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
libMesh::FE::reinit
virtual void reinit(const Elem *elem, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr) override
Definition: fe.C:129
libMesh::FEMap::dxyzdxi_map
std::vector< RealGradient > dxyzdxi_map
Definition: fe_map.h:624
libMesh::FE::n_shape_functions
virtual unsigned int n_shape_functions() const override
Definition: fe.C:50
libMesh::TypeVector::unit
TypeVector< T > unit() const
Definition: type_vector.C:37
libMesh::FEMap::d2xyzdeta2_map
std::vector< RealGradient > d2xyzdeta2_map
Definition: fe_map.h:654
libmesh_common.h
libMesh::Elem::n_nodes
virtual unsigned int n_nodes() const =0
libMesh::FEMap::d2xyzdxi2_map
std::vector< RealGradient > d2xyzdxi2_map
Definition: fe_map.h:642
libMesh::FEMap::calculate_d2xyz
bool calculate_d2xyz
Definition: fe_map.h:892
libmesh_logging.h
libMesh::FEMap::calculate_xyz
bool calculate_xyz
Definition: fe_map.h:882
libMesh::TypeVector::cross
TypeVector< typename CompareTypes< T, T2 >::supertype > cross(const TypeVector< T2 > &v) const
Definition: type_vector.h:882
libMesh::FEMap::curvatures
std::vector< Real > curvatures
Definition: fe_map.h:861
libMesh::FEMap::d2psidxi2_map
std::vector< std::vector< Real > > d2psidxi2_map
Definition: fe_map.h:830
libMesh::FEMap::JxW
std::vector< Real > JxW
Definition: fe_map.h:871
libMesh::FEMap::normals
std::vector< Point > normals
Definition: fe_map.h:854
libMesh::FEMap::xyz
std::vector< Point > xyz
Definition: fe_map.h:618
libMesh::FE::map_eta
static Point map_eta(const Elem *elem, const Point &reference_point)
Definition: fe_map.C:2040
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2050
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:131
libMesh::FEMap::dydeta_map
Real dydeta_map(const unsigned int p) const
Definition: fe_map.h:581
libMesh::FEMap::dpsidxi_map
std::vector< std::vector< Real > > dpsidxi_map
Definition: fe_map.h:817
libMesh::FEMap::tangents
std::vector< std::vector< Point > > tangents
Definition: fe_map.h:849
libMesh::FEMap::dxyzdeta_map
std::vector< RealGradient > dxyzdeta_map
Definition: fe_map.h:630
libMesh::FEMap::dxdxi_map
Real dxdxi_map(const unsigned int p) const
Definition: fe_map.h:549
libMesh::FEMap::determine_calculations
void determine_calculations()
Definition: fe_map.h:527
libMesh::FEMap::dydxi_map
Real dydxi_map(const unsigned int p) const
Definition: fe_map.h:557
libMesh::FEMap::compute_face_map
virtual void compute_face_map(int dim, const std::vector< Real > &qw, const Elem *side)
Definition: fe_boundary.C:574
libMesh::Elem::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i)=0
libMesh::Elem::default_order
virtual Order default_order() const =0
libMesh::FE::edge_reinit
virtual void edge_reinit(const Elem *elem, const unsigned int edge, const Real tolerance=TOLERANCE, const std::vector< Point > *const pts=nullptr, const std::vector< Real > *const weights=nullptr) override
Definition: fe_boundary.C:260
libMesh::Elem::type
virtual ElemType type() const =0
libMesh::Point
A geometric point in (x,y,z) space.
Definition: point.h:38
libMesh::Elem::node_id
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1914
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition: elem.h:1892
libMesh::FE::shape_second_deriv
static OutputShape shape_second_deriv(const ElemType t, const Order o, const unsigned int i, const unsigned int j, const Point &p)
libMesh::FEMap::init_face_shape_functions
void init_face_shape_functions(const std::vector< Point > &qp, const Elem *side)
Definition: fe_boundary.C:409
libMesh::FE::inverse_map
static Point inverse_map(const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: fe_map.C:1576
generated by doxygen