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