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