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 : // Local includes
21 : #include "libmesh/fe.h"
22 : #include "libmesh/elem.h"
23 : #include "libmesh/libmesh_logging.h"
24 : #include "libmesh/fe_interface.h"
25 : #include "libmesh/fe_macro.h"
26 : #include "libmesh/fe_map.h"
27 : #include "libmesh/fe_xyz_map.h"
28 : #include "libmesh/inf_fe_map.h"
29 : #include "libmesh/mesh_subdivision_support.h"
30 : #include "libmesh/dense_matrix.h"
31 : #include "libmesh/dense_vector.h"
32 : #include "libmesh/tensor_value.h"
33 : #include "libmesh/enum_elem_type.h"
34 : #include "libmesh/int_range.h"
35 :
36 : // C++ includes
37 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
38 : #include <cmath> // for std::sqrt, std::abs
39 : #include <memory>
40 :
41 : namespace libMesh
42 : {
43 : Threads::spin_mutex FEMap::_point_inv_err_mutex;
44 :
45 : FEFamily
46 1206931761 : FEMap::map_fe_type(const Elem & elem)
47 : {
48 1206931761 : switch (elem.mapping_type())
49 : {
50 631683 : case RATIONAL_BERNSTEIN_MAP:
51 631683 : return RATIONAL_BERNSTEIN;
52 1199952626 : case LAGRANGE_MAP:
53 1199952626 : return LAGRANGE;
54 0 : default:
55 0 : libmesh_error_msg("Unknown mapping type " << elem.mapping_type());
56 : }
57 : return LAGRANGE;
58 : }
59 :
60 :
61 :
62 : // Constructor
63 17285617 : FEMap::FEMap(Real jtol) :
64 15398621 : calculations_started(false),
65 15398621 : calculate_xyz(false),
66 15398621 : calculate_dxyz(false),
67 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
68 15398621 : calculate_d2xyz(false),
69 : #endif
70 18228569 : jacobian_tolerance(jtol)
71 17285617 : {}
72 :
73 :
74 :
75 17205914 : std::unique_ptr<FEMap> FEMap::build( FEType fe_type )
76 : {
77 17205914 : switch( fe_type.family )
78 : {
79 859376 : case XYZ:
80 859376 : return std::make_unique<FEXYZMap>();
81 :
82 16346538 : default:
83 16346538 : return std::make_unique<FEMap>();
84 : }
85 : }
86 :
87 :
88 :
89 24799413 : void FEMap::add_calculations()
90 : {
91 24799413 : this->calculations_started = false;
92 22558015 : this->phi_map.clear();
93 22558015 : this->dphidxi_map.clear();
94 22558015 : this->dphideta_map.clear();
95 22558015 : this->dphidzeta_map.clear();
96 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
97 22558015 : this->d2phidxi2_map.clear();
98 22558015 : this->d2phidxideta_map.clear();
99 22558015 : this->d2phideta2_map.clear();
100 22558015 : this->d2phidxidzeta_map.clear();
101 22558015 : this->d2phidetadzeta_map.clear();
102 22558015 : this->d2phidzeta2_map.clear();
103 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
104 24799413 : }
105 :
106 :
107 :
108 : template<unsigned int Dim>
109 48044778 : void FEMap::init_reference_to_physical_map(const std::vector<Point> & qp,
110 : const Elem * elem)
111 : {
112 : // Start logging the reference->physical map initialization
113 8462308 : LOG_SCOPE("init_reference_to_physical_map()", "FEMap");
114 :
115 : // We're calculating now!
116 4231154 : this->determine_calculations();
117 :
118 : // The number of quadrature points.
119 8461841 : const std::size_t n_qp = qp.size();
120 :
121 : // The element type and order to use in
122 : // the map
123 48044778 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
124 48044778 : const FEType map_fe_type(elem->default_order(), mapping_family);
125 :
126 : // Number of shape functions used to construct the map
127 : // (Lagrange shape functions are used for mapping)
128 : // Do not consider the Elem::p_level(), if any, when computing the
129 : // number of shape functions.
130 : const unsigned int n_mapping_shape_functions =
131 48044778 : FEInterface::n_shape_functions (map_fe_type, /*extra_order=*/0, elem);
132 :
133 : // Maybe we already have correctly-sized data? Check data sizes,
134 : // and get ready to break out of a "loop" if all these resize()
135 : // calls are redundant.
136 4231154 : unsigned int old_n_qp = 0;
137 :
138 : do
139 : {
140 48044778 : if (calculate_xyz)
141 : {
142 30667156 : if (this->phi_map.size() == n_mapping_shape_functions)
143 : {
144 2976537 : old_n_qp = n_mapping_shape_functions ? this->phi_map[0].size() : 0;
145 264086 : break;
146 : }
147 25157972 : this->phi_map.resize (n_mapping_shape_functions);
148 : }
149 : if (Dim > 0)
150 : {
151 44869461 : if (calculate_dxyz)
152 : {
153 39688871 : if (this->dphidxi_map.size() == n_mapping_shape_functions)
154 : {
155 7256312 : old_n_qp = n_mapping_shape_functions ? this->dphidxi_map[0].size() : 0;
156 659727 : break;
157 : }
158 29202645 : this->dphidxi_map.resize (n_mapping_shape_functions);
159 : }
160 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
161 : // Don't waste time considering early break here; if we
162 : // asked for d2xyz then we surely asked for dxyz too and
163 : // considered early break above.
164 37613149 : if (calculate_d2xyz)
165 1899924 : this->d2phidxi2_map.resize (n_mapping_shape_functions);
166 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
167 : }
168 :
169 : if (Dim > 1)
170 : {
171 34921511 : if (calculate_dxyz)
172 28831246 : this->dphideta_map.resize (n_mapping_shape_functions);
173 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
174 34921511 : if (calculate_d2xyz)
175 : {
176 1614571 : this->d2phidxideta_map.resize (n_mapping_shape_functions);
177 1614571 : this->d2phideta2_map.resize (n_mapping_shape_functions);
178 : }
179 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
180 : }
181 :
182 : if (Dim > 2)
183 : {
184 18889640 : if (calculate_dxyz)
185 16578132 : this->dphidzeta_map.resize (n_mapping_shape_functions);
186 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
187 18889640 : if (calculate_d2xyz)
188 : {
189 1455611 : this->d2phidxidzeta_map.resize (n_mapping_shape_functions);
190 1455611 : this->d2phidetadzeta_map.resize (n_mapping_shape_functions);
191 1455611 : this->d2phidzeta2_map.resize (n_mapping_shape_functions);
192 : }
193 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
194 : }
195 : }
196 : while (false);
197 :
198 48044778 : if (old_n_qp != n_qp)
199 394690868 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
200 : {
201 356851750 : if (calculate_xyz)
202 272447496 : this->phi_map[i].resize (n_qp);
203 : if (Dim > 0)
204 : {
205 356652970 : if (calculate_dxyz)
206 310651177 : this->dphidxi_map[i].resize (n_qp);
207 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
208 356652970 : if (calculate_d2xyz)
209 : {
210 17234763 : this->d2phidxi2_map[i].resize (n_qp);
211 : if (Dim > 1)
212 : {
213 16642094 : this->d2phidxideta_map[i].resize (n_qp);
214 16642094 : this->d2phideta2_map[i].resize (n_qp);
215 : }
216 : if (Dim > 2)
217 : {
218 15244901 : this->d2phidxidzeta_map[i].resize (n_qp);
219 15244901 : this->d2phidetadzeta_map[i].resize (n_qp);
220 15244901 : this->d2phidzeta2_map[i].resize (n_qp);
221 : }
222 : }
223 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
224 :
225 351130072 : if (Dim > 1 && calculate_dxyz)
226 309780094 : this->dphideta_map[i].resize (n_qp);
227 :
228 256142339 : if (Dim > 2 && calculate_dxyz)
229 237535591 : this->dphidzeta_map[i].resize (n_qp);
230 : }
231 : }
232 :
233 : // Optimize for the *linear* geometric elements case:
234 48044778 : bool is_linear = elem->is_linear();
235 :
236 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
237 48044778 : FEInterface::shape_deriv_function(map_fe_type, elem);
238 :
239 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
240 : FEInterface::shape_second_deriv_ptr shape_second_deriv_ptr =
241 48044778 : FEInterface::shape_second_deriv_function(map_fe_type, elem);
242 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
243 :
244 48044778 : if (calculate_xyz)
245 28134509 : FEInterface::all_shapes(Dim, map_fe_type, elem, qp, this->phi_map, false);
246 :
247 : switch (Dim)
248 : {
249 : //------------------------------------------------------------
250 : // 0D
251 : case 0:
252 : {
253 : // No mapping derivatives here
254 19049 : break;
255 : }
256 :
257 : //------------------------------------------------------------
258 : // 1D
259 2512411 : case 1:
260 : {
261 : // Compute gradients of mapping shape functions i at quadrature points p
262 :
263 : // TODO: optimizations for the RATIONAL_BERNSTEIN+linear case
264 2739907 : if (is_linear)
265 : {
266 7683882 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
267 : {
268 5122588 : if (calculate_dxyz)
269 : {
270 583196 : this->dphidxi_map[i][0] =
271 571494 : shape_deriv_ptr(map_fe_type, elem, i, 0, qp[0], false);
272 2237844 : for (std::size_t p=1; p<n_qp; p++)
273 1697496 : this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
274 : }
275 :
276 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
277 5122588 : if (calculate_d2xyz)
278 : {
279 577880 : this->d2phidxi2_map[i][0] =
280 566564 : shape_second_deriv_ptr(map_fe_type, elem, i, 0, qp[0], false);
281 2228492 : for (std::size_t p=1; p<n_qp; p++)
282 1692740 : this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
283 : }
284 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
285 : }
286 : }
287 : else
288 : {
289 178613 : if (calculate_dxyz)
290 : {
291 177635 : std::vector<std::vector<Real>> * comps[3]
292 177635 : { &this->dphidxi_map, &this->dphideta_map, &this->dphidzeta_map };
293 133921 : FEInterface::all_shape_derivs(Dim, map_fe_type, elem, qp, comps, false);
294 : }
295 :
296 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
297 178613 : if (calculate_d2xyz)
298 85780 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
299 219867 : for (std::size_t p=0; p<n_qp; p++)
300 167742 : this->d2phidxi2_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 0, qp[p], false);
301 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
302 : }
303 :
304 227496 : break;
305 : }
306 : //------------------------------------------------------------
307 : // 2D
308 17719886 : case 2:
309 : {
310 : // Compute gradients of mapping shape functions i at quadrature points p
311 :
312 : // TODO: optimizations for the RATIONAL_BERNSTEIN+linear case
313 19582844 : if (is_linear)
314 : {
315 13140204 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
316 : {
317 9855153 : if (calculate_dxyz)
318 : {
319 9735225 : this->dphidxi_map[i][0] = shape_deriv_ptr (map_fe_type, elem, i, 0, qp[0], false);
320 9735225 : this->dphideta_map[i][0] = shape_deriv_ptr (map_fe_type, elem, i, 1, qp[0], false);
321 10697289 : for (std::size_t p=1; p<n_qp; p++)
322 : {
323 1033578 : this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
324 1033578 : this->dphideta_map[i][p] = this->dphideta_map[i][0];
325 : }
326 : }
327 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
328 9855153 : if (calculate_d2xyz)
329 : {
330 145017 : this->d2phidxi2_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 0, qp[0], false);
331 145017 : this->d2phidxideta_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 1, qp[0], false);
332 145017 : this->d2phideta2_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 2, qp[0], false);
333 152145 : for (std::size_t p=1; p<n_qp; p++)
334 : {
335 7776 : this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
336 7776 : this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
337 7776 : this->d2phideta2_map[i][p] = this->d2phideta2_map[i][0];
338 : }
339 : }
340 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
341 : }
342 : }
343 : else
344 : {
345 16297793 : if (calculate_dxyz)
346 : {
347 17555096 : std::vector<std::vector<Real>> * comps[3]
348 17555096 : { &this->dphidxi_map, &this->dphideta_map, &this->dphidzeta_map };
349 12479834 : FEInterface::all_shape_derivs(Dim, map_fe_type, elem, qp, comps, false);
350 : }
351 :
352 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
353 16297793 : if (calculate_d2xyz)
354 7621862 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
355 37886886 : for (std::size_t p=0; p<n_qp; p++)
356 : {
357 33600542 : this->d2phidxi2_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 0, qp[p], false);
358 33600542 : this->d2phidxideta_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 1, qp[p], false);
359 33600542 : this->d2phideta2_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 2, qp[p], false);
360 : }
361 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
362 : }
363 :
364 1862958 : break;
365 : }
366 :
367 : //------------------------------------------------------------
368 : // 3D
369 23401596 : case 3:
370 : {
371 : // Compute gradients of mapping shape functions i at quadrature points p
372 :
373 : // TODO: optimizations for the RATIONAL_BERNSTEIN+linear case
374 25523247 : if (is_linear)
375 : {
376 784570 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
377 : {
378 627656 : if (calculate_dxyz)
379 : {
380 599096 : this->dphidxi_map[i][0] = shape_deriv_ptr (map_fe_type, elem, i, 0, qp[0], false);
381 599096 : this->dphideta_map[i][0] = shape_deriv_ptr (map_fe_type, elem, i, 1, qp[0], false);
382 599096 : this->dphidzeta_map[i][0] = shape_deriv_ptr (map_fe_type, elem, i, 2, qp[0], false);
383 :
384 1409512 : for (std::size_t p=1; p<n_qp; p++)
385 : {
386 872968 : this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
387 872968 : this->dphideta_map[i][p] = this->dphideta_map[i][0];
388 872968 : this->dphidzeta_map[i][p] = this->dphidzeta_map[i][0];
389 : }
390 : }
391 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
392 627656 : if (calculate_d2xyz)
393 : {
394 482724 : this->d2phidxi2_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 0, qp[0], false);
395 482724 : this->d2phidxideta_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 1, qp[0], false);
396 482724 : this->d2phideta2_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 2, qp[0], false);
397 482724 : this->d2phidxidzeta_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 3, qp[0], false);
398 482724 : this->d2phidetadzeta_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 4, qp[0], false);
399 482724 : this->d2phidzeta2_map[i][0] = shape_second_deriv_ptr (map_fe_type, elem, i, 5, qp[0], false);
400 :
401 560340 : for (std::size_t p=1; p<n_qp; p++)
402 : {
403 79968 : this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
404 79968 : this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
405 79968 : this->d2phideta2_map[i][p] = this->d2phideta2_map[i][0];
406 79968 : this->d2phidxidzeta_map[i][p] = this->d2phidxidzeta_map[i][0];
407 79968 : this->d2phidetadzeta_map[i][p] = this->d2phidetadzeta_map[i][0];
408 79968 : this->d2phidzeta2_map[i][p] = this->d2phidzeta2_map[i][0];
409 : }
410 : }
411 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
412 : }
413 : }
414 : else
415 : {
416 25366333 : if (calculate_dxyz)
417 : {
418 30427295 : std::vector<std::vector<Real>> * comps[3]
419 30427295 : { &this->dphidxi_map, &this->dphideta_map, &this->dphidzeta_map };
420 22834363 : FEInterface::all_shape_derivs(Dim, map_fe_type, elem, qp, comps, false);
421 : }
422 :
423 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
424 25366333 : if (calculate_d2xyz)
425 90802106 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
426 253291639 : for (std::size_t p=0; p<n_qp; p++)
427 : {
428 181905157 : this->d2phidxi2_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 0, qp[p], false);
429 181905157 : this->d2phidxideta_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 1, qp[p], false);
430 181905157 : this->d2phideta2_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 2, qp[p], false);
431 181905157 : this->d2phidxidzeta_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 3, qp[p], false);
432 181905157 : this->d2phidetadzeta_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 4, qp[p], false);
433 181905157 : this->d2phidzeta2_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 5, qp[p], false);
434 : }
435 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
436 : }
437 :
438 2121651 : break;
439 : }
440 :
441 : default:
442 : libmesh_error_msg("Invalid Dim = " << Dim);
443 : }
444 48044778 : }
445 :
446 :
447 :
448 190806164 : void FEMap::compute_single_point_map(const unsigned int dim,
449 : const std::vector<Real> & qw,
450 : const Elem * elem,
451 : unsigned int p,
452 : const std::vector<const Node *> & elem_nodes,
453 : bool compute_second_derivatives)
454 : {
455 16018629 : libmesh_assert(elem);
456 16018629 : libmesh_assert(calculations_started);
457 : #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
458 : libmesh_assert(!compute_second_derivatives);
459 : #endif
460 :
461 : // this parameter is only accounted for if LIBMESH_DIM==2.
462 16018629 : libmesh_ignore(compute_second_derivatives);
463 :
464 190806164 : if (calculate_xyz)
465 3086044 : libmesh_assert_equal_to(phi_map.size(), elem_nodes.size());
466 :
467 : auto check_for_degenerate_map =
468 147521985 : [this, elem, p]
469 14831005 : (Real det_J) {
470 176928448 : if (det_J <= jacobian_tolerance)
471 : {
472 : // Don't call get_info() recursively if we're already
473 : // failing. get_info() calls Elem::volume() which may
474 : // call FE::reinit() and trigger the same failure again.
475 : static bool failing = false;
476 0 : if (!failing)
477 : {
478 0 : failing = true;
479 0 : std::string elem_info = elem->get_info();
480 0 : failing = false;
481 0 : if (calculate_xyz)
482 : {
483 0 : libmesh_degenerate_mapping_msg
484 : ("Jacobian " << jac[p] << " under tolerance " <<
485 : jacobian_tolerance << " at point " << xyz[p] <<
486 : " in element: " << elem_info);
487 : }
488 : else
489 : {
490 : // In this case xyz[p] is not defined, so don't
491 : // try to print it out.
492 0 : libmesh_degenerate_mapping_msg
493 : ("Jacobian " << jac[p] << " under tolerance " <<
494 : jacobian_tolerance << " at point index " << p <<
495 : " in element: " << elem_info);
496 : }
497 : }
498 : else
499 : {
500 : // We were already failing when we called this, so just
501 : // stop the current computation and return with
502 : // incomplete results.
503 0 : return;
504 : }
505 : }
506 190806164 : };
507 :
508 190806164 : switch (dim)
509 : {
510 : //--------------------------------------------------------------------
511 : // 0D
512 198780 : case 0:
513 : {
514 19049 : libmesh_assert(elem_nodes[0]);
515 198780 : if (calculate_xyz)
516 0 : xyz[p] = *elem_nodes[0];
517 198780 : if (calculate_dxyz)
518 : {
519 198780 : jac[p] = 1.0;
520 236878 : JxW[p] = qw[p];
521 : }
522 19049 : break;
523 : }
524 :
525 : //--------------------------------------------------------------------
526 : // 1D
527 44980629 : case 1:
528 : {
529 : // Clear the entities that will be summed
530 44980629 : if (calculate_xyz)
531 126662 : xyz[p].zero();
532 44980629 : if (calculate_dxyz)
533 42646625 : dxyzdxi_map[p].zero();
534 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
535 44980629 : if (calculate_d2xyz)
536 : {
537 42505990 : d2xyzdxi2_map[p].zero();
538 : // Inverse map second derivatives
539 45132000 : d2xidxyz2_map[p].assign(6, 0.);
540 : }
541 : #endif
542 :
543 : // compute x, dx, d2x at the quadrature point
544 135173839 : for (auto i : index_range(elem_nodes)) // sum over the nodes
545 : {
546 : // Reference to the point, helps eliminate
547 : // excessive temporaries in the inner loop
548 6134904 : libmesh_assert(elem_nodes[i]);
549 90193210 : const Point & elem_point = *elem_nodes[i];
550 :
551 90193210 : if (calculate_xyz)
552 440420 : xyz[p].add_scaled (elem_point, phi_map[i][p] );
553 90193210 : if (calculate_dxyz)
554 96046306 : dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p]);
555 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
556 90193210 : if (calculate_d2xyz)
557 95558237 : d2xyzdxi2_map[p].add_scaled(elem_point, d2phidxi2_map[i][p]);
558 : #endif
559 : }
560 :
561 : // Compute the jacobian
562 : //
563 : // 1D elements can live in 2D or 3D space.
564 : // The transformation matrix from local->global
565 : // coordinates is
566 : //
567 : // T = | dx/dxi |
568 : // | dy/dxi |
569 : // | dz/dxi |
570 : //
571 : // The generalized determinant of T (from the
572 : // so-called "normal" eqns.) is
573 : // jac = "det(T)" = sqrt(det(T'T))
574 : //
575 : // where T'= transpose of T, so
576 : //
577 : // jac = sqrt( (dx/dxi)^2 + (dy/dxi)^2 + (dz/dxi)^2 )
578 :
579 44980629 : if (calculate_dxyz)
580 : {
581 45284300 : jac[p] = dxyzdxi_map[p].norm();
582 :
583 42646625 : check_for_degenerate_map(jac[p]);
584 :
585 : // The inverse Jacobian entries also come from the
586 : // generalized inverse of T (see also the 2D element
587 : // living in 3D code).
588 45284300 : const Real jacm2 = 1./jac[p]/jac[p];
589 45284300 : dxidx_map[p] = jacm2*dxdxi_map(p);
590 : #if LIBMESH_DIM > 1
591 45284300 : dxidy_map[p] = jacm2*dydxi_map(p);
592 : #endif
593 : #if LIBMESH_DIM > 2
594 42646625 : dxidz_map[p] = jacm2*dzdxi_map(p);
595 : #endif
596 :
597 47921975 : JxW[p] = jac[p]*qw[p];
598 : }
599 :
600 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
601 :
602 44980629 : if (calculate_d2xyz)
603 : {
604 : #if LIBMESH_DIM == 1
605 : // Compute inverse map second derivatives for 1D-element-living-in-1D case
606 : this->compute_inverse_map_second_derivs(p);
607 : #elif LIBMESH_DIM == 2
608 : // Compute inverse map second derivatives for 1D-element-living-in-2D case
609 : // See JWP notes for details
610 :
611 : // numer = x_xi*x_{xi xi} + y_xi*y_{xi xi}
612 : Real numer =
613 : dxyzdxi_map[p](0)*d2xyzdxi2_map[p](0) +
614 : dxyzdxi_map[p](1)*d2xyzdxi2_map[p](1);
615 :
616 : // denom = (x_xi)^2 + (y_xi)^2 must be >= 0.0
617 : Real denom =
618 : dxyzdxi_map[p](0)*dxyzdxi_map[p](0) +
619 : dxyzdxi_map[p](1)*dxyzdxi_map[p](1);
620 :
621 : if (denom <= 0.0)
622 : {
623 : // Don't call print_info() recursively if we're already
624 : // failing. print_info() calls Elem::volume() which may
625 : // call FE::reinit() and trigger the same failure again.
626 : static bool failing = false;
627 : if (!failing)
628 : {
629 : failing = true;
630 : elem->print_info(libMesh::err);
631 : failing = false;
632 : libmesh_error_msg("Encountered invalid 1D element!");
633 : }
634 : else
635 : {
636 : // We were already failing when we called this, so just
637 : // stop the current computation and return with
638 : // incomplete results.
639 : return;
640 : }
641 : }
642 :
643 : // xi_{x x}
644 : d2xidxyz2_map[p][0] = -numer * dxidx_map[p]*dxidx_map[p] / denom;
645 :
646 : // xi_{x y}
647 : d2xidxyz2_map[p][1] = -numer * dxidx_map[p]*dxidy_map[p] / denom;
648 :
649 : // xi_{y y}
650 : d2xidxyz2_map[p][3] = -numer * dxidy_map[p]*dxidy_map[p] / denom;
651 :
652 : #elif LIBMESH_DIM == 3
653 : // Compute inverse map second derivatives for 1D-element-living-in-3D case
654 : // See JWP notes for details
655 :
656 : // numer = x_xi*x_{xi xi} + y_xi*y_{xi xi} + z_xi*z_{xi xi}
657 : Real numer =
658 45132000 : dxyzdxi_map[p](0)*d2xyzdxi2_map[p](0) +
659 42505990 : dxyzdxi_map[p](1)*d2xyzdxi2_map[p](1) +
660 42505990 : dxyzdxi_map[p](2)*d2xyzdxi2_map[p](2);
661 :
662 : // denom = (x_xi)^2 + (y_xi)^2 + (z_xi)^2 must be >= 0.0
663 : Real denom =
664 45132000 : dxyzdxi_map[p](0)*dxyzdxi_map[p](0) +
665 42505990 : dxyzdxi_map[p](1)*dxyzdxi_map[p](1) +
666 42505990 : dxyzdxi_map[p](2)*dxyzdxi_map[p](2);
667 :
668 42505990 : if (denom <= 0.0)
669 : {
670 : // Don't call print_info() recursively if we're already
671 : // failing. print_info() calls Elem::volume() which may
672 : // call FE::reinit() and trigger the same failure again.
673 : static bool failing = false;
674 0 : if (!failing)
675 : {
676 0 : failing = true;
677 0 : elem->print_info(libMesh::err);
678 0 : failing = false;
679 0 : libmesh_error_msg("Encountered invalid 1D element!");
680 : }
681 : else
682 : {
683 : // We were already failing when we called this, so just
684 : // stop the current computation and return with
685 : // incomplete results.
686 0 : return;
687 : }
688 : }
689 :
690 : // xi_{x x}
691 45132000 : d2xidxyz2_map[p][0] = -numer * dxidx_map[p]*dxidx_map[p] / denom;
692 :
693 : // xi_{x y}
694 42505990 : d2xidxyz2_map[p][1] = -numer * dxidx_map[p]*dxidy_map[p] / denom;
695 :
696 : // xi_{x z}
697 42505990 : d2xidxyz2_map[p][2] = -numer * dxidx_map[p]*dxidz_map[p] / denom;
698 :
699 : // xi_{y y}
700 42505990 : d2xidxyz2_map[p][3] = -numer * dxidy_map[p]*dxidy_map[p] / denom;
701 :
702 : // xi_{y z}
703 42505990 : d2xidxyz2_map[p][4] = -numer * dxidy_map[p]*dxidz_map[p] / denom;
704 :
705 : // xi_{z z}
706 42505990 : d2xidxyz2_map[p][5] = -numer * dxidz_map[p]*dxidz_map[p] / denom;
707 : #endif //LIBMESH_DIM == 3
708 : } // calculate_d2xyz
709 :
710 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
711 :
712 : // done computing the map
713 3057986 : break;
714 : }
715 :
716 :
717 : //--------------------------------------------------------------------
718 : // 2D
719 88076592 : case 2:
720 : {
721 : //------------------------------------------------------------------
722 : // Compute the (x,y) values at the quadrature points,
723 : // the Jacobian at the quadrature points
724 :
725 88076592 : if (calculate_xyz)
726 17313551 : xyz[p].zero();
727 :
728 88076592 : if (calculate_dxyz)
729 : {
730 80826915 : dxyzdxi_map[p].zero();
731 15016145 : dxyzdeta_map[p].zero();
732 : }
733 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
734 88076592 : if (calculate_d2xyz)
735 : {
736 1357749 : d2xyzdxi2_map[p].zero();
737 230462 : d2xyzdxideta_map[p].zero();
738 230462 : d2xyzdeta2_map[p].zero();
739 : // Inverse map second derivatives
740 1472980 : d2xidxyz2_map[p].assign(6, 0.);
741 1472980 : d2etadxyz2_map[p].assign(6, 0.);
742 : }
743 : #endif
744 :
745 :
746 : // compute (x,y) at the quadrature points, derivatives once
747 521884361 : for (auto i : index_range(elem_nodes)) // sum over the nodes
748 : {
749 : // Reference to the point, helps eliminate
750 : // excessive temporaries in the inner loop
751 39921520 : libmesh_assert(elem_nodes[i]);
752 433807769 : const Point & elem_point = *elem_nodes[i];
753 :
754 433807769 : if (calculate_xyz)
755 131931158 : xyz[p].add_scaled (elem_point, phi_map[i][p] );
756 :
757 433807769 : if (calculate_dxyz)
758 : {
759 444611849 : dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
760 444611849 : dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p]);
761 : }
762 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
763 433807769 : if (calculate_d2xyz)
764 : {
765 10968110 : d2xyzdxi2_map[p].add_scaled (elem_point, d2phidxi2_map[i][p]);
766 10968110 : d2xyzdxideta_map[p].add_scaled (elem_point, d2phidxideta_map[i][p]);
767 10968110 : d2xyzdeta2_map[p].add_scaled (elem_point, d2phideta2_map[i][p]);
768 : }
769 : #endif
770 : }
771 :
772 88076592 : if (calculate_dxyz)
773 : {
774 : // compute the jacobian once
775 7536834 : const Real dx_dxi = dxdxi_map(p),
776 7536834 : dx_deta = dxdeta_map(p),
777 7536834 : dy_dxi = dydxi_map(p),
778 7536834 : dy_deta = dydeta_map(p);
779 :
780 : #if LIBMESH_DIM == 2
781 : // Compute the Jacobian. This assumes the 2D face
782 : // lives in 2D space
783 : //
784 : // Symbolically, the matrix determinant is
785 : //
786 : // | dx/dxi dx/deta |
787 : // jac = | dy/dxi dy/deta |
788 : //
789 : // jac = dx/dxi*dy/deta - dx/deta*dy/dxi
790 : jac[p] = (dx_dxi*dy_deta - dx_deta*dy_dxi);
791 :
792 : check_for_degenerate_map(jac[p]);
793 :
794 : JxW[p] = jac[p]*qw[p];
795 :
796 : // Compute the shape function derivatives wrt x,y at the
797 : // quadrature points
798 : const Real inv_jac = 1./jac[p];
799 :
800 : dxidx_map[p] = dy_deta*inv_jac; //dxi/dx = (1/J)*dy/deta
801 : dxidy_map[p] = -dx_deta*inv_jac; //dxi/dy = -(1/J)*dx/deta
802 : detadx_map[p] = -dy_dxi* inv_jac; //deta/dx = -(1/J)*dy/dxi
803 : detady_map[p] = dx_dxi* inv_jac; //deta/dy = (1/J)*dx/dxi
804 :
805 : dxidz_map[p] = detadz_map[p] = 0.;
806 :
807 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
808 : if (compute_second_derivatives)
809 : this->compute_inverse_map_second_derivs(p);
810 : #endif
811 : #else // LIBMESH_DIM == 3
812 :
813 7536834 : const Real dz_dxi = dzdxi_map(p),
814 7536834 : dz_deta = dzdeta_map(p);
815 :
816 : // Compute the Jacobian. This assumes a 2D face in
817 : // 3D space.
818 : //
819 : // The transformation matrix T from local to global
820 : // coordinates is
821 : //
822 : // | dx/dxi dx/deta |
823 : // T = | dy/dxi dy/deta |
824 : // | dz/dxi dz/deta |
825 : // note det(T' T) = det(T')det(T) = det(T)det(T)
826 : // so det(T) = std::sqrt(det(T' T))
827 : //
828 : //----------------------------------------------
829 : // Notes:
830 : //
831 : // dX = R dXi -> R'dX = R'R dXi
832 : // (R^-1)dX = dXi [(R'R)^-1 R']dX = dXi
833 : //
834 : // so R^-1 = (R'R)^-1 R'
835 : //
836 : // and R^-1 R = (R'R)^-1 R'R = I.
837 : //
838 110859205 : const Real g11 = (dx_dxi*dx_dxi +
839 80826915 : dy_dxi*dy_dxi +
840 80826915 : dz_dxi*dz_dxi);
841 :
842 110859205 : const Real g12 = (dx_dxi*dx_deta +
843 80826915 : dy_dxi*dy_deta +
844 80826915 : dz_dxi*dz_deta);
845 :
846 7536834 : const Real g21 = g12;
847 :
848 110859205 : const Real g22 = (dx_deta*dx_deta +
849 80826915 : dy_deta*dy_deta +
850 80826915 : dz_deta*dz_deta);
851 :
852 80826915 : const Real det = (g11*g22 - g12*g21);
853 :
854 80826915 : check_for_degenerate_map(det);
855 :
856 80826915 : const Real inv_det = 1./det;
857 80826915 : jac[p] = std::sqrt(det);
858 :
859 88306226 : JxW[p] = jac[p]*qw[p];
860 :
861 80826915 : const Real g11inv = g22*inv_det;
862 80826915 : const Real g12inv = -g12*inv_det;
863 7536834 : const Real g21inv = -g21*inv_det;
864 80826915 : const Real g22inv = g11*inv_det;
865 :
866 80826915 : dxidx_map[p] = g11inv*dx_dxi + g12inv*dx_deta;
867 80826915 : dxidy_map[p] = g11inv*dy_dxi + g12inv*dy_deta;
868 80826915 : dxidz_map[p] = g11inv*dz_dxi + g12inv*dz_deta;
869 :
870 80826915 : detadx_map[p] = g21inv*dx_dxi + g22inv*dx_deta;
871 80826915 : detady_map[p] = g21inv*dy_dxi + g22inv*dy_deta;
872 80826915 : detadz_map[p] = g21inv*dz_dxi + g22inv*dz_deta;
873 :
874 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
875 :
876 80826915 : if (calculate_d2xyz)
877 : {
878 : // Compute inverse map second derivative values for
879 : // 2D-element-living-in-3D case. We pursue a least-squares
880 : // solution approach for this "non-square" case, see JWP notes
881 : // for details.
882 :
883 : // A = [ x_{xi xi} x_{eta eta} ]
884 : // [ y_{xi xi} y_{eta eta} ]
885 : // [ z_{xi xi} z_{eta eta} ]
886 1588211 : DenseMatrix<Real> A(3,2);
887 1588211 : A(0,0) = d2xyzdxi2_map[p](0); A(0,1) = d2xyzdeta2_map[p](0);
888 1472980 : A(1,0) = d2xyzdxi2_map[p](1); A(1,1) = d2xyzdeta2_map[p](1);
889 1472980 : A(2,0) = d2xyzdxi2_map[p](2); A(2,1) = d2xyzdeta2_map[p](2);
890 :
891 : // J^T, the transpose of the Jacobian matrix
892 1588211 : DenseMatrix<Real> JT(2,3);
893 1357749 : JT(0,0) = dx_dxi; JT(0,1) = dy_dxi; JT(0,2) = dz_dxi;
894 1472980 : JT(1,0) = dx_deta; JT(1,1) = dy_deta; JT(1,2) = dz_deta;
895 :
896 : // (J^T J)^(-1), this has already been computed for us above...
897 1588211 : DenseMatrix<Real> JTJinv(2,2);
898 1357749 : JTJinv(0,0) = g11inv; JTJinv(0,1) = g12inv;
899 1357749 : JTJinv(1,0) = g21inv; JTJinv(1,1) = g22inv;
900 :
901 : // Some helper variables
902 : RealVectorValue
903 460924 : dxi (dxidx_map[p], dxidy_map[p], dxidz_map[p]),
904 460924 : deta (detadx_map[p], detady_map[p], detadz_map[p]);
905 :
906 : // To be filled in below
907 1357749 : DenseVector<Real> tmp1(2);
908 1357749 : DenseVector<Real> tmp2(3);
909 1357749 : DenseVector<Real> tmp3(2);
910 :
911 : // For (s,t) in {(x,x), (x,y), (x,z), (y,y), (y,z), (z,z)}, compute the
912 : // vector of inverse map second derivatives [xi_{s t}, eta_{s t}]
913 115231 : unsigned ctr=0;
914 5430996 : for (unsigned s=0; s<3; ++s)
915 12219741 : for (unsigned t=s; t<3; ++t)
916 : {
917 : // Construct tmp1 = [xi_s*xi_t, eta_s*eta_t]
918 8146494 : tmp1(0) = dxi(s)*dxi(t);
919 8146494 : tmp1(1) = deta(s)*deta(t);
920 :
921 : // Compute tmp2 = A * tmp1
922 8146494 : A.vector_mult(tmp2, tmp1);
923 :
924 : // Compute scalar value "alpha"
925 8146494 : Real alpha = dxi(s)*deta(t) + deta(s)*dxi(t);
926 :
927 : // Compute tmp2 <- tmp2 + alpha * x_{xi eta}
928 32585976 : for (unsigned i=0; i<3; ++i)
929 28587798 : tmp2(i) += alpha*d2xyzdxideta_map[p](i);
930 :
931 : // Compute tmp3 = J^T * tmp2
932 8146494 : JT.vector_mult(tmp3, tmp2);
933 :
934 : // Compute tmp1 = (J^T J)^(-1) * tmp3. tmp1 is available for us to reuse.
935 8146494 : JTJinv.vector_mult(tmp1, tmp3);
936 :
937 : // Fill in appropriate entries, don't forget to multiply by -1!
938 8837880 : d2xidxyz2_map[p][ctr] = -tmp1(0);
939 8837880 : d2etadxyz2_map[p][ctr] = -tmp1(1);
940 :
941 : // Increment the counter
942 8146494 : ctr++;
943 : }
944 1127287 : }
945 :
946 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
947 :
948 : #endif // LIBMESH_DIM == 3
949 : }
950 : // done computing the map
951 8154262 : break;
952 : }
953 :
954 :
955 :
956 : //--------------------------------------------------------------------
957 : // 3D
958 57550163 : case 3:
959 : {
960 : //------------------------------------------------------------------
961 : // Compute the (x,y,z) values at the quadrature points,
962 : // the Jacobian at the quadrature point
963 :
964 : // Clear the entities that will be summed
965 57550163 : if (calculate_xyz)
966 16684379 : xyz[p].zero ();
967 57550163 : if (calculate_dxyz)
968 : {
969 53454908 : dxyzdxi_map[p].zero ();
970 8906224 : dxyzdeta_map[p].zero ();
971 8906224 : dxyzdzeta_map[p].zero ();
972 : }
973 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
974 57550163 : if (calculate_d2xyz)
975 : {
976 6135761 : d2xyzdxi2_map[p].zero();
977 1001996 : d2xyzdxideta_map[p].zero();
978 1001996 : d2xyzdxidzeta_map[p].zero();
979 1001996 : d2xyzdeta2_map[p].zero();
980 1001996 : d2xyzdetadzeta_map[p].zero();
981 1001996 : d2xyzdzeta2_map[p].zero();
982 : // Inverse map second derivatives
983 6636759 : d2xidxyz2_map[p].assign(6, 0.);
984 6636759 : d2etadxyz2_map[p].assign(6, 0.);
985 6636759 : d2zetadxyz2_map[p].assign(6, 0.);
986 : }
987 : #endif
988 :
989 :
990 : // compute (x,y,z) at the quadrature points,
991 : // dxdxi, dydxi, dzdxi,
992 : // dxdeta, dydeta, dzdeta,
993 : // dxdzeta, dydzeta, dzdzeta all once
994 750485475 : for (auto i : index_range(elem_nodes)) // sum over the nodes
995 : {
996 : // Reference to the point, helps eliminate
997 : // excessive temporaries in the inner loop
998 57604551 : libmesh_assert(elem_nodes[i]);
999 692935312 : const Point & elem_point = *elem_nodes[i];
1000 :
1001 692935312 : if (calculate_xyz)
1002 277460631 : xyz[p].add_scaled (elem_point, phi_map[i][p] );
1003 692935312 : if (calculate_dxyz)
1004 : {
1005 731964550 : dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
1006 731964550 : dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p] );
1007 731964550 : dxyzdzeta_map[p].add_scaled (elem_point, dphidzeta_map[i][p]);
1008 : }
1009 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1010 692935312 : if (calculate_d2xyz)
1011 : {
1012 100131912 : d2xyzdxi2_map[p].add_scaled (elem_point,
1013 23004636 : d2phidxi2_map[i][p]);
1014 100131912 : d2xyzdxideta_map[p].add_scaled (elem_point,
1015 23004636 : d2phidxideta_map[i][p]);
1016 100131912 : d2xyzdxidzeta_map[p].add_scaled (elem_point,
1017 23004636 : d2phidxidzeta_map[i][p]);
1018 100131912 : d2xyzdeta2_map[p].add_scaled (elem_point,
1019 23004636 : d2phideta2_map[i][p]);
1020 100131912 : d2xyzdetadzeta_map[p].add_scaled (elem_point,
1021 23004636 : d2phidetadzeta_map[i][p]);
1022 100131912 : d2xyzdzeta2_map[p].add_scaled (elem_point,
1023 23004636 : d2phidzeta2_map[i][p]);
1024 : }
1025 : #endif
1026 : }
1027 :
1028 57550163 : if (calculate_dxyz)
1029 : {
1030 : // compute the jacobian
1031 : const Real
1032 4447752 : dx_dxi = dxdxi_map(p), dy_dxi = dydxi_map(p), dz_dxi = dzdxi_map(p),
1033 4447752 : dx_deta = dxdeta_map(p), dy_deta = dydeta_map(p), dz_deta = dzdeta_map(p),
1034 4447752 : dx_dzeta = dxdzeta_map(p), dy_dzeta = dydzeta_map(p), dz_dzeta = dzdzeta_map(p);
1035 :
1036 : // Symbolically, the matrix determinant is
1037 : //
1038 : // | dx/dxi dy/dxi dz/dxi |
1039 : // jac = | dx/deta dy/deta dz/deta |
1040 : // | dx/dzeta dy/dzeta dz/dzeta |
1041 : //
1042 : // jac = dx/dxi*(dy/deta*dz/dzeta - dz/deta*dy/dzeta) +
1043 : // dy/dxi*(dz/deta*dx/dzeta - dx/deta*dz/dzeta) +
1044 : // dz/dxi*(dx/deta*dy/dzeta - dy/deta*dx/dzeta)
1045 :
1046 80173580 : jac[p] = (dx_dxi*(dy_deta*dz_dzeta - dz_deta*dy_dzeta) +
1047 62361132 : dy_dxi*(dz_deta*dx_dzeta - dx_deta*dz_dzeta) +
1048 53454908 : dz_dxi*(dx_deta*dy_dzeta - dy_deta*dx_dzeta));
1049 :
1050 53454908 : check_for_degenerate_map(jac[p]);
1051 :
1052 62371852 : JxW[p] = jac[p]*qw[p];
1053 :
1054 : // Compute the shape function derivatives wrt x,y at the
1055 : // quadrature points
1056 53454908 : const Real inv_jac = 1./jac[p];
1057 :
1058 53454908 : dxidx_map[p] = (dy_deta*dz_dzeta - dz_deta*dy_dzeta)*inv_jac;
1059 53454908 : dxidy_map[p] = (dz_deta*dx_dzeta - dx_deta*dz_dzeta)*inv_jac;
1060 53454908 : dxidz_map[p] = (dx_deta*dy_dzeta - dy_deta*dx_dzeta)*inv_jac;
1061 :
1062 53454908 : detadx_map[p] = (dz_dxi*dy_dzeta - dy_dxi*dz_dzeta )*inv_jac;
1063 53454908 : detady_map[p] = (dx_dxi*dz_dzeta - dz_dxi*dx_dzeta )*inv_jac;
1064 53454908 : detadz_map[p] = (dy_dxi*dx_dzeta - dx_dxi*dy_dzeta )*inv_jac;
1065 :
1066 53454908 : dzetadx_map[p] = (dy_dxi*dz_deta - dz_dxi*dy_deta )*inv_jac;
1067 53454908 : dzetady_map[p] = (dz_dxi*dx_deta - dx_dxi*dz_deta )*inv_jac;
1068 57913380 : dzetadz_map[p] = (dx_dxi*dy_deta - dy_dxi*dx_deta )*inv_jac;
1069 : }
1070 :
1071 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1072 57550163 : if (compute_second_derivatives)
1073 281008 : this->compute_inverse_map_second_derivs(p);
1074 : #endif
1075 : // done computing the map
1076 4787332 : break;
1077 : }
1078 :
1079 0 : default:
1080 0 : libmesh_error_msg("Invalid dim = " << dim);
1081 : }
1082 : }
1083 :
1084 :
1085 :
1086 164140341 : void FEMap::resize_quadrature_map_vectors(const unsigned int dim, unsigned int n_qp)
1087 : {
1088 : // We're calculating now!
1089 13804954 : this->determine_calculations();
1090 :
1091 : // Resize the vectors to hold data at the quadrature points
1092 164140341 : if (calculate_xyz)
1093 33748006 : xyz.resize(n_qp);
1094 164140341 : if (calculate_dxyz)
1095 : {
1096 150483677 : dxyzdxi_map.resize(n_qp);
1097 150483677 : dxidx_map.resize(n_qp);
1098 150483677 : dxidy_map.resize(n_qp); // 1D element may live in 2D ...
1099 150483677 : dxidz_map.resize(n_qp); // ... or 3D
1100 : }
1101 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1102 164140341 : if (calculate_d2xyz)
1103 : {
1104 49817844 : d2xyzdxi2_map.resize(n_qp);
1105 :
1106 : // Inverse map second derivatives
1107 49817844 : d2xidxyz2_map.resize(n_qp);
1108 245892721 : for (auto i : index_range(d2xidxyz2_map))
1109 208733091 : d2xidxyz2_map[i].assign(6, 0.);
1110 : }
1111 : #endif
1112 164140341 : if (dim > 1)
1113 : {
1114 118963074 : if (calculate_dxyz)
1115 : {
1116 107640414 : dxyzdeta_map.resize(n_qp);
1117 107640414 : detadx_map.resize(n_qp);
1118 107640414 : detady_map.resize(n_qp);
1119 107640414 : detadz_map.resize(n_qp);
1120 : }
1121 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1122 118963074 : if (calculate_d2xyz)
1123 : {
1124 7311854 : d2xyzdxideta_map.resize(n_qp);
1125 7311854 : d2xyzdeta2_map.resize(n_qp);
1126 :
1127 : // Inverse map second derivatives
1128 7311854 : d2etadxyz2_map.resize(n_qp);
1129 33413335 : for (auto i : index_range(d2etadxyz2_map))
1130 28259881 : d2etadxyz2_map[i].assign(6, 0.);
1131 : }
1132 : #endif
1133 118963074 : if (dim > 2)
1134 : {
1135 32733126 : if (calculate_dxyz)
1136 : {
1137 28660143 : dxyzdzeta_map.resize (n_qp);
1138 28660143 : dzetadx_map.resize (n_qp);
1139 28660143 : dzetady_map.resize (n_qp);
1140 28660143 : dzetadz_map.resize (n_qp);
1141 : }
1142 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1143 32733126 : if (calculate_d2xyz)
1144 : {
1145 6113849 : d2xyzdxidzeta_map.resize(n_qp);
1146 6113849 : d2xyzdetadzeta_map.resize(n_qp);
1147 6113849 : d2xyzdzeta2_map.resize(n_qp);
1148 :
1149 : // Inverse map second derivatives
1150 6113849 : d2zetadxyz2_map.resize(n_qp);
1151 23924943 : for (auto i : index_range(d2zetadxyz2_map))
1152 19276489 : d2zetadxyz2_map[i].assign(6, 0.);
1153 : }
1154 : #endif
1155 : }
1156 : }
1157 :
1158 164140341 : if (calculate_dxyz)
1159 : {
1160 150483677 : jac.resize(n_qp);
1161 150483677 : JxW.resize(n_qp);
1162 : }
1163 164140341 : }
1164 :
1165 :
1166 :
1167 160879314 : void FEMap::compute_affine_map(const unsigned int dim,
1168 : const std::vector<Real> & qw,
1169 : const Elem * elem)
1170 : {
1171 : // Start logging the map computation.
1172 27079078 : LOG_SCOPE("compute_affine_map()", "FEMap");
1173 :
1174 13539539 : libmesh_assert(elem);
1175 :
1176 26812091 : const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1177 :
1178 : // Resize the vectors to hold data at the quadrature points
1179 160879314 : this->resize_quadrature_map_vectors(dim, n_qp);
1180 :
1181 : // Determine the nodes contributing to element elem
1182 160879314 : unsigned int n_nodes = elem->n_nodes();
1183 160879314 : _elem_nodes.resize(elem->n_nodes());
1184 1076188219 : for (unsigned int i=0; i<n_nodes; i++)
1185 1071407841 : _elem_nodes[i] = elem->node_ptr(i);
1186 :
1187 : // Compute map at quadrature point 0
1188 160879314 : this->compute_single_point_map(dim, qw, elem, 0, _elem_nodes, /*compute_second_derivatives=*/false);
1189 :
1190 : // Compute xyz at all other quadrature points
1191 160879314 : if (calculate_xyz)
1192 182949447 : for (unsigned int p=1; p<n_qp; p++)
1193 : {
1194 149383435 : xyz[p].zero();
1195 1838254421 : for (auto i : index_range(phi_map)) // sum over the nodes
1196 1986187598 : xyz[p].add_scaled (*_elem_nodes[i], phi_map[i][p]);
1197 : }
1198 :
1199 : // Copy other map data from quadrature point 0
1200 160879314 : if (calculate_dxyz)
1201 631013579 : for (unsigned int p=1; p<n_qp; p++) // for each extra quadrature point
1202 : {
1203 483782757 : dxyzdxi_map[p] = dxyzdxi_map[0];
1204 483782757 : dxidx_map[p] = dxidx_map[0];
1205 483782757 : dxidy_map[p] = dxidy_map[0];
1206 483782757 : dxidz_map[p] = dxidz_map[0];
1207 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1208 : // The map should be affine, so second derivatives are zero
1209 483782757 : if (calculate_d2xyz)
1210 19458462 : d2xyzdxi2_map[p] = 0.;
1211 : #endif
1212 483782757 : if (dim > 1)
1213 : {
1214 356223559 : dxyzdeta_map[p] = dxyzdeta_map[0];
1215 356223559 : detadx_map[p] = detadx_map[0];
1216 356223559 : detady_map[p] = detady_map[0];
1217 356223559 : detadz_map[p] = detadz_map[0];
1218 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1219 356223559 : if (calculate_d2xyz)
1220 : {
1221 3084342 : d2xyzdxideta_map[p] = 0.;
1222 3084342 : d2xyzdeta2_map[p] = 0.;
1223 : }
1224 : #endif
1225 356223559 : if (dim > 2)
1226 : {
1227 99144165 : dxyzdzeta_map[p] = dxyzdzeta_map[0];
1228 99144165 : dzetadx_map[p] = dzetadx_map[0];
1229 99144165 : dzetady_map[p] = dzetady_map[0];
1230 99144165 : dzetadz_map[p] = dzetadz_map[0];
1231 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1232 99144165 : if (calculate_d2xyz)
1233 : {
1234 1928794 : d2xyzdxidzeta_map[p] = 0.;
1235 1928794 : d2xyzdetadzeta_map[p] = 0.;
1236 1928794 : d2xyzdzeta2_map[p] = 0.;
1237 : }
1238 : #endif
1239 : }
1240 : }
1241 483782757 : jac[p] = jac[0];
1242 563604769 : JxW[p] = JxW[0] / qw[0] * qw[p];
1243 : }
1244 160879314 : }
1245 :
1246 :
1247 :
1248 0 : void FEMap::compute_null_map(const unsigned int dim,
1249 : const std::vector<Real> & qw)
1250 : {
1251 : // Start logging the map computation.
1252 0 : LOG_SCOPE("compute_null_map()", "FEMap");
1253 :
1254 0 : const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1255 :
1256 : // Resize the vectors to hold data at the quadrature points
1257 0 : this->resize_quadrature_map_vectors(dim, n_qp);
1258 :
1259 : // Compute "fake" xyz
1260 0 : for (unsigned int p=1; p<n_qp; p++)
1261 : {
1262 0 : if (calculate_xyz)
1263 0 : xyz[p].zero();
1264 :
1265 0 : if (calculate_dxyz)
1266 : {
1267 0 : dxyzdxi_map[p] = 0;
1268 0 : dxidx_map[p] = 0;
1269 0 : dxidy_map[p] = 0;
1270 0 : dxidz_map[p] = 0;
1271 : }
1272 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1273 0 : if (calculate_d2xyz)
1274 : {
1275 0 : d2xyzdxi2_map[p] = 0;
1276 : }
1277 : #endif
1278 0 : if (dim > 1)
1279 : {
1280 0 : if (calculate_dxyz)
1281 : {
1282 0 : dxyzdeta_map[p] = 0;
1283 0 : detadx_map[p] = 0;
1284 0 : detady_map[p] = 0;
1285 0 : detadz_map[p] = 0;
1286 : }
1287 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1288 0 : if (calculate_d2xyz)
1289 : {
1290 0 : d2xyzdxideta_map[p] = 0.;
1291 0 : d2xyzdeta2_map[p] = 0.;
1292 : }
1293 : #endif
1294 0 : if (dim > 2)
1295 : {
1296 0 : if (calculate_dxyz)
1297 : {
1298 0 : dxyzdzeta_map[p] = 0;
1299 0 : dzetadx_map[p] = 0;
1300 0 : dzetady_map[p] = 0;
1301 0 : dzetadz_map[p] = 0;
1302 : }
1303 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1304 0 : if (calculate_d2xyz)
1305 : {
1306 0 : d2xyzdxidzeta_map[p] = 0;
1307 0 : d2xyzdetadzeta_map[p] = 0;
1308 0 : d2xyzdzeta2_map[p] = 0;
1309 : }
1310 : #endif
1311 : }
1312 : }
1313 0 : if (calculate_dxyz)
1314 : {
1315 0 : jac[p] = 1;
1316 0 : JxW[p] = qw[p];
1317 : }
1318 : }
1319 0 : }
1320 :
1321 :
1322 :
1323 164140341 : void FEMap::compute_map(const unsigned int dim,
1324 : const std::vector<Real> & qw,
1325 : const Elem * elem,
1326 : bool calculate_d2phi)
1327 : {
1328 164140341 : if (!elem)
1329 : {
1330 0 : compute_null_map(dim, qw);
1331 13539539 : return;
1332 : }
1333 :
1334 164140341 : if (elem->has_affine_map())
1335 : {
1336 160879314 : compute_affine_map(dim, qw, elem);
1337 160879314 : return;
1338 : }
1339 : #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
1340 : libmesh_assert(!calculate_d2phi);
1341 : #endif
1342 :
1343 : // Start logging the map computation.
1344 530830 : LOG_SCOPE("compute_map()", "FEMap");
1345 :
1346 265415 : libmesh_assert(elem);
1347 :
1348 532434 : const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1349 :
1350 : // Resize the vectors to hold data at the quadrature points
1351 3261027 : this->resize_quadrature_map_vectors(dim, n_qp);
1352 :
1353 : // Determine the nodes contributing to element elem
1354 3261027 : if (elem->type() == TRI3SUBDIVISION)
1355 : {
1356 : // Subdivision surface FE require the 1-ring around elem
1357 256 : libmesh_assert_equal_to (dim, 2);
1358 256 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
1359 3072 : MeshTools::Subdivision::find_one_ring(sd_elem, _elem_nodes);
1360 : }
1361 : else
1362 : {
1363 : // All other FE use only the nodes of elem itself
1364 3257955 : _elem_nodes.resize(elem->n_nodes(), nullptr);
1365 30021308 : for (auto i : elem->node_index_range())
1366 31101601 : _elem_nodes[i] = elem->node_ptr(i);
1367 : }
1368 :
1369 : // Compute map at all quadrature points
1370 33187877 : for (unsigned int p=0; p!=n_qp; p++)
1371 29926850 : this->compute_single_point_map(dim, qw, elem, p, _elem_nodes, calculate_d2phi);
1372 : }
1373 :
1374 :
1375 :
1376 0 : void FEMap::print_JxW(std::ostream & os) const
1377 : {
1378 0 : for (auto i : index_range(JxW))
1379 0 : os << " [" << i << "]: " << JxW[i] << std::endl;
1380 0 : }
1381 :
1382 :
1383 :
1384 0 : void FEMap::print_xyz(std::ostream & os) const
1385 : {
1386 0 : for (auto i : index_range(xyz))
1387 0 : os << " [" << i << "]: " << xyz[i];
1388 0 : }
1389 :
1390 :
1391 :
1392 281008 : void FEMap::compute_inverse_map_second_derivs(unsigned p)
1393 : {
1394 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1395 : // Only certain second derivatives are valid depending on the
1396 : // dimension...
1397 35742 : std::set<unsigned> valid_indices;
1398 :
1399 : // Construct J^{-1}, A, and B matrices (see JWP's notes for details)
1400 : // for cases in which the element dimension matches LIBMESH_DIM.
1401 : #if LIBMESH_DIM==1
1402 : RealTensor
1403 : Jinv(dxidx_map[p], 0., 0.,
1404 : 0., 0., 0.,
1405 : 0., 0., 0.),
1406 :
1407 : A(d2xyzdxi2_map[p](0), 0., 0.,
1408 : 0., 0., 0.,
1409 : 0., 0., 0.),
1410 :
1411 : B(0., 0., 0.,
1412 : 0., 0., 0.,
1413 : 0., 0., 0.);
1414 :
1415 : RealVectorValue
1416 : dxi (dxidx_map[p], 0., 0.),
1417 : deta (0., 0., 0.),
1418 : dzeta(0., 0., 0.);
1419 :
1420 : // In 1D, we have only the xx second derivative
1421 : valid_indices.insert(0);
1422 :
1423 : #elif LIBMESH_DIM==2
1424 : RealTensor
1425 : Jinv(dxidx_map[p], dxidy_map[p], 0.,
1426 : detadx_map[p], detady_map[p], 0.,
1427 : 0., 0., 0.),
1428 :
1429 : A(d2xyzdxi2_map[p](0), d2xyzdeta2_map[p](0), 0.,
1430 : d2xyzdxi2_map[p](1), d2xyzdeta2_map[p](1), 0.,
1431 : 0., 0., 0.),
1432 :
1433 : B(d2xyzdxideta_map[p](0), 0., 0.,
1434 : d2xyzdxideta_map[p](1), 0., 0.,
1435 : 0., 0., 0.);
1436 :
1437 : RealVectorValue
1438 : dxi (dxidx_map[p], dxidy_map[p], 0.),
1439 : deta (detadx_map[p], detady_map[p], 0.),
1440 : dzeta(0., 0., 0.);
1441 :
1442 : // In 2D, we have xx, xy, and yy second derivatives
1443 : const unsigned tmp[3] = {0,1,3};
1444 : valid_indices.insert(tmp, tmp+3);
1445 :
1446 : #elif LIBMESH_DIM==3
1447 : RealTensor
1448 71484 : Jinv(dxidx_map[p], dxidy_map[p], dxidz_map[p],
1449 71484 : detadx_map[p], detady_map[p], detadz_map[p],
1450 352492 : dzetadx_map[p], dzetady_map[p], dzetadz_map[p]),
1451 :
1452 17871 : A(d2xyzdxi2_map[p](0), d2xyzdeta2_map[p](0), d2xyzdzeta2_map[p](0),
1453 17871 : d2xyzdxi2_map[p](1), d2xyzdeta2_map[p](1), d2xyzdzeta2_map[p](1),
1454 107226 : d2xyzdxi2_map[p](2), d2xyzdeta2_map[p](2), d2xyzdzeta2_map[p](2)),
1455 :
1456 17871 : B(d2xyzdxideta_map[p](0), d2xyzdxidzeta_map[p](0), d2xyzdetadzeta_map[p](0),
1457 17871 : d2xyzdxideta_map[p](1), d2xyzdxidzeta_map[p](1), d2xyzdetadzeta_map[p](1),
1458 107226 : d2xyzdxideta_map[p](2), d2xyzdxidzeta_map[p](2), d2xyzdetadzeta_map[p](2));
1459 :
1460 : RealVectorValue
1461 17871 : dxi (dxidx_map[p], dxidy_map[p], dxidz_map[p]),
1462 17871 : deta (detadx_map[p], detady_map[p], detadz_map[p]),
1463 17871 : dzeta(dzetadx_map[p], dzetady_map[p], dzetadz_map[p]);
1464 :
1465 : // In 3D, we have xx, xy, xz, yy, yz, and zz second derivatives
1466 281008 : const unsigned tmp[6] = {0,1,2,3,4,5};
1467 17871 : valid_indices.insert(tmp, tmp+6);
1468 :
1469 : #endif
1470 :
1471 : // For (s,t) in {(x,x), (x,y), (x,z), (y,y), (y,z), (z,z)}, compute the
1472 : // vector of inverse map second derivatives [xi_{s t}, eta_{s t}, zeta_{s t}]
1473 263137 : unsigned ctr=0;
1474 1124032 : for (unsigned s=0; s<3; ++s)
1475 2529072 : for (unsigned t=s; t<3; ++t)
1476 : {
1477 107226 : if (valid_indices.count(ctr))
1478 : {
1479 : RealVectorValue
1480 1686048 : v1(dxi(s)*dxi(t),
1481 1686048 : deta(s)*deta(t),
1482 1793274 : dzeta(s)*dzeta(t)),
1483 :
1484 1686048 : v2(dxi(s)*deta(t) + deta(s)*dxi(t),
1485 1686048 : dxi(s)*dzeta(t) + dzeta(s)*dxi(t),
1486 1793274 : deta(s)*dzeta(t) + dzeta(s)*deta(t));
1487 :
1488 : // Compute the inverse map second derivatives
1489 1686048 : RealVectorValue v3 = -Jinv*(A*v1 + B*v2);
1490 :
1491 : // Store them in the appropriate locations in the class data structures
1492 1793274 : d2xidxyz2_map[p][ctr] = v3(0);
1493 :
1494 : if (LIBMESH_DIM > 1)
1495 1793274 : d2etadxyz2_map[p][ctr] = v3(1);
1496 :
1497 : if (LIBMESH_DIM > 2)
1498 1900500 : d2zetadxyz2_map[p][ctr] = v3(2);
1499 : }
1500 :
1501 : // Increment the counter
1502 1686048 : ctr++;
1503 : }
1504 : #else
1505 : // to avoid compiler warnings:
1506 : libmesh_ignore(p);
1507 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1508 281008 : }
1509 :
1510 :
1511 :
1512 106910770 : Point FEMap::inverse_map (const unsigned int dim,
1513 : const Elem * elem,
1514 : const Point & physical_point,
1515 : const Real tolerance,
1516 : const bool secure,
1517 : const bool
1518 : extra_checks
1519 : )
1520 : {
1521 7605987 : libmesh_assert(elem);
1522 7605987 : libmesh_assert_greater_equal (tolerance, 0.);
1523 :
1524 7605987 : libmesh_ignore(extra_checks);
1525 :
1526 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1527 :
1528 : // TODO: possibly use the extra_checks parameter in InfFEMap::inverse_map() as well.
1529 :
1530 20908027 : if (elem->infinite())
1531 716 : return InfFEMap::inverse_map(dim, elem, physical_point, tolerance,
1532 711 : secure);
1533 : #endif
1534 :
1535 : // Start logging the map inversion.
1536 15211262 : LOG_SCOPE("inverse_map()", "FEMap");
1537 :
1538 : // How much did the point on the reference
1539 : // element change by in this Newton step?
1540 7605631 : Real inverse_map_error = 0.;
1541 :
1542 : // The point on the reference element. This is
1543 : // the "initial guess" for Newton's method. The
1544 : // centroid seems like a good idea, but computing
1545 : // it is a little more intensive than, say taking
1546 : // the zero point.
1547 : //
1548 : // Convergence should be insensitive of this choice
1549 : // for "good" elements.
1550 7605631 : Point p; // the zero point. No computation required
1551 :
1552 : // The number of iterations in the map inversion process.
1553 7605631 : unsigned int cnt = 0;
1554 :
1555 : // The number of iterations after which we give up and declare
1556 : // divergence
1557 7605631 : const unsigned int max_cnt = 10;
1558 :
1559 : // The distance (in master element space) beyond which we give up
1560 : // and declare divergence. This is no longer used...
1561 : // Real max_step_length = 4.;
1562 :
1563 :
1564 :
1565 : // Newton iteration loop.
1566 7395486 : do
1567 : {
1568 : // Where our current iterate \p p maps to.
1569 212470572 : const Point physical_guess = map(dim, elem, p);
1570 :
1571 : // How far our current iterate is from the actual point.
1572 15001117 : const Point delta = physical_point - physical_guess;
1573 :
1574 : // Increment in current iterate \p p, will be computed.
1575 15001117 : Point dp;
1576 :
1577 :
1578 : // The form of the map and how we invert it depends
1579 : // on the dimension that we are in.
1580 212470572 : switch (dim)
1581 : {
1582 : // ------------------------------------------------------------------
1583 : // 0D map inversion is trivial
1584 1070 : case 0:
1585 : {
1586 1070 : break;
1587 : }
1588 :
1589 : // ------------------------------------------------------------------
1590 : // 1D map inversion
1591 : //
1592 : // Here we find the point on a 1D reference element that maps to
1593 : // the point \p physical_point in the domain. This is a bit tricky
1594 : // since we do not want to assume that the point \p physical_point
1595 : // is also in a 1D domain. In particular, this method might get
1596 : // called on the edge of a 3D element, in which case
1597 : // \p physical_point actually lives in 3D.
1598 8694003 : case 1:
1599 : {
1600 8694003 : const Point dxi = map_deriv (dim, elem, 0, p);
1601 :
1602 : // Newton's method in this case looks like
1603 : //
1604 : // {X} - {X_n} = [J]*dp
1605 : //
1606 : // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x1 matrix
1607 : // d(x,y,z)/dxi, and we seek dp, a scalar. Since the above
1608 : // system is either overdetermined or rank-deficient, we will
1609 : // solve the normal equations for this system
1610 : //
1611 : // [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
1612 : //
1613 : // which involves the trivial inversion of the scalar
1614 : // G = [J]^T [J]
1615 888354 : const Real G = dxi*dxi;
1616 :
1617 8694003 : if (secure && G <= 0)
1618 0 : libmesh_degenerate_mapping_msg
1619 : ("inverse_map found a singular Jacobian " <<
1620 : " at master point " << p << " in element " <<
1621 : elem->id());
1622 :
1623 8694003 : const Real Ginv = 1./G;
1624 :
1625 888354 : const Real dxidelta = dxi*delta;
1626 :
1627 8694003 : dp(0) = Ginv*dxidelta;
1628 :
1629 : // No master elements have radius > 4, but sometimes we
1630 : // can take a step that big while still converging
1631 : // if (secure)
1632 : // libmesh_assert_less (dp.size(), max_step_length);
1633 :
1634 888354 : break;
1635 : }
1636 :
1637 :
1638 :
1639 : // ------------------------------------------------------------------
1640 : // 2D map inversion
1641 : //
1642 : // Here we find the point on a 2D reference element that maps to
1643 : // the point \p physical_point in the domain. This is a bit tricky
1644 : // since we do not want to assume that the point \p physical_point
1645 : // is also in a 2D domain. In particular, this method might get
1646 : // called on the face of a 3D element, in which case
1647 : // \p physical_point actually lives in 3D.
1648 99921643 : case 2:
1649 : {
1650 99921643 : const Point dxi = map_deriv (dim, elem, 0, p);
1651 99921643 : const Point deta = map_deriv (dim, elem, 1, p);
1652 :
1653 : // Newton's method in this case looks like
1654 : //
1655 : // {X} - {X_n} = [J]*{dp}
1656 : //
1657 : // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x2 matrix
1658 : // d(x,y,z)/d(xi,eta), and we seek {dp}, a 2x1 vector. Since
1659 : // the above system is either over-determined or rank-deficient,
1660 : // we will solve the normal equations for this system
1661 : //
1662 : // [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
1663 : //
1664 : // which involves the inversion of the 2x2 matrix
1665 : // [G] = [J]^T [J]
1666 : const Real
1667 8449216 : G11 = dxi*dxi, G12 = dxi*deta,
1668 8449216 : G21 = dxi*deta, G22 = deta*deta;
1669 :
1670 :
1671 99921643 : const Real det = (G11*G22 - G12*G21);
1672 :
1673 99921643 : if (secure && det == 0)
1674 0 : libmesh_degenerate_mapping_msg
1675 : ("inverse_map found a singular Jacobian " <<
1676 : " at master point " << p << " in element " <<
1677 : elem->id());
1678 :
1679 99921643 : const Real inv_det = 1./det;
1680 :
1681 : const Real
1682 99921643 : Ginv11 = G22*inv_det,
1683 99921643 : Ginv12 = -G12*inv_det,
1684 :
1685 8449216 : Ginv21 = -G21*inv_det,
1686 99921643 : Ginv22 = G11*inv_det;
1687 :
1688 :
1689 8449216 : const Real dxidelta = dxi*delta;
1690 8449216 : const Real detadelta = deta*delta;
1691 :
1692 99921643 : dp(0) = (Ginv11*dxidelta + Ginv12*detadelta);
1693 99921643 : dp(1) = (Ginv21*dxidelta + Ginv22*detadelta);
1694 :
1695 : // No master elements have radius > 4, but sometimes we
1696 : // can take a step that big while still converging
1697 : // if (secure)
1698 : // libmesh_assert_less (dp.size(), max_step_length);
1699 :
1700 8449216 : break;
1701 : }
1702 :
1703 :
1704 :
1705 : // ------------------------------------------------------------------
1706 : // 3D map inversion
1707 : //
1708 : // Here we find the point in a 3D reference element that maps to
1709 : // the point \p physical_point in a 3D domain. Nothing special
1710 : // has to happen here, since (unless the map is singular because
1711 : // you have a BAD element) the map will be invertible and we can
1712 : // apply Newton's method directly.
1713 103837444 : case 3:
1714 : {
1715 103837444 : const Point dxi = map_deriv (dim, elem, 0, p);
1716 103837444 : const Point deta = map_deriv (dim, elem, 1, p);
1717 103837444 : const Point dzeta = map_deriv (dim, elem, 2, p);
1718 :
1719 : // Newton's method in this case looks like
1720 : //
1721 : // {X} = {X_n} + [J]*{dp}
1722 : //
1723 : // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x3 matrix
1724 : // d(x,y,z)/d(xi,eta,zeta), and we seek {dp}, a 3x1 vector.
1725 : // Since the above system is nonsingular for invertible maps
1726 : // we will solve
1727 : //
1728 : // {dp} = [J]^-1 ({X} - {X_n})
1729 : //
1730 : // which involves the inversion of the 3x3 matrix [J]
1731 : libmesh_try
1732 : {
1733 103866897 : RealTensorValue(dxi(0), deta(0), dzeta(0),
1734 : dxi(1), deta(1), dzeta(1),
1735 103837444 : dxi(2), deta(2), dzeta(2)).solve(delta, dp);
1736 : }
1737 198138 : libmesh_catch (ConvergenceFailure &)
1738 : {
1739 : // If we encountered a singular Jacobian, but are at
1740 : // a singular node, return said singular node.
1741 : const unsigned int local_singular_node =
1742 187738 : elem->local_singular_node(physical_point, tolerance);
1743 187738 : if (local_singular_node != invalid_uint)
1744 : {
1745 240 : libmesh_assert_less(local_singular_node, elem->n_nodes());
1746 5136 : return elem->master_point(local_singular_node);
1747 : }
1748 :
1749 : // We encountered a singular Jacobian. The value of
1750 : // dp is zero, since it was never changed during the
1751 : // call to RealTensorValue::solve(). We don't want to
1752 : // continue iterating until max_cnt since there is no
1753 : // update to the Newton iterate, and we don't want to
1754 : // print the inverse_map_error value since it will
1755 : // confusingly be 0. Therefore, in the secure case we
1756 : // need to throw an error message while in the !secure
1757 : // case we can just return a far away point.
1758 182602 : if (secure)
1759 : {
1760 0 : libmesh_degenerate_mapping_msg(
1761 : "inverse_map found a singular Jacobian" <<
1762 : " at master point " << p << " in element " <<
1763 : elem->id());
1764 : }
1765 : else
1766 : {
1767 730408 : for (unsigned int i=0; i != dim; ++i)
1768 547806 : p(i) = 1e6;
1769 182602 : return p;
1770 : }
1771 177338 : }
1772 :
1773 : // No master elements have radius > 4, but sometimes we
1774 : // can take a step that big while still converging
1775 : // if (secure)
1776 : // libmesh_assert_less (dp.size(), max_step_length);
1777 :
1778 103649706 : break;
1779 : }
1780 :
1781 :
1782 : // Some other dimension?
1783 0 : default:
1784 0 : libmesh_error_msg("Invalid dim = " << dim);
1785 : } // end switch(Dim), dp now computed
1786 :
1787 :
1788 :
1789 : // ||P_n+1 - P_n||
1790 197438972 : inverse_map_error = dp.norm();
1791 :
1792 : // P_n+1 = P_n + dp
1793 14995917 : p.add (dp);
1794 :
1795 : // Increment the iteration count.
1796 212282834 : cnt++;
1797 :
1798 : // Watch for divergence of Newton's
1799 : // method. Here's how it goes:
1800 : // (1) For good elements, we expect convergence in 10
1801 : // iterations, with no too-large steps.
1802 : // - If called with (secure == true) and we have not yet converged
1803 : // print out a warning message.
1804 : // - If called with (secure == true) and we have not converged in
1805 : // 20 iterations abort
1806 : // (2) This method may be called in cases when the target point is not
1807 : // inside the element and we have no business expecting convergence.
1808 : // For these cases if we have not converged in 10 iterations forget
1809 : // about it.
1810 212282834 : if (cnt > max_cnt)
1811 : {
1812 : // Warn about divergence when secure is true - this
1813 : // shouldn't happen
1814 77576 : if (secure)
1815 : {
1816 : // Print every time in devel/dbg modes
1817 : #ifndef NDEBUG
1818 0 : libmesh_here();
1819 0 : libMesh::err << "WARNING: Newton scheme has not converged in "
1820 0 : << cnt << " iterations:" << std::endl
1821 0 : << " physical_point="
1822 0 : << physical_point
1823 0 : << " physical_guess="
1824 0 : << physical_guess
1825 0 : << " dp="
1826 0 : << dp
1827 0 : << " p="
1828 0 : << p
1829 0 : << " error=" << inverse_map_error
1830 0 : << " in element " << elem->id()
1831 0 : << std::endl;
1832 :
1833 0 : elem->print_info(libMesh::err);
1834 : #else
1835 : // In optimized mode, just print once that an inverse_map() call
1836 : // had trouble converging its Newton iteration.
1837 0 : libmesh_do_once(libMesh::err << "WARNING: At least one element took more than "
1838 : << max_cnt
1839 : << " iterations to converge in inverse_map()...\n"
1840 : << "Rerun in devel/dbg mode for more details."
1841 : << std::endl;);
1842 :
1843 : #endif // NDEBUG
1844 :
1845 0 : if (cnt > 2*max_cnt)
1846 : {
1847 : // Whether or not this is a degenerate map in the
1848 : // "Jacobian becomes singular or negative" sense,
1849 : // it's at least degenerate in the "straightforward
1850 : // Newton is failing to invert it" sense.
1851 0 : libmesh_degenerate_mapping_msg(
1852 : "inverse_map Newton FAILED to converge in " <<
1853 : cnt << " iterations in element " << elem->id() <<
1854 : " for physical point = " << physical_point <<
1855 : std::endl << elem->get_info());
1856 : }
1857 : }
1858 : // Return a far off point when secure is false - this
1859 : // should only happen when we're trying to map a point
1860 : // that's outside the element
1861 : else
1862 : {
1863 308737 : for (unsigned int i=0; i != dim; ++i)
1864 231161 : p(i) = 1e6;
1865 :
1866 77576 : return p;
1867 : }
1868 : }
1869 : }
1870 212205258 : while (inverse_map_error > tolerance);
1871 :
1872 :
1873 :
1874 : // If we are in debug mode and the user requested it, do two extra sanity checks.
1875 : #ifdef DEBUG
1876 :
1877 7598726 : if (extra_checks)
1878 : {
1879 : // Make sure the point \p p on the reference element actually
1880 : // does map to the point \p physical_point within a tolerance.
1881 :
1882 3598693 : const Point check = map (dim, elem, p);
1883 3598693 : const Point diff = physical_point - check;
1884 :
1885 3598693 : if (diff.norm() > tolerance)
1886 : {
1887 0 : libmesh_here();
1888 0 : libMesh::err << "WARNING: diff is "
1889 0 : << diff.norm()
1890 0 : << std::endl
1891 0 : << " point="
1892 0 : << physical_point;
1893 0 : libMesh::err << " local=" << check;
1894 0 : libMesh::err << " lref= " << p;
1895 :
1896 0 : elem->print_info(libMesh::err);
1897 : }
1898 :
1899 : // Make sure the point \p p on the reference element actually
1900 : // is
1901 :
1902 3598693 : if (!elem->on_reference_element(p, 2*tolerance))
1903 : {
1904 0 : libmesh_here();
1905 0 : libMesh::err << "WARNING: inverse_map of physical point "
1906 0 : << physical_point
1907 0 : << " is not on element." << '\n';
1908 0 : elem->print_info(libMesh::err);
1909 : }
1910 : }
1911 :
1912 : #endif
1913 :
1914 106644029 : return p;
1915 : }
1916 :
1917 :
1918 :
1919 2081723 : void FEMap::inverse_map (const unsigned int dim,
1920 : const Elem * elem,
1921 : const std::vector<Point> & physical_points,
1922 : std::vector<Point> & reference_points,
1923 : const Real tolerance,
1924 : const bool secure,
1925 : const bool extra_checks)
1926 : {
1927 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1928 699747 : if (elem->infinite())
1929 : {
1930 : // TODO: possibly use the extra_checks parameter in InfFEMap::inverse_map() as well.
1931 0 : libmesh_ignore(extra_checks);
1932 :
1933 0 : return InfFEMap::inverse_map(dim, elem, physical_points, reference_points, tolerance, secure);
1934 : }
1935 : #endif
1936 :
1937 : // The number of points to find the
1938 : // inverse map of
1939 468523 : const std::size_t n_points = physical_points.size();
1940 :
1941 : // Resize the vector to hold the points
1942 : // on the reference element
1943 2081723 : reference_points.resize(n_points);
1944 :
1945 : // Find the coordinates on the reference
1946 : // element of each point in physical space
1947 5213520 : for (std::size_t p=0; p<n_points; p++)
1948 3465341 : reference_points[p] =
1949 3465341 : inverse_map (dim, elem, physical_points[p], tolerance, secure, extra_checks);
1950 1381976 : }
1951 :
1952 :
1953 :
1954 242819514 : Point FEMap::map (const unsigned int dim,
1955 : const Elem * elem,
1956 : const Point & reference_point)
1957 : {
1958 20711328 : libmesh_assert(elem);
1959 20711328 : libmesh_ignore(dim);
1960 :
1961 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1962 50319717 : if (elem->infinite())
1963 80 : return InfFEMap::map(dim, elem, reference_point);
1964 : #endif
1965 :
1966 20711296 : Point p;
1967 :
1968 242819434 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
1969 242819434 : const FEType fe_type (elem->default_order(), mapping_family);
1970 :
1971 : // Do not consider the Elem::p_level(), if any, when computing the
1972 : // number of shape functions.
1973 242819434 : const unsigned int n_sf = FEInterface::n_shape_functions(fe_type, /*extra_order=*/0, elem);
1974 :
1975 : FEInterface::shape_ptr shape_ptr =
1976 242819434 : FEInterface::shape_function(fe_type, elem);
1977 :
1978 : // Lagrange basis functions are used for mapping
1979 2718446817 : for (unsigned int i=0; i<n_sf; i++)
1980 2475627383 : p.add_scaled (elem->point(i),
1981 2659297986 : shape_ptr(fe_type, elem, i, reference_point, false));
1982 :
1983 242819434 : return p;
1984 : }
1985 :
1986 :
1987 :
1988 565117454 : Point FEMap::map_deriv (const unsigned int dim,
1989 : const Elem * elem,
1990 : const unsigned int j,
1991 : const Point & reference_point)
1992 : {
1993 38893870 : libmesh_assert(elem);
1994 38893870 : libmesh_ignore(dim);
1995 :
1996 108236792 : if (elem->infinite())
1997 0 : libmesh_not_implemented();
1998 :
1999 38893870 : Point p;
2000 :
2001 565117454 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
2002 565117454 : const FEType fe_type (elem->default_order(), mapping_family);
2003 :
2004 : // Do not consider the Elem::p_level(), if any, when computing the
2005 : // number of shape functions.
2006 : const unsigned int n_sf =
2007 565117454 : FEInterface::n_shape_functions(fe_type, /*total_order=*/0, elem);
2008 :
2009 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
2010 565117454 : FEInterface::shape_deriv_function(fe_type, elem);
2011 :
2012 : // Lagrange basis functions are used for mapping
2013 6928365727 : for (unsigned int i=0; i<n_sf; i++)
2014 6363248273 : p.add_scaled (elem->point(i),
2015 6757585831 : shape_deriv_ptr(fe_type, elem, i, j, reference_point,
2016 : /*add_p_level=*/false));
2017 :
2018 604011324 : return p;
2019 : }
2020 :
2021 :
2022 :
2023 : // Explicit instantiation of FEMap member functions
2024 : template LIBMESH_EXPORT void FEMap::init_reference_to_physical_map<0>( const std::vector<Point> &, const Elem *);
2025 : template LIBMESH_EXPORT void FEMap::init_reference_to_physical_map<1>( const std::vector<Point> &, const Elem *);
2026 : template LIBMESH_EXPORT void FEMap::init_reference_to_physical_map<2>( const std::vector<Point> &, const Elem *);
2027 : template LIBMESH_EXPORT void FEMap::init_reference_to_physical_map<3>( const std::vector<Point> &, const Elem *);
2028 :
2029 : // subdivision elements are implemented only for 2D meshes & reimplement
2030 : // the inverse_maps method separately
2031 : INSTANTIATE_SUBDIVISION_MAPS;
2032 :
2033 : } // namespace libMesh
|
