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