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