fe_boundary.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
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  libMesh::err << "ERROR: This method makes no sense for low-D elements!" \
48  << std::endl; \
49  libmesh_error(); \
50 }
51 
52 #define SIDEMAP_ERROR(_dim, _type, _func) \
53 template <> \
54 void FE<_dim,_type>::_func(const Elem*, \
55  const Elem*, \
56  const unsigned int, \
57  const std::vector<Point>&, \
58  std::vector<Point>&) \
59 { \
60  libMesh::err << "ERROR: This method makes no sense for low-D elements!" \
61  << std::endl; \
62  libmesh_error(); \
63 }
64 
65 #define FACE_EDGE_SHAPE_ERROR(_dim, _func) \
66 template <> \
67  void FEMap::_func<_dim>(const std::vector<Point>&, \
68  const Elem* ) \
69 { \
70  libMesh::err << "ERROR: This method makes no sense for low-D elements!" \
71  << std::endl; \
72  libmesh_error(); \
73 }
74 
75 
76 REINIT_ERROR(0, CLOUGH, reinit)
77 REINIT_ERROR(0, CLOUGH, edge_reinit)
78 SIDEMAP_ERROR(0, CLOUGH, side_map)
79 REINIT_ERROR(0, HERMITE, reinit)
80 REINIT_ERROR(0, HERMITE, edge_reinit)
81 SIDEMAP_ERROR(0, HERMITE, side_map)
82 REINIT_ERROR(0, HIERARCHIC, reinit)
83 REINIT_ERROR(0, HIERARCHIC, edge_reinit)
84 SIDEMAP_ERROR(0, HIERARCHIC, side_map)
85 REINIT_ERROR(0, L2_HIERARCHIC, reinit)
86 REINIT_ERROR(0, L2_HIERARCHIC, edge_reinit)
87 SIDEMAP_ERROR(0, L2_HIERARCHIC, side_map)
88 REINIT_ERROR(0, LAGRANGE, reinit)
89 REINIT_ERROR(0, LAGRANGE, edge_reinit)
90 SIDEMAP_ERROR(0, LAGRANGE, side_map)
91 REINIT_ERROR(0, LAGRANGE_VEC, reinit)
92 REINIT_ERROR(0, LAGRANGE_VEC, edge_reinit)
93 SIDEMAP_ERROR(0, LAGRANGE_VEC, side_map)
94 REINIT_ERROR(0, L2_LAGRANGE, reinit)
95 REINIT_ERROR(0, L2_LAGRANGE, edge_reinit)
96 SIDEMAP_ERROR(0, L2_LAGRANGE, side_map)
97 REINIT_ERROR(0, MONOMIAL, reinit)
98 REINIT_ERROR(0, MONOMIAL, edge_reinit)
99 SIDEMAP_ERROR(0, MONOMIAL, side_map)
100 REINIT_ERROR(0, SCALAR, reinit)
101 REINIT_ERROR(0, SCALAR, edge_reinit)
102 SIDEMAP_ERROR(0, SCALAR, side_map)
103 REINIT_ERROR(0, XYZ, reinit)
104 REINIT_ERROR(0, XYZ, edge_reinit)
105 SIDEMAP_ERROR(0, XYZ, side_map)
106 REINIT_ERROR(0, NEDELEC_ONE, reinit)
107 REINIT_ERROR(0, NEDELEC_ONE, edge_reinit)
108 SIDEMAP_ERROR(0, NEDELEC_ONE, side_map)
109 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
110 REINIT_ERROR(0, BERNSTEIN, reinit)
111 REINIT_ERROR(0, BERNSTEIN, edge_reinit)
112 SIDEMAP_ERROR(0, BERNSTEIN, side_map)
113 REINIT_ERROR(0, SZABAB, reinit)
114 REINIT_ERROR(0, SZABAB, edge_reinit)
115 SIDEMAP_ERROR(0, SZABAB, side_map)
116 #endif
117 
118 REINIT_ERROR(1, CLOUGH, edge_reinit)
119 REINIT_ERROR(1, HERMITE, edge_reinit)
120 REINIT_ERROR(1, HIERARCHIC, edge_reinit)
121 REINIT_ERROR(1, L2_HIERARCHIC, edge_reinit)
122 REINIT_ERROR(1, LAGRANGE, edge_reinit)
123 REINIT_ERROR(1, LAGRANGE_VEC, edge_reinit)
124 REINIT_ERROR(1, L2_LAGRANGE, edge_reinit)
125 REINIT_ERROR(1, XYZ, edge_reinit)
126 REINIT_ERROR(1, MONOMIAL, edge_reinit)
127 REINIT_ERROR(1, SCALAR, edge_reinit)
128 REINIT_ERROR(1, NEDELEC_ONE, reinit)
129 REINIT_ERROR(1, NEDELEC_ONE, edge_reinit)
130 SIDEMAP_ERROR(1, NEDELEC_ONE, side_map)
131 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
132 REINIT_ERROR(1, BERNSTEIN, edge_reinit)
133 REINIT_ERROR(1, SZABAB, edge_reinit)
134 #endif
135 
136 
137 //-------------------------------------------------------
138 // Methods for 2D, 3D
139 template <unsigned int Dim, FEFamily T>
140 void FE<Dim,T>::reinit(const Elem* elem,
141  const unsigned int s,
142  const Real /* tolerance */,
143  const std::vector<Point>* const pts,
144  const std::vector<Real>* const weights)
145 {
146  libmesh_assert(elem);
147  libmesh_assert (this->qrule != NULL || pts != NULL);
148  // We now do this for 1D elements!
149  // libmesh_assert_not_equal_to (Dim, 1);
150 
151  // Build the side of interest
152  const AutoPtr<Elem> side(elem->build_side(s));
153 
154  // Find the max p_level to select
155  // the right quadrature rule for side integration
156  unsigned int side_p_level = elem->p_level();
157  if (elem->neighbor(s) != NULL)
158  side_p_level = std::max(side_p_level, elem->neighbor(s)->p_level());
159 
160  // Initialize the shape functions at the user-specified
161  // points
162  if (pts != NULL)
163  {
164  // The shape functions do not correspond to the qrule
165  this->shapes_on_quadrature = false;
166 
167  // Initialize the face shape functions
168  this->_fe_map->template init_face_shape_functions<Dim>(*pts, side.get());
169 
170  // Compute the Jacobian*Weight on the face for integration
171  if (weights != NULL)
172  {
173  this->_fe_map->compute_face_map (Dim, *weights, side.get());
174  }
175  else
176  {
177  std::vector<Real> dummy_weights (pts->size(), 1.);
178  this->_fe_map->compute_face_map (Dim, dummy_weights, side.get());
179  }
180  }
181  // If there are no user specified points, we use the
182  // quadrature rule
183  else
184  {
185  // initialize quadrature rule
186  this->qrule->init(side->type(), side_p_level);
187 
188  if(this->qrule->shapes_need_reinit())
189  this->shapes_on_quadrature = false;
190 
191  // FIXME - could this break if the same FE object was used
192  // for both volume and face integrals? - RHS
193  // We might not need to reinitialize the shape functions
194  if ((this->get_type() != elem->type()) ||
195  (side->type() != last_side) ||
196  (this->get_p_level() != side_p_level) ||
197  this->shapes_need_reinit() ||
198  !this->shapes_on_quadrature)
199  {
200  // Set the element type and p_level
201  this->elem_type = elem->type();
202 
203  // Set the last_side
204  last_side = side->type();
205 
206  // Set the last p level
207  this->_p_level = side_p_level;
208 
209  // Initialize the face shape functions
210  this->_fe_map->template init_face_shape_functions<Dim>(this->qrule->get_points(), side.get());
211  }
212 
213  // Compute the Jacobian*Weight on the face for integration
214  this->_fe_map->compute_face_map (Dim, this->qrule->get_weights(), side.get());
215 
216  // The shape functions correspond to the qrule
217  this->shapes_on_quadrature = true;
218  }
219 
220  // make a copy of the Jacobian for integration
221  const std::vector<Real> JxW_int(this->_fe_map->get_JxW());
222 
223  // make a copy of shape on quadrature info
224  bool shapes_on_quadrature_side = this->shapes_on_quadrature;
225 
226  // Find where the integration points are located on the
227  // full element.
228  const std::vector<Point>* ref_qp;
229  if (pts != NULL)
230  ref_qp = pts;
231  else
232  ref_qp = &this->qrule->get_points();
233 
234  std::vector<Point> qp;
235  this->side_map(elem, side.get(), s, *ref_qp, qp);
236 
237  // compute the shape function and derivative values
238  // at the points qp
239  this->reinit (elem, &qp);
240 
241  this->shapes_on_quadrature = shapes_on_quadrature_side;
242 
243  // copy back old data
244  this->_fe_map->get_JxW() = JxW_int;
245 }
246 
247 
248 
249 template <unsigned int Dim, FEFamily T>
250 void FE<Dim,T>::edge_reinit(const Elem* elem,
251  const unsigned int e,
252  const Real tolerance,
253  const std::vector<Point>* const pts,
254  const std::vector<Real>* const weights)
255 {
256  libmesh_assert(elem);
257  libmesh_assert (this->qrule != NULL || pts != NULL);
258  // We don't do this for 1D elements!
259  libmesh_assert_not_equal_to (Dim, 1);
260 
261  // Build the side of interest
262  const AutoPtr<Elem> edge(elem->build_edge(e));
263 
264  // Initialize the shape functions at the user-specified
265  // points
266  if (pts != NULL)
267  {
268  // The shape functions do not correspond to the qrule
269  this->shapes_on_quadrature = false;
270 
271  // Initialize the edge shape functions
272  this->_fe_map->template init_edge_shape_functions<Dim> (*pts, edge.get());
273 
274  // Compute the Jacobian*Weight on the face for integration
275  if (weights != NULL)
276  {
277  this->_fe_map->compute_edge_map (Dim, *weights, edge.get());
278  }
279  else
280  {
281  std::vector<Real> dummy_weights (pts->size(), 1.);
282  this->_fe_map->compute_edge_map (Dim, dummy_weights, edge.get());
283  }
284  }
285  // If there are no user specified points, we use the
286  // quadrature rule
287  else
288  {
289  // initialize quadrature rule
290  this->qrule->init(edge->type(), elem->p_level());
291 
292  if(this->qrule->shapes_need_reinit())
293  this->shapes_on_quadrature = false;
294 
295  // We might not need to reinitialize the shape functions
296  if ((this->get_type() != elem->type()) ||
297  (edge->type() != static_cast<int>(last_edge)) || // Comparison between enum and unsigned, cast the unsigned to int
298  this->shapes_need_reinit() ||
299  !this->shapes_on_quadrature)
300  {
301  // Set the element type
302  this->elem_type = elem->type();
303 
304  // Set the last_edge
305  last_edge = edge->type();
306 
307  // Initialize the edge shape functions
308  this->_fe_map->template init_edge_shape_functions<Dim> (this->qrule->get_points(), edge.get());
309  }
310 
311  // Compute the Jacobian*Weight on the face for integration
312  this->_fe_map->compute_edge_map (Dim, this->qrule->get_weights(), edge.get());
313 
314  // The shape functions correspond to the qrule
315  this->shapes_on_quadrature = true;
316  }
317 
318  // make a copy of the Jacobian for integration
319  const std::vector<Real> JxW_int(this->_fe_map->get_JxW());
320 
321  // Find where the integration points are located on the
322  // full element.
323  std::vector<Point> qp;
324  this->inverse_map (elem, this->_fe_map->get_xyz(), qp, tolerance);
325 
326  // compute the shape function and derivative values
327  // at the points qp
328  this->reinit (elem, &qp);
329 
330  // copy back old data
331  this->_fe_map->get_JxW() = JxW_int;
332 }
333 
334 template <unsigned int Dim, FEFamily T>
335 void FE<Dim,T>::side_map (const Elem* elem,
336  const Elem* side,
337  const unsigned int s,
338  const std::vector<Point>& reference_side_points,
339  std::vector<Point>& reference_points)
340 {
341  unsigned int side_p_level = elem->p_level();
342  if (elem->neighbor(s) != NULL)
343  side_p_level = std::max(side_p_level, elem->neighbor(s)->p_level());
344 
345  if (side->type() != last_side ||
346  side_p_level != this->_p_level ||
347  !this->shapes_on_quadrature)
348  {
349  // Set the element type
350  this->elem_type = elem->type();
351  this->_p_level = side_p_level;
352 
353  // Set the last_side
354  last_side = side->type();
355 
356  // Initialize the face shape functions
357  this->_fe_map->template init_face_shape_functions<Dim>(reference_side_points, side);
358  }
359 
360  const unsigned int n_points =
361  libmesh_cast_int<unsigned int>(reference_side_points.size());
362  reference_points.resize(n_points);
363  for (unsigned int i = 0; i < n_points; i++)
364  reference_points[i].zero();
365 
366  std::vector<unsigned int> elem_nodes_map;
367  elem_nodes_map.resize(side->n_nodes());
368  for (unsigned int j = 0; j < side->n_nodes(); j++)
369  for (unsigned int i = 0; i < elem->n_nodes(); i++)
370  if (side->node(j) == elem->node(i))
371  elem_nodes_map[j] = i;
372  std::vector<Point> refspace_nodes;
373  this->get_refspace_nodes(elem->type(), refspace_nodes);
374 
375  const std::vector<std::vector<Real> >& psi_map = this->_fe_map->get_psi();
376 
377  for (unsigned int i=0; i<psi_map.size(); i++) // sum over the nodes
378  {
379  const Point& side_node = refspace_nodes[elem_nodes_map[i]];
380  for (unsigned int p=0; p<n_points; p++)
381  reference_points[p].add_scaled (side_node, psi_map[i][p]);
382  }
383 }
384 
385 template<unsigned int Dim>
386 void FEMap::init_face_shape_functions(const std::vector<Point>& qp,
387  const Elem* side)
388 {
389  libmesh_assert(side);
390 
394  START_LOG("init_face_shape_functions()", "FEMap");
395 
396  // The element type and order to use in
397  // the map
398  const Order mapping_order (side->default_order());
399  const ElemType mapping_elem_type (side->type());
400 
401  // The number of quadrature points.
402  const unsigned int n_qp = libmesh_cast_int<unsigned int>(qp.size());
403 
404  const unsigned int n_mapping_shape_functions =
405  FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
406  mapping_order);
407 
408  // resize the vectors to hold current data
409  // Psi are the shape functions used for the FE mapping
410  this->psi_map.resize (n_mapping_shape_functions);
411 
412  if (Dim > 1)
413  {
414  this->dpsidxi_map.resize (n_mapping_shape_functions);
415  this->d2psidxi2_map.resize (n_mapping_shape_functions);
416  }
417 
418  if (Dim == 3)
419  {
420  this->dpsideta_map.resize (n_mapping_shape_functions);
421  this->d2psidxideta_map.resize (n_mapping_shape_functions);
422  this->d2psideta2_map.resize (n_mapping_shape_functions);
423  }
424 
425  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
426  {
427  // Allocate space to store the values of the shape functions
428  // and their first and second derivatives at the quadrature points.
429  this->psi_map[i].resize (n_qp);
430  if (Dim > 1)
431  {
432  this->dpsidxi_map[i].resize (n_qp);
433  this->d2psidxi2_map[i].resize (n_qp);
434  }
435  if (Dim == 3)
436  {
437  this->dpsideta_map[i].resize (n_qp);
438  this->d2psidxideta_map[i].resize (n_qp);
439  this->d2psideta2_map[i].resize (n_qp);
440  }
441 
442  // Compute the value of shape function i, and its first and
443  // second derivatives at quadrature point p
444  // (Lagrange shape functions are used for the mapping)
445  for (unsigned int p=0; p<n_qp; p++)
446  {
447  this->psi_map[i][p] = FE<Dim-1,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
448  if (Dim > 1)
449  {
450  this->dpsidxi_map[i][p] = FE<Dim-1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
451  this->d2psidxi2_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 0, qp[p]);
452  }
453  // libMesh::out << "this->d2psidxi2_map["<<i<<"][p]=" << d2psidxi2_map[i][p] << std::endl;
454 
455  // If we are in 3D, then our sides are 2D faces.
456  // For the second derivatives, we must also compute the cross
457  // derivative d^2() / dxi deta
458  if (Dim == 3)
459  {
460  this->dpsideta_map[i][p] = FE<Dim-1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
461  this->d2psidxideta_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 1, qp[p]);
462  this->d2psideta2_map[i][p] = FE<Dim-1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 2, qp[p]);
463  }
464  }
465  }
466 
467 
471  STOP_LOG("init_face_shape_functions()", "FEMap");
472 }
473 
474 template<unsigned int Dim>
475 void FEMap::init_edge_shape_functions(const std::vector<Point>& qp,
476  const Elem* edge)
477 {
478  libmesh_assert(edge);
479 
483  START_LOG("init_edge_shape_functions()", "FEMap");
484 
485  // The element type and order to use in
486  // the map
487  const Order mapping_order (edge->default_order());
488  const ElemType mapping_elem_type (edge->type());
489 
490  // The number of quadrature points.
491  const unsigned int n_qp = libmesh_cast_int<unsigned int>(qp.size());
492 
493  const unsigned int n_mapping_shape_functions =
494  FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
495  mapping_order);
496 
497  // resize the vectors to hold current data
498  // Psi are the shape functions used for the FE mapping
499  this->psi_map.resize (n_mapping_shape_functions);
500  this->dpsidxi_map.resize (n_mapping_shape_functions);
501  this->d2psidxi2_map.resize (n_mapping_shape_functions);
502 
503  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
504  {
505  // Allocate space to store the values of the shape functions
506  // and their first and second derivatives at the quadrature points.
507  this->psi_map[i].resize (n_qp);
508  this->dpsidxi_map[i].resize (n_qp);
509  this->d2psidxi2_map[i].resize (n_qp);
510 
511  // Compute the value of shape function i, and its first and
512  // second derivatives at quadrature point p
513  // (Lagrange shape functions are used for the mapping)
514  for (unsigned int p=0; p<n_qp; p++)
515  {
516  this->psi_map[i][p] = FE<1,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
517  this->dpsidxi_map[i][p] = FE<1,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
518  this->d2psidxi2_map[i][p] = FE<1,LAGRANGE>::shape_second_deriv(mapping_elem_type, mapping_order, i, 0, qp[p]);
519  }
520  }
521 
525  STOP_LOG("init_edge_shape_functions()", "FEMap");
526 }
527 
528 
529 
530 void FEMap::compute_face_map(int dim, const std::vector<Real>& qw,
531  const Elem* side)
532 {
533  libmesh_assert(side);
534 
535  START_LOG("compute_face_map()", "FEMap");
536 
537  // The number of quadrature points.
538  const unsigned int n_qp = libmesh_cast_int<unsigned int>(qw.size());
539 
540  switch (dim)
541  {
542  case 1:
543  {
544  // A 1D finite element, currently assumed to be in 1D space
545  // This means the boundary is a "0D finite element", a
546  // NODEELEM.
547 
548  // Resize the vectors to hold data at the quadrature points
549  {
550  this->xyz.resize(n_qp);
551  normals.resize(n_qp);
552 
553  this->JxW.resize(n_qp);
554  }
555 
556  // If we have no quadrature points, there's nothing else to do
557  if (!n_qp)
558  break;
559 
560  // We need to look back at the full edge to figure out the normal
561  // vector
562  const Elem *elem = side->parent();
563  libmesh_assert (elem);
564  if (side->node(0) == elem->node(0))
565  normals[0] = Point(-1.);
566  else
567  {
568  libmesh_assert_equal_to (side->node(0), elem->node(1));
569  normals[0] = Point(1.);
570  }
571 
572  // Calculate x at the point
573  libmesh_assert_equal_to (this->psi_map.size(), 1);
574  // In the unlikely event we have multiple quadrature
575  // points, they'll be in the same place
576  for (unsigned int p=0; p<n_qp; p++)
577  {
578  this->xyz[p].zero();
579  this->xyz[p].add_scaled (side->point(0), this->psi_map[0][p]);
580  normals[p] = normals[0];
581  this->JxW[p] = 1.0*qw[p];
582  }
583 
584  // done computing the map
585  break;
586  }
587 
588  case 2:
589  {
590  // A 2D finite element living in either 2D or 3D space.
591  // This means the boundary is a 1D finite element, i.e.
592  // and EDGE2 or EDGE3.
593  // Resize the vectors to hold data at the quadrature points
594  {
595  this->xyz.resize(n_qp);
596  this->dxyzdxi_map.resize(n_qp);
597  this->d2xyzdxi2_map.resize(n_qp);
598  this->tangents.resize(n_qp);
599  this->normals.resize(n_qp);
600  this->curvatures.resize(n_qp);
601 
602  this->JxW.resize(n_qp);
603  }
604 
605  // Clear the entities that will be summed
606  // Compute the tangent & normal at the quadrature point
607  for (unsigned int p=0; p<n_qp; p++)
608  {
609  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
610  this->xyz[p].zero();
611  this->dxyzdxi_map[p].zero();
612  this->d2xyzdxi2_map[p].zero();
613  }
614 
615  // compute x, dxdxi at the quadrature points
616  for (unsigned int i=0; i<this->psi_map.size(); i++) // sum over the nodes
617  {
618  const Point& side_point = side->point(i);
619 
620  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
621  {
622  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
623  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
624  this->d2xyzdxi2_map[p].add_scaled(side_point, this->d2psidxi2_map[i][p]);
625  }
626  }
627 
628  // Compute the tangent & normal at the quadrature point
629  for (unsigned int p=0; p<n_qp; p++)
630  {
631  // The first tangent comes from just the edge's Jacobian
632  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
633 
634 #if LIBMESH_DIM == 2
635  // For a 2D element living in 2D, the normal is given directly
636  // from the entries in the edge Jacobian.
637  this->normals[p] = (Point(this->dxyzdxi_map[p](1), -this->dxyzdxi_map[p](0), 0.)).unit();
638 
639 #elif LIBMESH_DIM == 3
640  // For a 2D element living in 3D, there is a second tangent.
641  // For the second tangent, we need to refer to the full
642  // element's (not just the edge's) Jacobian.
643  const Elem *elem = side->parent();
644  libmesh_assert(elem);
645 
646  // Inverse map xyz[p] to a reference point on the parent...
647  Point reference_point = FE<2,LAGRANGE>::inverse_map(elem, this->xyz[p]);
648 
649  // Get dxyz/dxi and dxyz/deta from the parent map.
650  Point dx_dxi = FE<2,LAGRANGE>::map_xi (elem, reference_point);
651  Point dx_deta = FE<2,LAGRANGE>::map_eta(elem, reference_point);
652 
653  // The second tangent vector is formed by crossing these vectors.
654  tangents[p][1] = dx_dxi.cross(dx_deta).unit();
655 
656  // Finally, the normal in this case is given by crossing these
657  // two tangents.
658  normals[p] = tangents[p][0].cross(tangents[p][1]).unit();
659 #endif
660 
661 
662  // The curvature is computed via the familiar Frenet formula:
663  // curvature = [d^2(x) / d (xi)^2] dot [normal]
664  // For a reference, see:
665  // F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill, p. 310
666  //
667  // Note: The sign convention here is different from the
668  // 3D case. Concave-upward curves (smiles) have a positive
669  // curvature. Concave-downward curves (frowns) have a
670  // negative curvature. Be sure to take that into account!
671  const Real numerator = this->d2xyzdxi2_map[p] * this->normals[p];
672  const Real denominator = this->dxyzdxi_map[p].size_sq();
673  libmesh_assert_not_equal_to (denominator, 0);
674  curvatures[p] = numerator / denominator;
675  }
676 
677  // compute the jacobian at the quadrature points
678  for (unsigned int p=0; p<n_qp; p++)
679  {
680  const Real the_jac = this->dxyzdxi_map[p].size();
681 
682  libmesh_assert_greater (the_jac, 0.);
683 
684  this->JxW[p] = the_jac*qw[p];
685  }
686 
687  // done computing the map
688  break;
689  }
690 
691 
692 
693  case 3:
694  {
695  // A 3D finite element living in 3D space.
696  // Resize the vectors to hold data at the quadrature points
697  {
698  this->xyz.resize(n_qp);
699  this->dxyzdxi_map.resize(n_qp);
700  this->dxyzdeta_map.resize(n_qp);
701  this->d2xyzdxi2_map.resize(n_qp);
702  this->d2xyzdxideta_map.resize(n_qp);
703  this->d2xyzdeta2_map.resize(n_qp);
704  this->tangents.resize(n_qp);
705  this->normals.resize(n_qp);
706  this->curvatures.resize(n_qp);
707 
708  this->JxW.resize(n_qp);
709  }
710 
711  // Clear the entities that will be summed
712  for (unsigned int p=0; p<n_qp; p++)
713  {
714  this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
715  this->xyz[p].zero();
716  this->dxyzdxi_map[p].zero();
717  this->dxyzdeta_map[p].zero();
718  this->d2xyzdxi2_map[p].zero();
719  this->d2xyzdxideta_map[p].zero();
720  this->d2xyzdeta2_map[p].zero();
721  }
722 
723  // compute x, dxdxi at the quadrature points
724  for (unsigned int i=0; i<this->psi_map.size(); i++) // sum over the nodes
725  {
726  const Point& side_point = side->point(i);
727 
728  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
729  {
730  this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
731  this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
732  this->dxyzdeta_map[p].add_scaled(side_point, this->dpsideta_map[i][p]);
733  this->d2xyzdxi2_map[p].add_scaled (side_point, this->d2psidxi2_map[i][p]);
734  this->d2xyzdxideta_map[p].add_scaled(side_point, this->d2psidxideta_map[i][p]);
735  this->d2xyzdeta2_map[p].add_scaled (side_point, this->d2psideta2_map[i][p]);
736  }
737  }
738 
739  // Compute the tangents, normal, and curvature at the quadrature point
740  for (unsigned int p=0; p<n_qp; p++)
741  {
742  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
743  this->normals[p] = n.unit();
744  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
745  this->tangents[p][1] = n.cross(this->dxyzdxi_map[p]).unit();
746 
747  // Compute curvature using the typical nomenclature
748  // of the first and second fundamental forms.
749  // For reference, see:
750  // 1) http://mathworld.wolfram.com/MeanCurvature.html
751  // (note -- they are using inward normal)
752  // 2) F.S. Merritt, Mathematics Manual, 1962, McGraw-Hill
753  const Real L = -this->d2xyzdxi2_map[p] * this->normals[p];
754  const Real M = -this->d2xyzdxideta_map[p] * this->normals[p];
755  const Real N = -this->d2xyzdeta2_map[p] * this->normals[p];
756  const Real E = this->dxyzdxi_map[p].size_sq();
757  const Real F = this->dxyzdxi_map[p] * this->dxyzdeta_map[p];
758  const Real G = this->dxyzdeta_map[p].size_sq();
759 
760  const Real numerator = E*N -2.*F*M + G*L;
761  const Real denominator = E*G - F*F;
762  libmesh_assert_not_equal_to (denominator, 0.);
763  curvatures[p] = 0.5*numerator/denominator;
764  }
765 
766  // compute the jacobian at the quadrature points, see
767  // http://sp81.msi.umn.edu:999/fluent/fidap/help/theory/thtoc.htm
768  for (unsigned int p=0; p<n_qp; p++)
769  {
770  const Real g11 = (dxdxi_map(p)*dxdxi_map(p) +
771  dydxi_map(p)*dydxi_map(p) +
772  dzdxi_map(p)*dzdxi_map(p));
773 
774  const Real g12 = (dxdxi_map(p)*dxdeta_map(p) +
775  dydxi_map(p)*dydeta_map(p) +
776  dzdxi_map(p)*dzdeta_map(p));
777 
778  const Real g21 = g12;
779 
780  const Real g22 = (dxdeta_map(p)*dxdeta_map(p) +
781  dydeta_map(p)*dydeta_map(p) +
782  dzdeta_map(p)*dzdeta_map(p));
783 
784 
785  const Real the_jac = std::sqrt(g11*g22 - g12*g21);
786 
787  libmesh_assert_greater (the_jac, 0.);
788 
789  this->JxW[p] = the_jac*qw[p];
790  }
791 
792  // done computing the map
793  break;
794  }
795 
796 
797  default:
798  libmesh_error();
799 
800  }
801  STOP_LOG("compute_face_map()", "FEMap");
802 }
803 
804 
805 
806 
808  const std::vector<Real>& qw,
809  const Elem* edge)
810 {
811  libmesh_assert(edge);
812 
813  if (dim == 2)
814  {
815  // A 2D finite element living in either 2D or 3D space.
816  // The edges here are the sides of the element, so the
817  // (misnamed) compute_face_map function does what we want
818  this->compute_face_map(dim, qw, edge);
819  return;
820  }
821 
822  libmesh_assert_equal_to (dim, 3); // 1D is unnecessary and currently unsupported
823 
824  START_LOG("compute_edge_map()", "FEMap");
825 
826  // The number of quadrature points.
827  const unsigned int n_qp = libmesh_cast_int<unsigned int>(qw.size());
828 
829  // Resize the vectors to hold data at the quadrature points
830  this->xyz.resize(n_qp);
831  this->dxyzdxi_map.resize(n_qp);
832  this->dxyzdeta_map.resize(n_qp);
833  this->d2xyzdxi2_map.resize(n_qp);
834  this->d2xyzdxideta_map.resize(n_qp);
835  this->d2xyzdeta2_map.resize(n_qp);
836  this->tangents.resize(n_qp);
837  this->normals.resize(n_qp);
838  this->curvatures.resize(n_qp);
839 
840  this->JxW.resize(n_qp);
841 
842  // Clear the entities that will be summed
843  for (unsigned int p=0; p<n_qp; p++)
844  {
845  this->tangents[p].resize(1);
846  this->xyz[p].zero();
847  this->dxyzdxi_map[p].zero();
848  this->dxyzdeta_map[p].zero();
849  this->d2xyzdxi2_map[p].zero();
850  this->d2xyzdxideta_map[p].zero();
851  this->d2xyzdeta2_map[p].zero();
852  }
853 
854  // compute x, dxdxi at the quadrature points
855  for (unsigned int i=0; i<this->psi_map.size(); i++) // sum over the nodes
856  {
857  const Point& edge_point = edge->point(i);
858 
859  for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
860  {
861  this->xyz[p].add_scaled (edge_point, this->psi_map[i][p]);
862  this->dxyzdxi_map[p].add_scaled (edge_point, this->dpsidxi_map[i][p]);
863  this->d2xyzdxi2_map[p].add_scaled (edge_point, this->d2psidxi2_map[i][p]);
864  }
865  }
866 
867  // Compute the tangents at the quadrature point
868  // FIXME: normals (plural!) and curvatures are uncalculated
869  for (unsigned int p=0; p<n_qp; p++)
870  {
871  const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
872  this->tangents[p][0] = this->dxyzdxi_map[p].unit();
873 
874  // compute the jacobian at the quadrature points
875  const Real the_jac = std::sqrt(this->dxdxi_map(p)*this->dxdxi_map(p) +
876  this->dydxi_map(p)*this->dydxi_map(p) +
877  this->dzdxi_map(p)*this->dzdxi_map(p));
878 
879  libmesh_assert_greater (the_jac, 0.);
880 
881  this->JxW[p] = the_jac*qw[p];
882  }
883 
884  STOP_LOG("compute_edge_map()", "FEMap");
885 }
886 
887 
888 // Explicit FEMap Instantiations
889 FACE_EDGE_SHAPE_ERROR(0,init_face_shape_functions)
890 template void FEMap::init_face_shape_functions<1>(const std::vector<Point>&,const Elem*);
891 template void FEMap::init_face_shape_functions<2>(const std::vector<Point>&,const Elem*);
892 template void FEMap::init_face_shape_functions<3>(const std::vector<Point>&,const Elem*);
893 
894 FACE_EDGE_SHAPE_ERROR(0,init_edge_shape_functions)
895 template void FEMap::init_edge_shape_functions<1>(const std::vector<Point>&, const Elem*);
896 template void FEMap::init_edge_shape_functions<2>(const std::vector<Point>&, const Elem*);
897 template void FEMap::init_edge_shape_functions<3>(const std::vector<Point>&, const Elem*);
898 
899 //--------------------------------------------------------------
900 // Explicit FE instantiations
901 template void FE<1,LAGRANGE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
902 template void FE<1,LAGRANGE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
903 template void FE<1,LAGRANGE_VEC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
904 template void FE<1,LAGRANGE_VEC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
905 template void FE<1,L2_LAGRANGE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
906 template void FE<1,L2_LAGRANGE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
907 template void FE<1,HIERARCHIC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
908 template void FE<1,HIERARCHIC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
909 template void FE<1,L2_HIERARCHIC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
910 template void FE<1,L2_HIERARCHIC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
911 template void FE<1,CLOUGH>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
912 template void FE<1,CLOUGH>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
913 template void FE<1,HERMITE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
914 template void FE<1,HERMITE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
915 template void FE<1,MONOMIAL>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
916 template void FE<1,MONOMIAL>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
917 template void FE<1,SCALAR>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
918 template void FE<1,SCALAR>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
919 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
920 template void FE<1,BERNSTEIN>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
921 template void FE<1,BERNSTEIN>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
922 template void FE<1,SZABAB>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
923 template void FE<1,SZABAB>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
924 #endif
925 template void FE<1,XYZ>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
926 template void FE<1,XYZ>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
927 
928 template void FE<2,LAGRANGE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
929 template void FE<2,LAGRANGE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
930 template void FE<2,LAGRANGE>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
931 template void FE<2,LAGRANGE_VEC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
932 template void FE<2,LAGRANGE_VEC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
933 template void FE<2,LAGRANGE_VEC>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
934 template void FE<2,L2_LAGRANGE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
935 template void FE<2,L2_LAGRANGE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
936 template void FE<2,L2_LAGRANGE>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
937 template void FE<2,HIERARCHIC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
938 template void FE<2,HIERARCHIC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
939 template void FE<2,HIERARCHIC>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
940 template void FE<2,L2_HIERARCHIC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
941 template void FE<2,L2_HIERARCHIC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
942 template void FE<2,L2_HIERARCHIC>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
943 template void FE<2,CLOUGH>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
944 template void FE<2,CLOUGH>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
945 template void FE<2,CLOUGH>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
946 template void FE<2,HERMITE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
947 template void FE<2,HERMITE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
948 template void FE<2,HERMITE>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
949 template void FE<2,MONOMIAL>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
950 template void FE<2,MONOMIAL>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
951 template void FE<2,MONOMIAL>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
952 template void FE<2,SCALAR>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
953 template void FE<2,SCALAR>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
954 template void FE<2,SCALAR>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
955 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
956 template void FE<2,BERNSTEIN>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
957 template void FE<2,BERNSTEIN>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
958 template void FE<2,BERNSTEIN>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
959 template void FE<2,SZABAB>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
960 template void FE<2,SZABAB>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
961 template void FE<2,SZABAB>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
962 #endif
963 template void FE<2,XYZ>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
964 template void FE<2,XYZ>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
965 template void FE<2,XYZ>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
966 template void FE<2,NEDELEC_ONE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
967 template void FE<2,NEDELEC_ONE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
968 template void FE<2,NEDELEC_ONE>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
969 
970 // Intel 9.1 complained it needed this in devel mode.
971 //template void FE<2,XYZ>::init_face_shape_functions(const std::vector<Point>&, const Elem*);
972 
973 template void FE<3,LAGRANGE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
974 template void FE<3,LAGRANGE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
975 template void FE<3,LAGRANGE>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
976 template void FE<3,LAGRANGE_VEC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
977 template void FE<3,LAGRANGE_VEC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
978 template void FE<3,LAGRANGE_VEC>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
979 template void FE<3,L2_LAGRANGE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
980 template void FE<3,L2_LAGRANGE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
981 template void FE<3,L2_LAGRANGE>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
982 template void FE<3,HIERARCHIC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
983 template void FE<3,HIERARCHIC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
984 template void FE<3,HIERARCHIC>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
985 template void FE<3,L2_HIERARCHIC>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
986 template void FE<3,L2_HIERARCHIC>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
987 template void FE<3,L2_HIERARCHIC>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
988 template void FE<3,CLOUGH>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
989 template void FE<3,CLOUGH>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
990 template void FE<3,CLOUGH>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
991 template void FE<3,HERMITE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
992 template void FE<3,HERMITE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
993 template void FE<3,HERMITE>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
994 template void FE<3,MONOMIAL>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
995 template void FE<3,MONOMIAL>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
996 template void FE<3,MONOMIAL>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
997 template void FE<3,SCALAR>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
998 template void FE<3,SCALAR>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
999 template void FE<3,SCALAR>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1000 #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1001 template void FE<3,BERNSTEIN>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1002 template void FE<3,BERNSTEIN>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
1003 template void FE<3,BERNSTEIN>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1004 template void FE<3,SZABAB>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1005 template void FE<3,SZABAB>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
1006 template void FE<3,SZABAB>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1007 #endif
1008 template void FE<3,XYZ>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1009 template void FE<3,XYZ>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
1010 template void FE<3,XYZ>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1011 template void FE<3,NEDELEC_ONE>::reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1012 template void FE<3,NEDELEC_ONE>::side_map(Elem const*, Elem const*, const unsigned int, const std::vector<Point>&, std::vector<Point>&);
1013 template void FE<3,NEDELEC_ONE>::edge_reinit(Elem const*, unsigned int, Real, const std::vector<Point>* const, const std::vector<Real>* const);
1014 
1015 // Intel 9.1 complained it needed this in devel mode.
1016 //template void FE<3,XYZ>::init_face_shape_functions(const std::vector<Point>&, const Elem*);
1017 
1018 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo