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 : #ifndef LIBMESH_FE_MAP_H
21 : #define LIBMESH_FE_MAP_H
22 :
23 : // libMesh includes
24 : #include "libmesh/reference_counted_object.h"
25 : #include "libmesh/point.h"
26 : #include "libmesh/vector_value.h"
27 : #include "libmesh/fe_type.h"
28 :
29 : // C++ includes
30 : #include <memory>
31 :
32 : namespace libMesh
33 : {
34 :
35 : // forward declarations
36 : class Elem;
37 : class Node;
38 :
39 : /**
40 : * Class contained in FE that encapsulates mapping (i.e. from physical
41 : * space to reference space and vice-versa) quantities and operations.
42 : *
43 : * \author Paul Bauman
44 : * \date 2012
45 : * \brief Computes finite element mapping function values, gradients, etc.
46 : */
47 : class FEMap
48 : {
49 : public:
50 :
51 : FEMap(Real jtol = 0);
52 44401686 : virtual ~FEMap() = default;
53 :
54 : static std::unique_ptr<FEMap> build(FEType fe_type);
55 :
56 : static FEFamily map_fe_type(const Elem & elem);
57 :
58 : template<unsigned int Dim>
59 : void init_reference_to_physical_map(const std::vector<Point> & qp,
60 : const Elem * elem);
61 :
62 : // Non-templated version for runtime selection
63 : void init_reference_to_physical_map(unsigned int dim,
64 : const std::vector<Point> & qp,
65 : const Elem * elem);
66 :
67 : /**
68 : * Compute the jacobian and some other additional data fields at the
69 : * single point with index p. Takes the integration weights as
70 : * input, along with a pointer to the element and a list of points
71 : * that contribute to the element. Also takes a boolean flag
72 : * telling whether second derivatives should actually be computed.
73 : */
74 : void compute_single_point_map(const unsigned int dim,
75 : const std::vector<Real> & qw,
76 : const Elem * elem,
77 : unsigned int p,
78 : const std::vector<const Node *> & elem_nodes,
79 : bool compute_second_derivatives);
80 :
81 : /**
82 : * Compute the jacobian and some other additional
83 : * data fields. Takes the integration weights
84 : * as input, along with a pointer to the element.
85 : * The element is assumed to have a constant Jacobian
86 : */
87 : virtual void compute_affine_map(const unsigned int dim,
88 : const std::vector<Real> & qw,
89 : const Elem * elem);
90 :
91 : /**
92 : * Assign a fake jacobian and some other additional data fields.
93 : * Takes the integration weights as input. For use on non-element
94 : * evaluations.
95 : */
96 : virtual void compute_null_map(const unsigned int dim,
97 : const std::vector<Real> & qw);
98 :
99 :
100 : /**
101 : * Compute the jacobian and some other additional
102 : * data fields. Takes the integration weights
103 : * as input, along with a pointer to the element.
104 : * Also takes a boolean parameter indicating whether second
105 : * derivatives need to be calculated, allowing us to potentially
106 : * skip unnecessary, expensive computations.
107 : */
108 : virtual void compute_map(const unsigned int dim,
109 : const std::vector<Real> & qw,
110 : const Elem * elem,
111 : bool calculate_d2phi);
112 :
113 : /**
114 : * Same as compute_map, but for a side. Useful for boundary integration.
115 : */
116 : virtual void compute_face_map(int dim,
117 : const std::vector<Real> & qw,
118 : const Elem * side);
119 :
120 : /**
121 : * Same as before, but for an edge. Useful for some projections.
122 : */
123 : void compute_edge_map(int dim,
124 : const std::vector<Real> & qw,
125 : const Elem * side);
126 :
127 : /**
128 : * Initializes the reference to physical element map for a side.
129 : * This is used for boundary integration.
130 : */
131 : template<unsigned int Dim>
132 : void init_face_shape_functions(const std::vector<Point> & qp,
133 : const Elem * side);
134 :
135 : /**
136 : * Same as before, but for an edge. This is used for some projection
137 : * operators.
138 : */
139 : template<unsigned int Dim>
140 : void init_edge_shape_functions(const std::vector<Point> & qp,
141 : const Elem * edge);
142 :
143 : /**
144 : * \returns The location (in physical space) of the point
145 : * \p p located on the reference element.
146 : */
147 : static Point map (const unsigned int dim,
148 : const Elem * elem,
149 : const Point & reference_point);
150 :
151 : /**
152 : * \returns component \p j of d(xyz)/d(xi eta zeta) (in physical
153 : * space) of the point \p p located on the reference element.
154 : */
155 : static Point map_deriv (const unsigned int dim,
156 : const Elem * elem,
157 : const unsigned int j,
158 : const Point & reference_point);
159 :
160 : /**
161 : * \returns The location (on the reference element) of the
162 : * point \p p located in physical space. This function requires
163 : * inverting the (possibly nonlinear) transformation map, so
164 : * it is not trivial. The optional parameter \p tolerance defines
165 : * how close is "good enough." The map inversion iteration
166 : * computes the sequence \f$ \{ p_n \} \f$, and the iteration is
167 : * terminated when \f$ \|p - p_n\| < \mbox{\texttt{tolerance}} \f$
168 : *
169 : * When secure == true, the following checks are enabled:
170 : *
171 : * In DEBUG mode only:
172 : * .) dim==1,2: throw an error if det(J) <= 0 for any Newton iteration.
173 : * .) Print warning for every iteration beyond max_cnt in which the Newton scheme has not converged.
174 : *
175 : * In !DEBUG mode only:
176 : * .) Print a _single_ warning (1 warning for the entire simulation)
177 : * if the Newton scheme ever requires more than max_cnt iterations.
178 : *
179 : * In both DEBUG and !DEBUG modes:
180 : * .) dim==3: Throw an exception for singular Jacobian.
181 : * .) Throw an error if the Newton iteration has not converged in 2*max_cnt iterations.
182 : *
183 : * In addition to the checks above, the "extra_checks" parameter can
184 : * be used to turn on some additional tests. In particular, when
185 : * extra_checks == true *and* compiled in DEBUG mode:
186 : * .) Print a warning if p != map(inverse_map(p)) to within tolerance.
187 : * .) Print a warning if the inverse-mapped point is not on the reference element to within tolerance.
188 : */
189 : static Point inverse_map (const unsigned int dim,
190 : const Elem * elem,
191 : const Point & p,
192 : const Real tolerance = TOLERANCE,
193 : const bool secure = true,
194 : const bool extra_checks = true);
195 :
196 : /**
197 : * Takes a number points in physical space (in the \p
198 : * physical_points vector) and finds their location on the reference
199 : * element for the input element \p elem. The values on the
200 : * reference element are returned in the vector \p
201 : * reference_points. The other parameters have the same meaning
202 : * as the single Point version of inverse_map() above.
203 : */
204 : static void inverse_map (unsigned int dim,
205 : const Elem * elem,
206 : const std::vector<Point> & physical_points,
207 : std::vector<Point> & reference_points,
208 : const Real tolerance = TOLERANCE,
209 : const bool secure = true,
210 : const bool extra_checks = true);
211 :
212 : /**
213 : * \returns The \p xyz spatial locations of the quadrature
214 : * points on the element.
215 : */
216 3680850 : const std::vector<Point> & get_xyz() const
217 3680850 : { libmesh_assert(!calculations_started || calculate_xyz);
218 36265817 : calculate_xyz = true; return xyz; }
219 :
220 : /**
221 : * \returns The element Jacobian for each quadrature point.
222 : */
223 239195 : const std::vector<Real> & get_jacobian() const
224 239195 : { libmesh_assert(!calculations_started || calculate_dxyz);
225 1120237 : calculate_dxyz = true; return jac; }
226 :
227 : /**
228 : * \returns The element Jacobian times the quadrature weight for
229 : * each quadrature point.
230 : */
231 24 : const std::vector<Real> & get_JxW() const
232 24 : { libmesh_assert(!calculations_started || calculate_dxyz);
233 852 : calculate_dxyz = true; return JxW; }
234 :
235 : /**
236 : * \returns The element tangents in xi-direction at the quadrature
237 : * points.
238 : */
239 173793 : const std::vector<RealGradient> & get_dxyzdxi() const
240 173793 : { libmesh_assert(!calculations_started || calculate_dxyz);
241 1938067 : calculate_dxyz = true; return dxyzdxi_map; }
242 :
243 : /**
244 : * \returns The element tangents in eta-direction at the quadrature
245 : * points.
246 : */
247 173793 : const std::vector<RealGradient> & get_dxyzdeta() const
248 173793 : { libmesh_assert(!calculations_started || calculate_dxyz);
249 173797 : calculate_dxyz = true; return dxyzdeta_map; }
250 :
251 : /**
252 : * \returns The element tangents in zeta-direction at the quadrature
253 : * points.
254 : */
255 131209 : const std::vector<RealGradient> & get_dxyzdzeta() const
256 131209 : { libmesh_assert(!calculations_started || calculate_dxyz);
257 131209 : calculate_dxyz = true; return dxyzdzeta_map; }
258 :
259 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
260 :
261 : /**
262 : * \returns The second partial derivatives in xi.
263 : */
264 10 : const std::vector<RealGradient> & get_d2xyzdxi2() const
265 10 : { libmesh_assert(!calculations_started || calculate_d2xyz);
266 14 : calculate_d2xyz = true; return d2xyzdxi2_map; }
267 :
268 : /**
269 : * \returns The second partial derivatives in eta.
270 : */
271 10 : const std::vector<RealGradient> & get_d2xyzdeta2() const
272 10 : { libmesh_assert(!calculations_started || calculate_d2xyz);
273 10 : calculate_d2xyz = true; return d2xyzdeta2_map; }
274 :
275 : /**
276 : * \returns The second partial derivatives in zeta.
277 : */
278 0 : const std::vector<RealGradient> & get_d2xyzdzeta2() const
279 0 : { libmesh_assert(!calculations_started || calculate_d2xyz);
280 0 : calculate_d2xyz = true; return d2xyzdzeta2_map; }
281 :
282 : /**
283 : * \returns The second partial derivatives in xi-eta.
284 : */
285 10 : const std::vector<RealGradient> & get_d2xyzdxideta() const
286 10 : { libmesh_assert(!calculations_started || calculate_d2xyz);
287 10 : calculate_d2xyz = true; return d2xyzdxideta_map; }
288 :
289 : /**
290 : * \returns The second partial derivatives in xi-zeta.
291 : */
292 0 : const std::vector<RealGradient> & get_d2xyzdxidzeta() const
293 0 : { libmesh_assert(!calculations_started || calculate_d2xyz);
294 0 : calculate_d2xyz = true; return d2xyzdxidzeta_map; }
295 :
296 : /**
297 : * \returns The second partial derivatives in eta-zeta.
298 : */
299 0 : const std::vector<RealGradient> & get_d2xyzdetadzeta() const
300 0 : { libmesh_assert(!calculations_started || calculate_d2xyz);
301 0 : calculate_d2xyz = true; return d2xyzdetadzeta_map; }
302 :
303 : #endif
304 :
305 : /**
306 : * \returns The dxi/dx entry in the transformation
307 : * matrix from physical to local coordinates.
308 : */
309 32091763 : const std::vector<Real> & get_dxidx() const
310 32091763 : { libmesh_assert(!calculations_started || calculate_dxyz);
311 116864570 : calculate_dxyz = true; return dxidx_map; }
312 :
313 : /**
314 : * \returns The dxi/dy entry in the transformation
315 : * matrix from physical to local coordinates.
316 : */
317 7978099 : const std::vector<Real> & get_dxidy() const
318 7978099 : { libmesh_assert(!calculations_started || calculate_dxyz);
319 7979848 : calculate_dxyz = true; return dxidy_map; }
320 :
321 : /**
322 : * \returns The dxi/dz entry in the transformation
323 : * matrix from physical to local coordinates.
324 : */
325 7958215 : const std::vector<Real> & get_dxidz() const
326 7958215 : { libmesh_assert(!calculations_started || calculate_dxyz);
327 7959964 : calculate_dxyz = true; return dxidz_map; }
328 :
329 : /**
330 : * \returns The deta/dx entry in the transformation
331 : * matrix from physical to local coordinates.
332 : */
333 7956939 : const std::vector<Real> & get_detadx() const
334 7956939 : { libmesh_assert(!calculations_started || calculate_dxyz);
335 7958688 : calculate_dxyz = true; return detadx_map; }
336 :
337 : /**
338 : * \returns The deta/dy entry in the transformation
339 : * matrix from physical to local coordinates.
340 : */
341 7956939 : const std::vector<Real> & get_detady() const
342 7956939 : { libmesh_assert(!calculations_started || calculate_dxyz);
343 7958688 : calculate_dxyz = true; return detady_map; }
344 :
345 : /**
346 : * \returns The deta/dz entry in the transformation
347 : * matrix from physical to local coordinates.
348 : */
349 7937055 : const std::vector<Real> & get_detadz() const
350 7937055 : { libmesh_assert(!calculations_started || calculate_dxyz);
351 7938804 : calculate_dxyz = true; return detadz_map; }
352 :
353 : /**
354 : * \returns The dzeta/dx entry in the transformation
355 : * matrix from physical to local coordinates.
356 : */
357 1195084 : const std::vector<Real> & get_dzetadx() const
358 1195084 : { libmesh_assert(!calculations_started || calculate_dxyz);
359 1195084 : calculate_dxyz = true; return dzetadx_map; }
360 :
361 : /**
362 : * \returns The dzeta/dy entry in the transformation
363 : * matrix from physical to local coordinates.
364 : */
365 1195084 : const std::vector<Real> & get_dzetady() const
366 1195084 : { libmesh_assert(!calculations_started || calculate_dxyz);
367 1195084 : calculate_dxyz = true; return dzetady_map; }
368 :
369 : /**
370 : * \returns The dzeta/dz entry in the transformation
371 : * matrix from physical to local coordinates.
372 : */
373 1195084 : const std::vector<Real> & get_dzetadz() const
374 1195084 : { libmesh_assert(!calculations_started || calculate_dxyz);
375 1195084 : calculate_dxyz = true; return dzetadz_map; }
376 :
377 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
378 : /**
379 : * Second derivatives of "xi" reference coordinate wrt physical coordinates.
380 : */
381 4622849 : const std::vector<std::vector<Real>> & get_d2xidxyz2() const
382 4622849 : { libmesh_assert(!calculations_started || calculate_d2xyz);
383 10670199 : calculate_d2xyz = true; return d2xidxyz2_map; }
384 :
385 : /**
386 : * Second derivatives of "eta" reference coordinate wrt physical coordinates.
387 : */
388 513883 : const std::vector<std::vector<Real>> & get_d2etadxyz2() const
389 513883 : { libmesh_assert(!calculations_started || calculate_d2xyz);
390 513883 : calculate_d2xyz = true; return d2etadxyz2_map; }
391 :
392 : /**
393 : * Second derivatives of "zeta" reference coordinate wrt physical coordinates.
394 : */
395 432120 : const std::vector<std::vector<Real>> & get_d2zetadxyz2() const
396 432120 : { libmesh_assert(!calculations_started || calculate_d2xyz);
397 432120 : calculate_d2xyz = true; return d2zetadxyz2_map; }
398 : #endif
399 :
400 : /**
401 : * \returns The reference to physical map for the side/edge
402 : */
403 : const std::vector<std::vector<Real>> & get_psi() const
404 : { return psi_map; }
405 :
406 : /**
407 : * \returns The reference to physical map for the element
408 : */
409 : const std::vector<std::vector<Real>> & get_phi_map() const
410 : { libmesh_assert(!calculations_started || calculate_xyz);
411 : calculate_xyz = true; return phi_map; }
412 :
413 : /**
414 : * \returns The reference to physical map derivative
415 : */
416 13680 : const std::vector<std::vector<Real>> & get_dphidxi_map() const
417 13680 : { libmesh_assert(!calculations_started || calculate_dxyz);
418 164160 : calculate_dxyz = true; return dphidxi_map; }
419 :
420 : /**
421 : * \returns The reference to physical map derivative
422 : */
423 13680 : const std::vector<std::vector<Real>> & get_dphideta_map() const
424 13680 : { libmesh_assert(!calculations_started || calculate_dxyz);
425 164160 : calculate_dxyz = true; return dphideta_map; }
426 :
427 : /**
428 : * \returns The reference to physical map derivative
429 : */
430 13680 : const std::vector<std::vector<Real>> & get_dphidzeta_map() const
431 13680 : { libmesh_assert(!calculations_started || calculate_dxyz);
432 164160 : calculate_dxyz = true; return dphidzeta_map; }
433 :
434 : /**
435 : * \returns The tangent vectors for face integration.
436 : */
437 2 : const std::vector<std::vector<Point>> & get_tangents() const
438 2 : { libmesh_assert(!calculations_started || calculate_dxyz);
439 2 : calculate_dxyz = true; return tangents; }
440 :
441 : /**
442 : * \returns The outward pointing normal vectors for face integration.
443 : */
444 266957 : const std::vector<Point> & get_normals() const
445 266957 : { libmesh_assert(!calculations_started || calculate_dxyz);
446 748193 : calculate_dxyz = true; return normals; }
447 :
448 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
449 :
450 : /**
451 : * \returns The curvatures for use in face integration.
452 : */
453 0 : const std::vector<Real> & get_curvatures() const
454 0 : { libmesh_assert(!calculations_started || calculate_d2xyz);
455 0 : calculate_d2xyz = true; return curvatures;}
456 :
457 : #endif
458 :
459 : /**
460 : * Prints the Jacobian times the weight for each quadrature point.
461 : */
462 : void print_JxW(std::ostream & os) const;
463 :
464 : /**
465 : * Prints the spatial location of each quadrature point
466 : * (on the physical element).
467 : */
468 : void print_xyz(std::ostream & os) const;
469 :
470 : /* FIXME: PB: The public functions below break encapsulation! I'm not happy
471 : about it, but I don't have time to redo the infinite element routines.
472 : These are needed in inf_fe_boundary.C and inf_fe.C. A proper implementation
473 : would subclass this class and hide all the radial function business. */
474 : /**
475 : * \returns The reference to physical map for the side/edge
476 : */
477 2243689 : std::vector<std::vector<Real>> & get_psi()
478 2243689 : { return psi_map; }
479 :
480 : /**
481 : * \returns The reference to physical map derivative for the side/edge
482 : */
483 : std::vector<std::vector<Real>> & get_dpsidxi()
484 : { libmesh_assert(!calculations_started || calculate_dxyz);
485 : calculate_dxyz = true; return dpsidxi_map; }
486 :
487 : const std::vector<std::vector<Real>> & get_dpsidxi() const
488 : { libmesh_assert(!calculations_started || calculate_dxyz);
489 : calculate_dxyz = true; return dpsidxi_map; }
490 :
491 : /**
492 : * \returns The reference to physical map derivative for the side/edge
493 : */
494 : std::vector<std::vector<Real>> & get_dpsideta()
495 : { libmesh_assert(!calculations_started || calculate_dxyz);
496 : calculate_dxyz = true; return dpsideta_map; }
497 :
498 : const std::vector<std::vector<Real>> & get_dpsideta() const
499 : { libmesh_assert(!calculations_started || calculate_dxyz);
500 : calculate_dxyz = true; return dpsideta_map; }
501 :
502 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
503 :
504 : /**
505 : * \returns The reference to physical map 2nd derivative for the side/edge
506 : */
507 : std::vector<std::vector<Real>> & get_d2psidxi2()
508 : { libmesh_assert(!calculations_started || calculate_d2xyz);
509 : calculate_d2xyz = true; return d2psidxi2_map; }
510 :
511 : /**
512 : * \returns const reference to physical map 2nd derivative for the side/edge
513 : */
514 : const std::vector<std::vector<Real>> & get_d2psidxi2() const
515 : { libmesh_assert(!calculations_started || calculate_d2xyz);
516 : calculate_d2xyz = true; return d2psidxi2_map; }
517 :
518 : /**
519 : * \returns The reference to physical map 2nd derivative for the side/edge
520 : */
521 : std::vector<std::vector<Real>> & get_d2psidxideta()
522 : { libmesh_assert(!calculations_started || calculate_d2xyz);
523 : calculate_d2xyz = true; return d2psidxideta_map; }
524 :
525 : /**
526 : * \returns const reference to physical map 2nd derivative for the side/edge
527 : */
528 : const std::vector<std::vector<Real>> & get_d2psidxideta() const
529 : { libmesh_assert(!calculations_started || calculate_d2xyz);
530 : calculate_d2xyz = true; return d2psidxideta_map; }
531 :
532 : /**
533 : * \returns The reference to physical map 2nd derivative for the side/edge
534 : */
535 : std::vector<std::vector<Real>> & get_d2psideta2()
536 : { libmesh_assert(!calculations_started || calculate_d2xyz);
537 : calculate_d2xyz = true; return d2psideta2_map; }
538 :
539 :
540 : /**
541 : * \returns const reference to physical map 2nd derivative for the side/edge
542 : */
543 : const std::vector<std::vector<Real>> & get_d2psideta2() const
544 : { libmesh_assert(!calculations_started || calculate_d2xyz);
545 : calculate_d2xyz = true; return d2psideta2_map; }
546 :
547 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
548 :
549 : /**
550 : * \returns The reference to physical map for the element
551 : */
552 1152 : std::vector<std::vector<Real>> & get_phi_map()
553 1152 : { libmesh_assert(!calculations_started || calculate_xyz);
554 3968 : calculate_xyz = true; return phi_map; }
555 :
556 : /**
557 : * \returns The reference to physical map derivative
558 : */
559 256 : std::vector<std::vector<Real>> & get_dphidxi_map()
560 256 : { libmesh_assert(!calculations_started || calculate_dxyz);
561 3072 : calculate_dxyz = true; return dphidxi_map; }
562 :
563 : /**
564 : * \returns The reference to physical map derivative
565 : */
566 256 : std::vector<std::vector<Real>> & get_dphideta_map()
567 256 : { libmesh_assert(!calculations_started || calculate_dxyz);
568 3072 : calculate_dxyz = true; return dphideta_map; }
569 :
570 : /**
571 : * \returns The reference to physical map derivative
572 : */
573 : std::vector<std::vector<Real>> & get_dphidzeta_map()
574 : { libmesh_assert(!calculations_started || calculate_dxyz);
575 : calculate_dxyz = true; return dphidzeta_map; }
576 :
577 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
578 : /**
579 : * \returns The reference to physical map 2nd derivative
580 : */
581 256 : std::vector<std::vector<Real>> & get_d2phidxi2_map()
582 256 : { libmesh_assert(!calculations_started || calculate_d2xyz);
583 3072 : calculate_d2xyz = true; return d2phidxi2_map; }
584 :
585 : /**
586 : * \returns The reference to physical map 2nd derivative
587 : */
588 256 : std::vector<std::vector<Real>> & get_d2phidxideta_map()
589 256 : { libmesh_assert(!calculations_started || calculate_d2xyz);
590 3072 : calculate_d2xyz = true; return d2phidxideta_map; }
591 :
592 : /**
593 : * \returns The reference to physical map 2nd derivative
594 : */
595 : std::vector<std::vector<Real>> & get_d2phidxidzeta_map()
596 : { libmesh_assert(!calculations_started || calculate_d2xyz);
597 : calculate_d2xyz = true; return d2phidxidzeta_map; }
598 :
599 : /**
600 : * \returns The reference to physical map 2nd derivative
601 : */
602 256 : std::vector<std::vector<Real>> & get_d2phideta2_map()
603 256 : { libmesh_assert(!calculations_started || calculate_d2xyz);
604 3072 : calculate_d2xyz = true; return d2phideta2_map; }
605 :
606 : /**
607 : * \returns The reference to physical map 2nd derivative
608 : */
609 : std::vector<std::vector<Real>> & get_d2phidetadzeta_map()
610 : { libmesh_assert(!calculations_started || calculate_d2xyz);
611 : calculate_d2xyz = true; return d2phidetadzeta_map; }
612 :
613 : /**
614 : * \returns The reference to physical map 2nd derivative
615 : */
616 : std::vector<std::vector<Real>> & get_d2phidzeta2_map()
617 : { libmesh_assert(!calculations_started || calculate_d2xyz);
618 : calculate_d2xyz = true; return d2phidzeta2_map; }
619 : #endif
620 :
621 : /* FIXME: PB: This function breaks encapsulation! Needed in FE<>::reinit and
622 : InfFE<>::reinit. Not sure yet if the algorithm can be redone to avoid this. */
623 : /**
624 : * \returns Writable reference to the element Jacobian times
625 : * the quadrature weight for each quadrature point.
626 : */
627 14100075 : std::vector<Real> & get_JxW()
628 14100075 : { libmesh_assert(!calculations_started || calculate_dxyz);
629 121953199 : calculate_dxyz = true; return JxW; }
630 :
631 : /**
632 : * Allows the user to prerequest additional calculations in between
633 : * two calls to reinit();
634 : */
635 : void add_calculations();
636 :
637 : /**
638 : * Set the Jacobian tolerance used for determining when the mapping fails. The mapping is
639 : * determined to fail if jac <= jacobian_tolerance.
640 : */
641 32449 : void set_jacobian_tolerance(Real tol) { jacobian_tolerance = tol; }
642 :
643 : protected:
644 :
645 : /**
646 : * Determine which values are to be calculated
647 : */
648 21369291 : void determine_calculations()
649 : {
650 249976929 : calculations_started = true;
651 :
652 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
653 : // Second derivative calculations currently have first derivative
654 : // calculations as a prerequisite
655 249976929 : if (calculate_d2xyz)
656 57651080 : calculate_dxyz = true;
657 : #endif
658 21369291 : }
659 :
660 : /**
661 : * A utility function for use by compute_*_map
662 : */
663 : void resize_quadrature_map_vectors(const unsigned int dim, unsigned int n_qp);
664 :
665 : /**
666 : * Used in \p FEMap::compute_map(), which should be
667 : * be usable in derived classes, and therefore protected.
668 : *
669 : * \returns The x value of the pth entry of the dxzydxi_map.
670 : */
671 235373189 : Real dxdxi_map(const unsigned int p) const { return dxyzdxi_map[p](0); }
672 :
673 : /**
674 : * Used in \p FEMap::compute_map(), which should be
675 : * be usable in derived classes, and therefore protected.
676 : *
677 : * \returns The y value of the pth entry of the dxzydxi_map.
678 : */
679 127872328 : Real dydxi_map(const unsigned int p) const { return dxyzdxi_map[p](1); }
680 :
681 : /**
682 : * Used in \p FEMap::compute_map(), which should be
683 : * be usable in derived classes, and therefore protected.
684 : *
685 : * \returns The z value of the pth entry of the dxzydxi_map.
686 : */
687 127872328 : Real dzdxi_map(const unsigned int p) const { return dxyzdxi_map[p](2); }
688 :
689 : /**
690 : * Used in \p FEMap::compute_map(), which should be
691 : * be usable in derived classes, and therefore protected.
692 : *
693 : * \returns The x value of the pth entry of the dxzydeta_map.
694 : */
695 189424953 : Real dxdeta_map(const unsigned int p) const { return dxyzdeta_map[p](0); }
696 :
697 : /**
698 : * Used in \p FEMap::compute_map(), which should be
699 : * be usable in derived classes, and therefore protected.
700 : *
701 : * \returns The y value of the pth entry of the dxzydeta_map.
702 : */
703 184071153 : Real dydeta_map(const unsigned int p) const { return dxyzdeta_map[p](1); }
704 :
705 : /**
706 : * Used in \p FEMap::compute_map(), which should be
707 : * be usable in derived classes, and therefore protected.
708 : *
709 : * \returns The z value of the pth entry of the dxzydeta_map.
710 : */
711 157772107 : Real dzdeta_map(const unsigned int p) const { return dxyzdeta_map[p](2); }
712 :
713 : /**
714 : * Used in \p FEMap::compute_map(), which should be
715 : * be usable in derived classes, and therefore protected.
716 : *
717 : * \returns The x value of the pth entry of the dxzydzeta_map.
718 : */
719 28693497 : Real dxdzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](0); }
720 :
721 : /**
722 : * Used in \p FEMap::compute_map(), which should be
723 : * be usable in derived classes, and therefore protected.
724 : *
725 : * \returns The y value of the pth entry of the dxzydzeta_map.
726 : */
727 28693497 : Real dydzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](1); }
728 :
729 : /**
730 : * Used in \p FEMap::compute_map(), which should be
731 : * be usable in derived classes, and therefore protected.
732 : *
733 : * \returns The z value of the pth entry of the dxzydzeta_map.
734 : */
735 28693497 : Real dzdzeta_map(const unsigned int p) const { return dxyzdzeta_map[p](2); }
736 :
737 : /**
738 : * The spatial locations of the quadrature points
739 : */
740 : std::vector<Point> xyz;
741 :
742 : /**
743 : * Vector of partial derivatives:
744 : * d(x)/d(xi), d(y)/d(xi), d(z)/d(xi)
745 : */
746 : std::vector<RealGradient> dxyzdxi_map;
747 :
748 : /**
749 : * Vector of partial derivatives:
750 : * d(x)/d(eta), d(y)/d(eta), d(z)/d(eta)
751 : */
752 : std::vector<RealGradient> dxyzdeta_map;
753 :
754 : /**
755 : * Vector of partial derivatives:
756 : * d(x)/d(zeta), d(y)/d(zeta), d(z)/d(zeta)
757 : */
758 : std::vector<RealGradient> dxyzdzeta_map;
759 :
760 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
761 :
762 : /**
763 : * Vector of second partial derivatives in xi:
764 : * d^2(x)/d(xi)^2, d^2(y)/d(xi)^2, d^2(z)/d(xi)^2
765 : */
766 : std::vector<RealGradient> d2xyzdxi2_map;
767 :
768 : /**
769 : * Vector of mixed second partial derivatives in xi-eta:
770 : * d^2(x)/d(xi)d(eta) d^2(y)/d(xi)d(eta) d^2(z)/d(xi)d(eta)
771 : */
772 : std::vector<RealGradient> d2xyzdxideta_map;
773 :
774 : /**
775 : * Vector of second partial derivatives in eta:
776 : * d^2(x)/d(eta)^2
777 : */
778 : std::vector<RealGradient> d2xyzdeta2_map;
779 :
780 : /**
781 : * Vector of second partial derivatives in xi-zeta:
782 : * d^2(x)/d(xi)d(zeta), d^2(y)/d(xi)d(zeta), d^2(z)/d(xi)d(zeta)
783 : */
784 : std::vector<RealGradient> d2xyzdxidzeta_map;
785 :
786 : /**
787 : * Vector of mixed second partial derivatives in eta-zeta:
788 : * d^2(x)/d(eta)d(zeta) d^2(y)/d(eta)d(zeta) d^2(z)/d(eta)d(zeta)
789 : */
790 : std::vector<RealGradient> d2xyzdetadzeta_map;
791 :
792 : /**
793 : * Vector of second partial derivatives in zeta:
794 : * d^2(x)/d(zeta)^2
795 : */
796 : std::vector<RealGradient> d2xyzdzeta2_map;
797 :
798 : #endif //LIBMESH_ENABLE_SECOND_DERIVATIVES
799 :
800 : /**
801 : * Map for partial derivatives:
802 : * d(xi)/d(x). Needed for the Jacobian.
803 : */
804 : std::vector<Real> dxidx_map;
805 :
806 : /**
807 : * Map for partial derivatives:
808 : * d(xi)/d(y). Needed for the Jacobian.
809 : */
810 : std::vector<Real> dxidy_map;
811 :
812 : /**
813 : * Map for partial derivatives:
814 : * d(xi)/d(z). Needed for the Jacobian.
815 : */
816 : std::vector<Real> dxidz_map;
817 :
818 :
819 : /**
820 : * Map for partial derivatives:
821 : * d(eta)/d(x). Needed for the Jacobian.
822 : */
823 : std::vector<Real> detadx_map;
824 :
825 : /**
826 : * Map for partial derivatives:
827 : * d(eta)/d(y). Needed for the Jacobian.
828 : */
829 : std::vector<Real> detady_map;
830 :
831 : /**
832 : * Map for partial derivatives:
833 : * d(eta)/d(z). Needed for the Jacobian.
834 : */
835 : std::vector<Real> detadz_map;
836 :
837 :
838 : /**
839 : * Map for partial derivatives:
840 : * d(zeta)/d(x). Needed for the Jacobian.
841 : */
842 : std::vector<Real> dzetadx_map;
843 :
844 : /**
845 : * Map for partial derivatives:
846 : * d(zeta)/d(y). Needed for the Jacobian.
847 : */
848 : std::vector<Real> dzetady_map;
849 :
850 : /**
851 : * Map for partial derivatives:
852 : * d(zeta)/d(z). Needed for the Jacobian.
853 : */
854 : std::vector<Real> dzetadz_map;
855 :
856 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
857 : /**
858 : * Second derivatives of "xi" reference coordinate wrt physical coordinates.
859 : * At each qp: (xi_{xx}, xi_{xy}, xi_{xz}, xi_{yy}, xi_{yz}, xi_{zz})
860 : */
861 : std::vector<std::vector<Real>> d2xidxyz2_map;
862 :
863 : /**
864 : * Second derivatives of "eta" reference coordinate wrt physical coordinates.
865 : * At each qp: (eta_{xx}, eta_{xy}, eta_{xz}, eta_{yy}, eta_{yz}, eta_{zz})
866 : */
867 : std::vector<std::vector<Real>> d2etadxyz2_map;
868 :
869 : /**
870 : * Second derivatives of "zeta" reference coordinate wrt physical coordinates.
871 : * At each qp: (zeta_{xx}, zeta_{xy}, zeta_{xz}, zeta_{yy}, zeta_{yz}, zeta_{zz})
872 : */
873 : std::vector<std::vector<Real>> d2zetadxyz2_map;
874 : #endif
875 :
876 : /**
877 : * Map for the shape function phi.
878 : */
879 : std::vector<std::vector<Real>> phi_map;
880 :
881 : /**
882 : * Map for the derivative, d(phi)/d(xi).
883 : */
884 : std::vector<std::vector<Real>> dphidxi_map;
885 :
886 : /**
887 : * Map for the derivative, d(phi)/d(eta).
888 : */
889 : std::vector<std::vector<Real>> dphideta_map;
890 :
891 : /**
892 : * Map for the derivative, d(phi)/d(zeta).
893 : */
894 : std::vector<std::vector<Real>> dphidzeta_map;
895 :
896 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
897 :
898 : /**
899 : * Map for the second derivative, d^2(phi)/d(xi)^2.
900 : */
901 : std::vector<std::vector<Real>> d2phidxi2_map;
902 :
903 : /**
904 : * Map for the second derivative, d^2(phi)/d(xi)d(eta).
905 : */
906 : std::vector<std::vector<Real>> d2phidxideta_map;
907 :
908 : /**
909 : * Map for the second derivative, d^2(phi)/d(xi)d(zeta).
910 : */
911 : std::vector<std::vector<Real>> d2phidxidzeta_map;
912 :
913 : /**
914 : * Map for the second derivative, d^2(phi)/d(eta)^2.
915 : */
916 : std::vector<std::vector<Real>> d2phideta2_map;
917 :
918 : /**
919 : * Map for the second derivative, d^2(phi)/d(eta)d(zeta).
920 : */
921 : std::vector<std::vector<Real>> d2phidetadzeta_map;
922 :
923 : /**
924 : * Map for the second derivative, d^2(phi)/d(zeta)^2.
925 : */
926 : std::vector<std::vector<Real>> d2phidzeta2_map;
927 :
928 : #endif
929 :
930 : /**
931 : * Map for the side shape functions, psi.
932 : */
933 : std::vector<std::vector<Real>> psi_map;
934 :
935 : /**
936 : * Map for the derivative of the side functions,
937 : * d(psi)/d(xi).
938 : */
939 : std::vector<std::vector<Real>> dpsidxi_map;
940 :
941 : /**
942 : * Map for the derivative of the side function,
943 : * d(psi)/d(eta).
944 : */
945 : std::vector<std::vector<Real>> dpsideta_map;
946 :
947 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
948 :
949 : /**
950 : * Map for the second derivatives (in xi) of the
951 : * side shape functions. Useful for computing
952 : * the curvature at the quadrature points.
953 : */
954 : std::vector<std::vector<Real>> d2psidxi2_map;
955 :
956 : /**
957 : * Map for the second (cross) derivatives in xi, eta
958 : * of the side shape functions. Useful for
959 : * computing the curvature at the quadrature points.
960 : */
961 : std::vector<std::vector<Real>> d2psidxideta_map;
962 :
963 : /**
964 : * Map for the second derivatives (in eta) of the
965 : * side shape functions. Useful for computing the
966 : * curvature at the quadrature points.
967 : */
968 : std::vector<std::vector<Real>> d2psideta2_map;
969 :
970 : #endif
971 :
972 : /**
973 : * Tangent vectors on boundary at quadrature points.
974 : */
975 : std::vector<std::vector<Point>> tangents;
976 :
977 : /**
978 : * Normal vectors on boundary at quadrature points
979 : */
980 : std::vector<Point> normals;
981 :
982 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
983 : /**
984 : * The mean curvature (= one half the sum of the principal
985 : * curvatures) on the boundary at the quadrature points.
986 : * The mean curvature is a scalar value.
987 : */
988 : std::vector<Real> curvatures;
989 :
990 : #endif
991 :
992 : /**
993 : * Jacobian values at quadrature points
994 : */
995 : std::vector<Real> jac;
996 :
997 : /**
998 : * Jacobian*Weight values at quadrature points
999 : */
1000 : std::vector<Real> JxW;
1001 :
1002 : /**
1003 : * Have calculations with this object already been started?
1004 : * Then all get_* functions should already have been called.
1005 : */
1006 : mutable bool calculations_started;
1007 :
1008 : /**
1009 : * Should we calculate physical point locations?
1010 : */
1011 : mutable bool calculate_xyz;
1012 :
1013 : /**
1014 : * Should we calculate mapping gradients?
1015 : */
1016 : mutable bool calculate_dxyz;
1017 :
1018 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1019 :
1020 : /**
1021 : * Should we calculate mapping hessians?
1022 : */
1023 : mutable bool calculate_d2xyz;
1024 :
1025 : #endif
1026 :
1027 : /**
1028 : * The Jacobian tolerance used for determining when the mapping fails. The mapping is
1029 : * determined to fail if jac <= jacobian_tolerance. If not set by the user, this number
1030 : * defaults to 0
1031 : */
1032 : Real jacobian_tolerance;
1033 :
1034 : private:
1035 : /**
1036 : * A helper function used by FEMap::compute_single_point_map() to
1037 : * compute second derivatives of the inverse map.
1038 : */
1039 : void compute_inverse_map_second_derivs(unsigned p);
1040 :
1041 : /**
1042 : * Work vector for compute_affine_map()
1043 : */
1044 : std::vector<const Node *> _elem_nodes;
1045 :
1046 : /**
1047 : * A mutex for locking the error stream for failed point inversions
1048 : */
1049 : static Threads::spin_mutex _point_inv_err_mutex;
1050 : };
1051 :
1052 :
1053 :
1054 : inline void
1055 1015839 : FEMap::init_reference_to_physical_map(unsigned int dim,
1056 : const std::vector<Point> & qp,
1057 : const Elem * elem)
1058 : {
1059 1015839 : switch (dim)
1060 : {
1061 0 : case 0:
1062 0 : this->init_reference_to_physical_map<0>(qp, elem);
1063 0 : break;
1064 0 : case 1:
1065 0 : this->init_reference_to_physical_map<1>(qp, elem);
1066 0 : break;
1067 27327 : case 2:
1068 27327 : this->init_reference_to_physical_map<2>(qp, elem);
1069 27327 : break;
1070 988512 : case 3:
1071 988512 : this->init_reference_to_physical_map<3>(qp, elem);
1072 988512 : break;
1073 0 : default:
1074 0 : libmesh_error();
1075 : }
1076 1015839 : }
1077 :
1078 : }
1079 :
1080 : #endif // LIBMESH_FE_MAP_H
|
