Line data Source code
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 :
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(L2_LAGRANGE_GLL)
83 0 : LIBMESH_ERRORS_IN_LOW_D(LAGRANGE_VEC)
84 0 : LIBMESH_ERRORS_IN_LOW_D(L2_LAGRANGE_VEC)
85 0 : LIBMESH_ERRORS_IN_LOW_D(MONOMIAL)
86 0 : LIBMESH_ERRORS_IN_LOW_D(MONOMIAL_VEC)
87 0 : LIBMESH_ERRORS_IN_LOW_D(NEDELEC_ONE)
88 0 : LIBMESH_ERRORS_IN_LOW_D(RAVIART_THOMAS)
89 0 : LIBMESH_ERRORS_IN_LOW_D(L2_RAVIART_THOMAS)
90 0 : LIBMESH_ERRORS_IN_LOW_D(SCALAR)
91 0 : LIBMESH_ERRORS_IN_LOW_D(XYZ)
92 :
93 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
94 0 : LIBMESH_ERRORS_IN_LOW_D(BERNSTEIN)
95 0 : LIBMESH_ERRORS_IN_LOW_D(SZABAB)
96 0 : LIBMESH_ERRORS_IN_LOW_D(RATIONAL_BERNSTEIN)
97 : #endif
98 :
99 : // Special error instantiations
100 0 : REINIT_ERROR(1, NEDELEC_ONE, reinit) { libmesh_error_msg("ERROR: Cannot reinit 1D NEDELEC_ONE elements!"); }
101 0 : SIDEMAP_ERROR(1, NEDELEC_ONE, side_map) { libmesh_error_msg("ERROR: Cannot side_map 1D NEDELEC_ONE elements!"); }
102 0 : REINIT_ERROR(1, RAVIART_THOMAS, reinit) { libmesh_error_msg("ERROR: Cannot reinit 1D RAVIART_THOMAS elements!"); }
103 0 : SIDEMAP_ERROR(1, RAVIART_THOMAS, side_map) { libmesh_error_msg("ERROR: Cannot side_map 1D RAVIART_THOMAS elements!"); }
104 0 : REINIT_ERROR(1, L2_RAVIART_THOMAS, reinit) { libmesh_error_msg("ERROR: Cannot reinit 1D L2_RAVIART_THOMAS elements!"); }
105 0 : SIDEMAP_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
109 : template <unsigned int Dim, FEFamily T>
110 9526172 : void 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 2508959 : libmesh_assert(elem);
117 2508959 : 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 9526172 : if (this->calculating_nothing())
125 : {
126 550655 : this->calculations_started = true; // Ironic
127 550655 : 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 8975517 : this->_fe_map->add_calculations();
134 2343246 : this->_fe_map->get_JxW();
135 2343246 : this->_fe_map->get_xyz();
136 8975517 : this->determine_calculations();
137 :
138 : // Build the side of interest
139 8975517 : 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 8975517 : unsigned int side_p_level = elem->p_level();
144 11318763 : if (elem->neighbor_ptr(s) != nullptr)
145 10173943 : 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 8975517 : if (pts != nullptr)
150 : {
151 : // The shape functions do not correspond to the qrule
152 1113 : this->shapes_on_quadrature = false;
153 :
154 : // Initialize the face shape functions
155 1113 : this->_fe_map->template init_face_shape_functions<Dim>(*pts, side.get());
156 :
157 : // Compute the Jacobian*Weight on the face for integration
158 1113 : if (weights != nullptr)
159 : {
160 0 : this->_fe_map->compute_face_map (Dim, *weights, side.get());
161 : }
162 : else
163 : {
164 1749 : std::vector<Real> dummy_weights (pts->size(), 1.);
165 1113 : 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 8974404 : this->qrule->init(*side, side_p_level);
174 :
175 8974404 : if (this->qrule->shapes_need_reinit())
176 0 : 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 15507799 : if ((this->get_type() != elem->type()) ||
182 8836594 : (elem->runtime_topology() &&
183 0 : this->_elem != elem) ||
184 8836594 : (side->type() != last_side) ||
185 13025020 : (this->_elem_p_level != side_p_level) ||
186 15494433 : this->shapes_need_reinit() ||
187 4525146 : !this->shapes_on_quadrature)
188 : {
189 : // Set the element
190 4449315 : this->_elem = elem;
191 4449315 : this->_elem_type = elem->type();
192 :
193 : // Set the last_side
194 4449315 : last_side = side->type();
195 :
196 : // Set the last p level
197 4449315 : this->_p_level = this->_add_p_level_in_reinit * side_p_level;
198 :
199 : // Initialize the face shape functions
200 5593556 : this->_fe_map->template init_face_shape_functions<Dim>(this->qrule->get_points(), side.get());
201 : }
202 : else
203 4525089 : this->_elem = elem;
204 :
205 : // Compute the Jacobian*Weight on the face for integration
206 11317332 : this->_fe_map->compute_face_map (Dim, this->qrule->get_weights(), side.get());
207 :
208 : // The shape functions correspond to the qrule
209 8974404 : this->shapes_on_quadrature = true;
210 : }
211 :
212 : // make a copy of the Jacobian for integration
213 11318763 : const std::vector<Real> JxW_int(this->_fe_map->get_JxW());
214 :
215 : // make a copy of shape on quadrature info
216 8975517 : 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 8975517 : if (pts != nullptr)
222 318 : ref_qp = pts;
223 : else
224 8974404 : ref_qp = &this->qrule->get_points();
225 :
226 4686492 : std::vector<Point> qp;
227 8975517 : 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 8975517 : this->reinit (elem, &qp);
232 :
233 8975517 : this->shapes_on_quadrature = shapes_on_quadrature_side;
234 :
235 8975517 : this->_elem_p_level = side_p_level;
236 :
237 : // copy back old data
238 8975517 : this->_fe_map->get_JxW() = JxW_int;
239 4289025 : }
240 :
241 :
242 :
243 : template <unsigned int Dim, FEFamily T>
244 50934 : void 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 13087 : libmesh_assert(elem);
251 13087 : 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 50934 : this->_fe_map->add_calculations();
259 13087 : this->_fe_map->get_JxW();
260 13087 : this->_fe_map->get_xyz();
261 50934 : this->determine_calculations();
262 :
263 : // Build the side of interest
264 50934 : 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 50934 : if (pts != nullptr)
269 : {
270 : // The shape functions do not correspond to the qrule
271 0 : this->shapes_on_quadrature = false;
272 :
273 : // Initialize the edge shape functions
274 0 : this->_fe_map->template init_edge_shape_functions<Dim> (*pts, edge.get());
275 :
276 : // Compute the Jacobian*Weight on the face for integration
277 0 : if (weights != nullptr)
278 : {
279 0 : this->_fe_map->compute_edge_map (Dim, *weights, edge.get());
280 : }
281 : else
282 : {
283 0 : std::vector<Real> dummy_weights (pts->size(), 1.);
284 0 : 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 64021 : this->qrule->init(*edge, elem->p_level());
293 :
294 50934 : if (this->qrule->shapes_need_reinit())
295 0 : this->shapes_on_quadrature = false;
296 :
297 : // We might not need to reinitialize the shape functions
298 87585 : if ((this->get_type() != elem->type()) ||
299 49280 : (elem->runtime_topology() &&
300 0 : this->_elem != elem) ||
301 73302 : (edge->type() != last_edge) ||
302 87585 : this->shapes_need_reinit() ||
303 10782 : !this->shapes_on_quadrature)
304 : {
305 : // Set the element
306 50934 : this->_elem = elem;
307 50934 : this->_elem_type = elem->type();
308 :
309 : // Set the last_edge
310 50934 : last_edge = edge->type();
311 :
312 : // Initialize the edge shape functions
313 64021 : this->_fe_map->template init_edge_shape_functions<Dim> (this->qrule->get_points(), edge.get());
314 : }
315 : else
316 0 : this->_elem = elem;
317 :
318 : // Compute the Jacobian*Weight on the face for integration
319 64021 : this->_fe_map->compute_edge_map (Dim, this->qrule->get_weights(), edge.get());
320 :
321 : // The shape functions correspond to the qrule
322 50934 : this->shapes_on_quadrature = true;
323 : }
324 :
325 : // make a copy of the Jacobian for integration
326 64021 : 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 50934 : if (pts != nullptr)
332 0 : ref_qp = pts;
333 : else
334 50934 : ref_qp = & this->qrule->get_points();
335 :
336 26174 : std::vector<Point> qp;
337 50934 : 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 50934 : this->reinit (elem, &qp);
342 :
343 : // copy back old data
344 50934 : this->_fe_map->get_JxW() = JxW_int;
345 50934 : }
346 :
347 : template <unsigned int Dim, FEFamily T>
348 8984797 : void 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 8984797 : this->calculate_phi = true;
356 8984797 : this->determine_calculations();
357 :
358 8984797 : unsigned int side_p_level = elem->p_level();
359 11330363 : if (elem->neighbor_ptr(s) != nullptr)
360 10185570 : side_p_level = std::max(side_p_level, elem->neighbor_ptr(s)->p_level());
361 :
362 15623208 : if (side->type() != last_side ||
363 8983649 : (elem->runtime_topology() &&
364 0 : this->_elem != elem) ||
365 20313400 : side_p_level != this->_elem_p_level ||
366 8982554 : !this->shapes_on_quadrature)
367 : {
368 : // Set the element type
369 11461 : this->_elem = elem;
370 11461 : this->_elem_type = elem->type();
371 11461 : this->_elem_p_level = side_p_level;
372 11461 : this->_p_level = this->_add_p_level_in_reinit * side_p_level;
373 :
374 : // Set the last_side
375 11461 : last_side = side->type();
376 :
377 : // Initialize the face shape functions
378 11461 : this->_fe_map->template init_face_shape_functions<Dim>(reference_side_points, side);
379 : }
380 : else
381 8973336 : this->_elem = elem;
382 :
383 : const unsigned int n_points =
384 4691132 : cast_int<unsigned int>(reference_side_points.size());
385 8984797 : reference_points.resize(n_points);
386 37257483 : for (unsigned int i = 0; i < n_points; i++)
387 28272686 : reference_points[i].zero();
388 :
389 4691132 : std::vector<Point> refspace_nodes;
390 8984797 : this->get_refspace_nodes(elem->type(), refspace_nodes);
391 :
392 2345566 : const std::vector<std::vector<Real>> & psi_map = this->_fe_map->get_psi();
393 :
394 : // sum over the nodes
395 50591827 : for (auto i : index_range(psi_map))
396 : {
397 41607030 : const Point & side_node = refspace_nodes[elem->local_side_node(s,i)];
398 199816427 : for (unsigned int p=0; p<n_points; p++)
399 238547269 : reference_points[p].add_scaled (side_node, psi_map[i][p]);
400 : }
401 8984797 : }
402 :
403 : template <unsigned int Dim, FEFamily T>
404 50934 : void 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 50934 : this->calculate_phi = true;
412 50934 : this->determine_calculations();
413 :
414 26174 : unsigned int edge_p_level = elem->p_level();
415 :
416 88781 : if (edge->type() != last_edge ||
417 50934 : (elem->runtime_topology() &&
418 0 : this->_elem != elem) ||
419 114955 : edge_p_level != this->_elem_p_level ||
420 50934 : !this->shapes_on_quadrature)
421 : {
422 : // Set the element type
423 0 : this->_elem = elem;
424 0 : this->_elem_type = elem->type();
425 0 : this->_elem_p_level = edge_p_level;
426 0 : this->_p_level = this->_add_p_level_in_reinit * edge_p_level;
427 :
428 : // Set the last_edge
429 0 : last_edge = edge->type();
430 :
431 : // Initialize the edge shape functions
432 0 : this->_fe_map->template init_edge_shape_functions<Dim>(reference_edge_points, edge);
433 : }
434 : else
435 50934 : this->_elem = elem;
436 :
437 : const unsigned int n_points =
438 26174 : cast_int<unsigned int>(reference_edge_points.size());
439 50934 : reference_points.resize(n_points);
440 245412 : for (unsigned int i = 0; i < n_points; i++)
441 194478 : reference_points[i].zero();
442 :
443 26174 : std::vector<Point> refspace_nodes;
444 50934 : this->get_refspace_nodes(elem->type(), refspace_nodes);
445 :
446 13087 : const std::vector<std::vector<Real>> & psi_map = this->_fe_map->get_psi();
447 :
448 : // sum over the nodes
449 203724 : for (auto i : index_range(psi_map))
450 : {
451 152790 : const Point & edge_node = refspace_nodes[elem->local_edge_node(e,i)];
452 736200 : for (unsigned int p=0; p<n_points; p++)
453 883206 : reference_points[p].add_scaled (edge_node, psi_map[i][p]);
454 : }
455 50934 : }
456 :
457 :
458 : template<unsigned int Dim>
459 5035561 : void FEMap::init_face_shape_functions(const std::vector<Point> & qp,
460 : const Elem * side)
461 : {
462 : // Start logging the shape function initialization
463 3872681 : LOG_SCOPE("init_face_shape_functions()", "FEMap");
464 :
465 1290893 : libmesh_assert(side);
466 :
467 : // We're calculating now!
468 1290893 : this->determine_calculations();
469 :
470 : // The element type and order to use in
471 : // the map
472 5035561 : const FEFamily mapping_family = FEMap::map_fe_type(*side);
473 5035561 : const FEType map_fe_type(side->default_order(), mapping_family);
474 :
475 : // The number of quadrature points.
476 2581788 : 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 5035561 : 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 5035561 : if (calculate_xyz)
485 5026281 : this->psi_map.resize (n_mapping_shape_functions);
486 :
487 : if (Dim > 1)
488 : {
489 5012166 : if (calculate_dxyz)
490 4952214 : this->dpsidxi_map.resize (n_mapping_shape_functions);
491 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
492 5012166 : if (calculate_d2xyz)
493 141075 : this->d2psidxi2_map.resize (n_mapping_shape_functions);
494 : #endif
495 : }
496 :
497 : if (Dim == 3)
498 : {
499 3128133 : if (calculate_dxyz)
500 3078981 : this->dpsideta_map.resize (n_mapping_shape_functions);
501 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
502 3128133 : if (calculate_d2xyz)
503 : {
504 97600 : this->d2psidxideta_map.resize (n_mapping_shape_functions);
505 97600 : this->d2psideta2_map.resize (n_mapping_shape_functions);
506 : }
507 : #endif
508 : }
509 :
510 : FEInterface::shape_ptr shape_ptr =
511 5035561 : FEInterface::shape_function(map_fe_type, side);
512 :
513 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
514 5035561 : 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 5035561 : FEInterface::shape_second_deriv_function(map_fe_type, side);
519 : #endif
520 :
521 32619564 : 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 27584003 : if (calculate_xyz)
526 34530594 : this->psi_map[i].resize (n_qp);
527 : if (Dim > 1)
528 : {
529 27560608 : if (calculate_dxyz)
530 34077156 : this->dpsidxi_map[i].resize (n_qp);
531 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
532 27560608 : if (calculate_d2xyz)
533 1052591 : this->d2psidxi2_map[i].resize (n_qp);
534 : #endif
535 : }
536 : if (Dim == 3)
537 : {
538 21925212 : if (calculate_dxyz)
539 26979462 : this->dpsideta_map[i].resize (n_qp);
540 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
541 21925212 : if (calculate_d2xyz)
542 : {
543 890000 : this->d2psidxideta_map[i].resize (n_qp);
544 890000 : 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 142760172 : for (unsigned int p=0; p<n_qp; p++)
553 : {
554 115176169 : if (calculate_xyz)
555 143877827 : this->psi_map[i][p] = shape_ptr (map_fe_type, side, i, qp[p], false);
556 : if (Dim > 1)
557 : {
558 115152774 : if (calculate_dxyz)
559 142437009 : this->dpsidxi_map[i][p] = shape_deriv_ptr (map_fe_type, side, i, 0, qp[p], false);
560 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
561 115152774 : if (calculate_d2xyz)
562 8397982 : 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 99956073 : if (calculate_dxyz)
573 123273259 : this->dpsideta_map[i][p] = shape_deriv_ptr (map_fe_type, side, i, 1, qp[p], false);
574 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
575 99956073 : if (calculate_d2xyz)
576 : {
577 8042800 : this->d2psidxideta_map[i][p] = shape_second_deriv_ptr(map_fe_type, side, i, 1, qp[p], false);
578 8042800 : 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 5035561 : }
585 :
586 : template<unsigned int Dim>
587 50934 : void FEMap::init_edge_shape_functions(const std::vector<Point> & qp,
588 : const Elem * edge)
589 : {
590 : // Start logging the shape function initialization
591 39261 : LOG_SCOPE("init_edge_shape_functions()", "FEMap");
592 :
593 13087 : libmesh_assert(edge);
594 :
595 : // We're calculating now!
596 13087 : this->determine_calculations();
597 :
598 : // The element type and order to use in
599 : // the map
600 50934 : const FEFamily mapping_family = FEMap::map_fe_type(*edge);
601 50934 : const FEType map_fe_type(edge->default_order(), mapping_family);
602 :
603 : // The number of quadrature points.
604 26174 : 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 50934 : 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 50934 : if (calculate_xyz)
613 50934 : this->psi_map.resize (n_mapping_shape_functions);
614 50934 : if (calculate_dxyz)
615 50934 : this->dpsidxi_map.resize (n_mapping_shape_functions);
616 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
617 50934 : if (calculate_d2xyz)
618 0 : this->d2psidxi2_map.resize (n_mapping_shape_functions);
619 : #endif
620 :
621 : FEInterface::shape_ptr shape_ptr =
622 50934 : FEInterface::shape_function(map_fe_type, edge);
623 :
624 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
625 50934 : 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 50934 : FEInterface::shape_second_deriv_function(map_fe_type, edge);
630 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
631 :
632 203724 : 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 152790 : if (calculate_xyz)
637 192048 : this->psi_map[i].resize (n_qp);
638 152790 : if (calculate_dxyz)
639 192048 : this->dpsidxi_map[i].resize (n_qp);
640 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
641 152790 : if (calculate_d2xyz)
642 0 : 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 736200 : for (unsigned int p=0; p<n_qp; p++)
648 : {
649 583410 : if (calculate_xyz)
650 733308 : this->psi_map[i][p] = shape_ptr (map_fe_type, edge, i, qp[p], false);
651 583410 : if (calculate_dxyz)
652 733308 : this->dpsidxi_map[i][p] = shape_deriv_ptr (map_fe_type, edge, i, 0, qp[p], false);
653 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
654 583410 : if (calculate_d2xyz)
655 0 : this->d2psidxi2_map[i][p] = shape_second_deriv_ptr(map_fe_type, edge, i, 0, qp[p], false);
656 : #endif
657 : }
658 : }
659 50934 : }
660 :
661 :
662 :
663 9549189 : void FEMap::compute_face_map(int dim, const std::vector<Real> & qw,
664 : const Elem * side)
665 : {
666 2486664 : libmesh_assert(side);
667 :
668 : // We're calculating now!
669 2486664 : this->determine_calculations();
670 :
671 7459992 : LOG_ONE_ENTITY_SCOPE("compute_face_map()", "FEMap");
672 :
673 : // The number of quadrature points.
674 4973328 : const unsigned int n_qp = cast_int<unsigned int>(qw.size());
675 :
676 9549189 : const FEFamily mapping_family = FEMap::map_fe_type(*side);
677 9549189 : const Order mapping_order (side->default_order());
678 2486664 : 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 9549189 : FEInterface::n_shape_functions(map_fe_type, /*extra_order=*/0, side);
683 :
684 9549189 : switch (dim)
685 : {
686 24218 : 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 24218 : if (calculate_xyz)
695 24218 : this->xyz.resize(n_qp);
696 24218 : if (calculate_dxyz)
697 24218 : normals.resize(n_qp);
698 :
699 24218 : if (calculate_dxyz)
700 24218 : this->JxW.resize(n_qp);
701 : }
702 :
703 : // If we have no quadrature points, there's nothing else to do
704 24218 : if (!n_qp)
705 0 : break;
706 :
707 : // We need to look back at the full edge to figure out the normal
708 : // vector
709 24218 : const Elem * elem = side->interior_parent();
710 6059 : libmesh_assert (elem);
711 24218 : if (calculate_dxyz)
712 : {
713 36336 : if (side->node_id(0) == elem->node_id(0))
714 : {
715 2185 : const Point reference_point = Point(-1.);
716 8730 : Point dx_dxi = FEMap::map_deriv (1, elem, 0, reference_point);
717 8730 : normals[0] = -dx_dxi.unit();
718 : }
719 : else
720 : {
721 3874 : libmesh_assert_equal_to (side->node_id(0),
722 : elem->node_id(1));
723 3874 : const Point reference_point = Point(1.);
724 15488 : Point dx_dxi = FEMap::map_deriv (1, elem, 0, reference_point);
725 15488 : normals[0] = dx_dxi.unit();
726 : }
727 : }
728 :
729 : // Calculate x at the point
730 24218 : if (calculate_xyz)
731 6059 : 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 48436 : for (unsigned int p=0; p<n_qp; p++)
735 : {
736 24218 : if (calculate_xyz)
737 : {
738 24218 : this->xyz[p].zero();
739 30277 : this->xyz[p].add_scaled (side->point(0), this->psi_map[0][p]);
740 : }
741 24218 : if (calculate_dxyz)
742 : {
743 24218 : normals[p] = normals[0];
744 36336 : this->JxW[p] = 1.0*qw[p];
745 : }
746 : }
747 :
748 : // done computing the map
749 6059 : break;
750 : }
751 :
752 4635017 : 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 4635017 : if (calculate_xyz)
760 4635017 : this->xyz.resize(n_qp);
761 4635017 : if (calculate_dxyz)
762 : {
763 4624217 : this->dxyzdxi_map.resize(n_qp);
764 4624217 : this->tangents.resize(n_qp);
765 4624217 : this->normals.resize(n_qp);
766 :
767 4624217 : this->JxW.resize(n_qp);
768 : }
769 :
770 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
771 4635017 : if (calculate_d2xyz)
772 : {
773 43475 : this->d2xyzdxi2_map.resize(n_qp);
774 43475 : 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 13455324 : for (unsigned int p=0; p<n_qp; p++)
782 : {
783 8820307 : if (calculate_xyz)
784 8820307 : this->xyz[p].zero();
785 8820307 : if (calculate_dxyz)
786 : {
787 11156294 : this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
788 4715174 : this->dxyzdxi_map[p].zero();
789 : }
790 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
791 8820307 : if (calculate_d2xyz)
792 94950 : this->d2xyzdxi2_map[p].zero();
793 : #endif
794 : }
795 :
796 : // compute x, dxdxi at the quadrature points
797 16426808 : for (unsigned int i=0; i<n_mapping_shape_functions; i++) // sum over the nodes
798 : {
799 6345676 : const Point & side_point = side->point(i);
800 :
801 36043484 : for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
802 : {
803 24251693 : if (calculate_xyz)
804 37207359 : this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
805 24251693 : if (calculate_dxyz)
806 37110159 : this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
807 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
808 24251693 : if (calculate_d2xyz)
809 425514 : 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 4635017 : if (calculate_dxyz)
816 : {
817 13422924 : for (unsigned int p=0; p<n_qp; p++)
818 : {
819 : // The first tangent comes from just the edge's Jacobian
820 11156294 : 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 8798707 : const Elem * elem = side->interior_parent();
832 2357587 : libmesh_assert(elem);
833 :
834 : // Inverse map xyz[p] to a reference point on the parent...
835 11156294 : Point reference_point = FEMap::inverse_map(2, elem, this->xyz[p]);
836 :
837 : // Get dxyz/dxi and dxyz/deta from the parent map.
838 8798707 : Point dx_dxi = FEMap::map_deriv (2, elem, 0, reference_point);
839 8798707 : Point dx_deta = FEMap::map_deriv (2, elem, 1, reference_point);
840 :
841 : // The second tangent vector is formed by crossing these vectors.
842 8798707 : 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 8798707 : 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 8798707 : if (calculate_d2xyz)
861 : {
862 46888 : const Real numerator = this->d2xyzdxi2_map[p] * this->normals[p];
863 46888 : const Real denominator = this->dxyzdxi_map[p].norm_sq();
864 23444 : libmesh_assert_not_equal_to (denominator, 0);
865 118394 : curvatures[p] = numerator / denominator;
866 : }
867 : #endif
868 : }
869 :
870 : // compute the jacobian at the quadrature points
871 13422924 : for (unsigned int p=0; p<n_qp; p++)
872 : {
873 8798707 : const Real the_jac = this->dxyzdxi_map[p].norm();
874 :
875 2357587 : libmesh_assert_greater (the_jac, 0.);
876 :
877 13513881 : this->JxW[p] = the_jac*qw[p];
878 : }
879 : }
880 :
881 : // done computing the map
882 1252996 : break;
883 : }
884 :
885 :
886 :
887 4889954 : 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 4889954 : if (calculate_xyz)
893 4889954 : this->xyz.resize(n_qp);
894 4889954 : if (calculate_dxyz)
895 : {
896 4840802 : this->dxyzdxi_map.resize(n_qp);
897 4840802 : this->dxyzdeta_map.resize(n_qp);
898 4840802 : this->tangents.resize(n_qp);
899 4840802 : this->normals.resize(n_qp);
900 4840802 : this->JxW.resize(n_qp);
901 : }
902 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
903 4889954 : if (calculate_d2xyz)
904 : {
905 97600 : this->d2xyzdxi2_map.resize(n_qp);
906 97600 : this->d2xyzdxideta_map.resize(n_qp);
907 97600 : this->d2xyzdeta2_map.resize(n_qp);
908 97600 : this->curvatures.resize(n_qp);
909 : }
910 : #endif
911 : }
912 :
913 : // Clear the entities that will be summed
914 26376691 : for (unsigned int p=0; p<n_qp; p++)
915 : {
916 21486737 : if (calculate_xyz)
917 21486737 : this->xyz[p].zero();
918 21486737 : if (calculate_dxyz)
919 : {
920 26646531 : this->tangents[p].resize(LIBMESH_DIM-1); // 1 Tangent in 2D, 2 in 3D
921 10712804 : this->dxyzdxi_map[p].zero();
922 10712804 : this->dxyzdeta_map[p].zero();
923 : }
924 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
925 21486737 : if (calculate_d2xyz)
926 : {
927 873920 : this->d2xyzdxi2_map[p].zero();
928 436960 : this->d2xyzdxideta_map[p].zero();
929 436960 : this->d2xyzdeta2_map[p].zero();
930 : }
931 : #endif
932 : }
933 :
934 : // compute x, dxdxi at the quadrature points
935 38380359 : for (unsigned int i=0; i<n_mapping_shape_functions; i++) // sum over the nodes
936 : {
937 16798252 : const Point & side_point = side->point(i);
938 :
939 181746947 : for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
940 : {
941 148256542 : if (calculate_xyz)
942 222788158 : this->xyz[p].add_scaled (side_point, this->psi_map[i][p]);
943 148256542 : if (calculate_dxyz)
944 : {
945 220723774 : this->dxyzdxi_map[p].add_scaled (side_point, this->dpsidxi_map[i][p]);
946 220723774 : this->dxyzdeta_map[p].add_scaled(side_point, this->dpsideta_map[i][p]);
947 : }
948 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
949 148256542 : if (calculate_d2xyz)
950 : {
951 9651360 : this->d2xyzdxi2_map[p].add_scaled (side_point, this->d2psidxi2_map[i][p]);
952 9651360 : this->d2xyzdxideta_map[p].add_scaled(side_point, this->d2psidxideta_map[i][p]);
953 9651360 : 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 4889954 : if (calculate_dxyz)
961 : {
962 26130931 : for (unsigned int p=0; p<n_qp; p++)
963 : {
964 26646531 : const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
965 21290129 : this->normals[p] = n.unit();
966 26646531 : this->tangents[p][0] = this->dxyzdxi_map[p].unit();
967 26646531 : this->tangents[p][1] = n.cross(this->dxyzdxi_map[p]).unit();
968 :
969 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
970 21290129 : 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 655440 : const Real L = -this->d2xyzdxi2_map[p] * this->normals[p];
979 436960 : const Real M = -this->d2xyzdxideta_map[p] * this->normals[p];
980 436960 : const Real N = -this->d2xyzdeta2_map[p] * this->normals[p];
981 436960 : const Real E = this->dxyzdxi_map[p].norm_sq();
982 436960 : const Real F = this->dxyzdxi_map[p] * this->dxyzdeta_map[p];
983 218480 : const Real G = this->dxyzdeta_map[p].norm_sq();
984 :
985 873920 : const Real numerator = E*N -2.*F*M + G*L;
986 873920 : const Real denominator = E*G - F*F;
987 218480 : libmesh_assert_not_equal_to (denominator, 0.);
988 1092400 : 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 26130931 : for (unsigned int p=0; p<n_qp; p++)
996 : {
997 21290129 : const Real g11 = (dxdxi_map(p)*dxdxi_map(p) +
998 21290129 : dydxi_map(p)*dydxi_map(p) +
999 21290129 : dzdxi_map(p)*dzdxi_map(p));
1000 :
1001 21290129 : const Real g12 = (dxdxi_map(p)*dxdeta_map(p) +
1002 21290129 : dydxi_map(p)*dydeta_map(p) +
1003 21290129 : dzdxi_map(p)*dzdeta_map(p));
1004 :
1005 5356402 : const Real g21 = g12;
1006 :
1007 26646531 : const Real g22 = (dxdeta_map(p)*dxdeta_map(p) +
1008 21290129 : dydeta_map(p)*dydeta_map(p) +
1009 21290129 : dzdeta_map(p)*dzdeta_map(p));
1010 :
1011 :
1012 21290129 : const Real the_jac = std::sqrt(g11*g22 - g12*g21);
1013 :
1014 5356402 : libmesh_assert_greater (the_jac, 0.);
1015 :
1016 32002933 : this->JxW[p] = the_jac*qw[p];
1017 : }
1018 : }
1019 :
1020 : // done computing the map
1021 1227609 : break;
1022 : }
1023 :
1024 :
1025 0 : default:
1026 0 : libmesh_error_msg("Invalid dimension dim = " << dim);
1027 : }
1028 9549189 : }
1029 :
1030 :
1031 :
1032 :
1033 50934 : void FEMap::compute_edge_map(int dim,
1034 : const std::vector<Real> & qw,
1035 : const Elem * edge)
1036 : {
1037 13087 : libmesh_assert(edge);
1038 :
1039 50934 : 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 0 : this->compute_face_map(dim, qw, edge);
1045 0 : return;
1046 : }
1047 :
1048 13087 : libmesh_assert_equal_to (dim, 3); // 1D is unnecessary and currently unsupported
1049 :
1050 39261 : LOG_ONE_ENTITY_SCOPE("compute_edge_map()", "FEMap");
1051 :
1052 : // We're calculating now!
1053 13087 : this->determine_calculations();
1054 :
1055 : // The number of quadrature points.
1056 26174 : const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1057 :
1058 : // Resize the vectors to hold data at the quadrature points
1059 50934 : if (calculate_xyz)
1060 50934 : this->xyz.resize(n_qp);
1061 50934 : if (calculate_dxyz)
1062 : {
1063 50934 : this->dxyzdxi_map.resize(n_qp);
1064 50934 : this->dxyzdeta_map.resize(n_qp);
1065 50934 : this->tangents.resize(n_qp);
1066 50934 : this->normals.resize(n_qp);
1067 50934 : this->JxW.resize(n_qp);
1068 : }
1069 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1070 50934 : if (calculate_d2xyz)
1071 : {
1072 0 : this->d2xyzdxi2_map.resize(n_qp);
1073 0 : this->d2xyzdxideta_map.resize(n_qp);
1074 0 : this->d2xyzdeta2_map.resize(n_qp);
1075 0 : this->curvatures.resize(n_qp);
1076 : }
1077 : #endif
1078 :
1079 : // Clear the entities that will be summed
1080 245412 : for (unsigned int p=0; p<n_qp; p++)
1081 : {
1082 194478 : if (calculate_xyz)
1083 194478 : this->xyz[p].zero();
1084 194478 : if (calculate_dxyz)
1085 : {
1086 244446 : this->tangents[p].resize(1);
1087 99936 : this->dxyzdxi_map[p].zero();
1088 99936 : this->dxyzdeta_map[p].zero();
1089 : }
1090 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1091 194478 : if (calculate_d2xyz)
1092 : {
1093 0 : this->d2xyzdxi2_map[p].zero();
1094 0 : this->d2xyzdxideta_map[p].zero();
1095 0 : this->d2xyzdeta2_map[p].zero();
1096 : }
1097 : #endif
1098 : }
1099 :
1100 : // compute x, dxdxi at the quadrature points
1101 177550 : for (unsigned int i=0,
1102 26174 : psi_map_size=cast_int<unsigned int>(psi_map.size());
1103 203724 : i != psi_map_size; i++) // sum over the nodes
1104 : {
1105 78516 : const Point & edge_point = edge->point(i);
1106 :
1107 736200 : for (unsigned int p=0; p<n_qp; p++) // for each quadrature point...
1108 : {
1109 583410 : if (calculate_xyz)
1110 733308 : this->xyz[p].add_scaled (edge_point, this->psi_map[i][p]);
1111 583410 : if (calculate_dxyz)
1112 883206 : this->dxyzdxi_map[p].add_scaled (edge_point, this->dpsidxi_map[i][p]);
1113 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1114 583410 : if (calculate_d2xyz)
1115 0 : 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 50934 : if (calculate_dxyz)
1123 245412 : for (unsigned int p=0; p<n_qp; p++)
1124 : {
1125 : // const Point n = this->dxyzdxi_map[p].cross(this->dxyzdeta_map[p]);
1126 244446 : this->tangents[p][0] = this->dxyzdxi_map[p].unit();
1127 :
1128 : // compute the jacobian at the quadrature points
1129 194478 : const Real the_jac = std::sqrt(this->dxdxi_map(p)*this->dxdxi_map(p) +
1130 194478 : this->dydxi_map(p)*this->dydxi_map(p) +
1131 194478 : this->dzdxi_map(p)*this->dzdxi_map(p));
1132 :
1133 49968 : libmesh_assert_greater (the_jac, 0.);
1134 :
1135 294414 : this->JxW[p] = the_jac*qw[p];
1136 : }
1137 : }
1138 :
1139 :
1140 : // Explicit FEMap Instantiations
1141 0 : FACE_EDGE_SHAPE_ERROR(0,init_face_shape_functions)
1142 : template LIBMESH_EXPORT void FEMap::init_face_shape_functions<1>(const std::vector<Point> &, const Elem *);
1143 : template LIBMESH_EXPORT void FEMap::init_face_shape_functions<2>(const std::vector<Point> &, const Elem *);
1144 : template LIBMESH_EXPORT void FEMap::init_face_shape_functions<3>(const std::vector<Point> &, const Elem *);
1145 :
1146 0 : FACE_EDGE_SHAPE_ERROR(0,init_edge_shape_functions)
1147 : template LIBMESH_EXPORT void FEMap::init_edge_shape_functions<1>(const std::vector<Point> &, const Elem *);
1148 : template LIBMESH_EXPORT void FEMap::init_edge_shape_functions<2>(const std::vector<Point> &, const Elem *);
1149 : template 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) \
1154 : template LIBMESH_EXPORT void FE<1,_type>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const); \
1155 : template LIBMESH_EXPORT void FE<1,_type>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &); \
1156 : template LIBMESH_EXPORT void FE<2,_type>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const); \
1157 : template LIBMESH_EXPORT void FE<2,_type>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &); \
1158 : template LIBMESH_EXPORT void FE<2,_type>::edge_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &); \
1159 : template LIBMESH_EXPORT void FE<2,_type>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const); \
1160 : template LIBMESH_EXPORT void FE<3,_type>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const); \
1161 : template LIBMESH_EXPORT void FE<3,_type>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &); \
1162 : template LIBMESH_EXPORT void FE<3,_type>::edge_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &); \
1163 : template LIBMESH_EXPORT void FE<3,_type>::edge_reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const)
1164 :
1165 : REINIT_AND_SIDE_MAPS(LAGRANGE);
1166 : REINIT_AND_SIDE_MAPS(LAGRANGE_VEC);
1167 : REINIT_AND_SIDE_MAPS(L2_LAGRANGE);
1168 : REINIT_AND_SIDE_MAPS(L2_LAGRANGE_GLL);
1169 : REINIT_AND_SIDE_MAPS(L2_LAGRANGE_VEC);
1170 : REINIT_AND_SIDE_MAPS(HIERARCHIC);
1171 : REINIT_AND_SIDE_MAPS(HIERARCHIC_VEC);
1172 : REINIT_AND_SIDE_MAPS(L2_HIERARCHIC);
1173 : REINIT_AND_SIDE_MAPS(L2_HIERARCHIC_VEC);
1174 : REINIT_AND_SIDE_MAPS(SIDE_HIERARCHIC);
1175 : REINIT_AND_SIDE_MAPS(CLOUGH);
1176 : REINIT_AND_SIDE_MAPS(HERMITE);
1177 : REINIT_AND_SIDE_MAPS(MONOMIAL);
1178 : REINIT_AND_SIDE_MAPS(MONOMIAL_VEC);
1179 : REINIT_AND_SIDE_MAPS(SCALAR);
1180 : REINIT_AND_SIDE_MAPS(XYZ);
1181 :
1182 : #ifdef LIBMESH_ENABLE_HIGHER_ORDER_SHAPES
1183 : REINIT_AND_SIDE_MAPS(BERNSTEIN);
1184 : REINIT_AND_SIDE_MAPS(SZABAB);
1185 : REINIT_AND_SIDE_MAPS(RATIONAL_BERNSTEIN);
1186 : #endif
1187 :
1188 : template LIBMESH_EXPORT void FE<2,SUBDIVISION>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1189 : template LIBMESH_EXPORT void FE<2,SUBDIVISION>::edge_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1190 :
1191 : template LIBMESH_EXPORT void FE<2,NEDELEC_ONE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1192 : template LIBMESH_EXPORT void FE<2,NEDELEC_ONE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1193 : template LIBMESH_EXPORT void FE<2,NEDELEC_ONE>::edge_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1194 : 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);
1195 :
1196 : template LIBMESH_EXPORT void FE<2,RAVIART_THOMAS>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1197 : template LIBMESH_EXPORT void FE<2,RAVIART_THOMAS>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1198 : template LIBMESH_EXPORT void FE<2,RAVIART_THOMAS>::edge_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1199 : 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);
1200 :
1201 : 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);
1202 : 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> &);
1203 : 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> &);
1204 : 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);
1205 :
1206 : template LIBMESH_EXPORT void FE<3,NEDELEC_ONE>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1207 : template LIBMESH_EXPORT void FE<3,NEDELEC_ONE>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1208 : template LIBMESH_EXPORT void FE<3,NEDELEC_ONE>::edge_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1209 : 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);
1210 :
1211 : template LIBMESH_EXPORT void FE<3,RAVIART_THOMAS>::reinit(Elem const *, unsigned int, Real, const std::vector<Point> * const, const std::vector<Real> * const);
1212 : template LIBMESH_EXPORT void FE<3,RAVIART_THOMAS>::side_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1213 : template LIBMESH_EXPORT void FE<3,RAVIART_THOMAS>::edge_map(Elem const *, Elem const *, const unsigned int, const std::vector<Point> &, std::vector<Point> &);
1214 : 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);
1215 :
1216 : 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);
1217 : 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> &);
1218 : 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> &);
1219 : 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);
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
|
