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 1206664026 : FEMap::map_fe_type(const Elem & elem)
47 : {
48 1206664026 : switch (elem.mapping_type())
49 : {
50 631683 : case RATIONAL_BERNSTEIN_MAP:
51 631683 : return RATIONAL_BERNSTEIN;
52 1199684891 : case LAGRANGE_MAP:
53 1199684891 : 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 15405430 : FEMap::FEMap(Real jtol) :
64 13653619 : calculations_started(false),
65 13653619 : calculate_xyz(false),
66 13653619 : calculate_dxyz(false),
67 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
68 13653619 : calculate_d2xyz(false),
69 : #endif
70 16281039 : jacobian_tolerance(jtol)
71 15405430 : {}
72 :
73 :
74 :
75 15326404 : std::unique_ptr<FEMap> FEMap::build( FEType fe_type )
76 : {
77 15326404 : switch( fe_type.family )
78 : {
79 859376 : case XYZ:
80 859376 : return std::make_unique<FEXYZMap>();
81 :
82 14467028 : default:
83 14467028 : return std::make_unique<FEMap>();
84 : }
85 : }
86 :
87 :
88 :
89 24798064 : void FEMap::add_calculations()
90 : {
91 24798064 : this->calculations_started = false;
92 22556692 : this->phi_map.clear();
93 22556692 : this->dphidxi_map.clear();
94 22556692 : this->dphideta_map.clear();
95 22556692 : this->dphidzeta_map.clear();
96 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
97 22556692 : this->d2phidxi2_map.clear();
98 22556692 : this->d2phidxideta_map.clear();
99 22556692 : this->d2phideta2_map.clear();
100 22556692 : this->d2phidxidzeta_map.clear();
101 22556692 : this->d2phidetadzeta_map.clear();
102 22556692 : this->d2phidzeta2_map.clear();
103 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
104 24798064 : }
105 :
106 :
107 :
108 : template<unsigned int Dim>
109 47779483 : 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 8437784 : LOG_SCOPE("init_reference_to_physical_map()", "FEMap");
114 :
115 : // We're calculating now!
116 4218892 : this->determine_calculations();
117 :
118 : // The number of quadrature points.
119 8437248 : const std::size_t n_qp = qp.size();
120 :
121 : // The element type and order to use in
122 : // the map
123 47779483 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
124 47779483 : 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 47779483 : 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 4218892 : unsigned int old_n_qp = 0;
137 :
138 : do
139 : {
140 47779483 : if (calculate_xyz)
141 : {
142 30665793 : if (this->phi_map.size() == n_mapping_shape_functions)
143 : {
144 2976548 : old_n_qp = n_mapping_shape_functions ? this->phi_map[0].size() : 0;
145 264076 : break;
146 : }
147 25156623 : this->phi_map.resize (n_mapping_shape_functions);
148 : }
149 : if (Dim > 0)
150 : {
151 44604155 : if (calculate_dxyz)
152 : {
153 39411513 : if (this->dphidxi_map.size() == n_mapping_shape_functions)
154 : {
155 7256450 : old_n_qp = n_mapping_shape_functions ? this->dphidxi_map[0].size() : 0;
156 659821 : break;
157 : }
158 28937440 : 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 37347705 : if (calculate_d2xyz)
165 1899923 : this->d2phidxi2_map.resize (n_mapping_shape_functions);
166 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
167 : }
168 :
169 : if (Dim > 1)
170 : {
171 34656345 : if (calculate_dxyz)
172 28566466 : this->dphideta_map.resize (n_mapping_shape_functions);
173 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
174 34656345 : 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 18654742 : if (calculate_dxyz)
185 16343234 : this->dphidzeta_map.resize (n_mapping_shape_functions);
186 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
187 18654742 : 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 47779483 : if (old_n_qp != n_qp)
199 391211307 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
200 : {
201 353637614 : if (calculate_xyz)
202 272438434 : this->phi_map[i].resize (n_qp);
203 : if (Dim > 0)
204 : {
205 353438834 : if (calculate_dxyz)
206 307277821 : this->dphidxi_map[i].resize (n_qp);
207 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
208 353438834 : if (calculate_d2xyz)
209 : {
210 17234788 : 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 347916206 : if (Dim > 1 && calculate_dxyz)
226 306407461 : this->dphideta_map[i].resize (n_qp);
227 :
228 253053167 : if (Dim > 2 && calculate_dxyz)
229 234290289 : this->dphidzeta_map[i].resize (n_qp);
230 : }
231 : }
232 :
233 : // Optimize for the *linear* geometric elements case:
234 47779483 : bool is_linear = elem->is_linear();
235 :
236 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
237 47779483 : 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 47779483 : FEInterface::shape_second_deriv_function(map_fe_type, elem);
242 : #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
243 :
244 47779483 : if (calculate_xyz)
245 28133171 : 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 2512095 : 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 2739798 : if (is_linear)
265 : {
266 7682364 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
267 : {
268 5121576 : if (calculate_dxyz)
269 : {
270 582112 : this->dphidxi_map[i][0] =
271 570482 : shape_deriv_ptr(map_fe_type, elem, i, 0, qp[0], false);
272 2235820 : for (std::size_t p=1; p<n_qp; p++)
273 1696412 : this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
274 : }
275 :
276 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
277 5121576 : 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 179010 : if (calculate_dxyz)
290 : {
291 178261 : std::vector<std::vector<Real>> * comps[3]
292 178261 : { &this->dphidxi_map, &this->dphideta_map, &this->dphidzeta_map };
293 134171 : FEInterface::all_shape_derivs(Dim, map_fe_type, elem, qp, comps, false);
294 : }
295 :
296 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
297 179010 : if (calculate_d2xyz)
298 85852 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
299 220269 : for (std::size_t p=0; p<n_qp; p++)
300 168114 : 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 227703 : break;
305 : }
306 : //------------------------------------------------------------
307 : // 2D
308 17691052 : 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 19552556 : if (is_linear)
314 : {
315 13086456 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
316 : {
317 9814842 : if (calculate_dxyz)
318 : {
319 9694914 : this->dphidxi_map[i][0] = shape_deriv_ptr (map_fe_type, elem, i, 0, qp[0], false);
320 9694914 : this->dphideta_map[i][0] = shape_deriv_ptr (map_fe_type, elem, i, 1, qp[0], false);
321 10536045 : for (std::size_t p=1; p<n_qp; p++)
322 : {
323 908055 : this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
324 908055 : this->dphideta_map[i][p] = this->dphideta_map[i][0];
325 : }
326 : }
327 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
328 9814842 : 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 16280942 : if (calculate_dxyz)
346 : {
347 17535265 : std::vector<std::vector<Real>> * comps[3]
348 17535265 : { &this->dphidxi_map, &this->dphideta_map, &this->dphidzeta_map };
349 12463369 : FEInterface::all_shape_derivs(Dim, map_fe_type, elem, qp, comps, false);
350 : }
351 :
352 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
353 16280942 : if (calculate_d2xyz)
354 7621841 : for (unsigned int i=0; i<n_mapping_shape_functions; i++)
355 37886220 : for (std::size_t p=0; p<n_qp; p++)
356 : {
357 33599894 : this->d2phidxi2_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 0, qp[p], false);
358 33599894 : this->d2phidxideta_map[i][p] = shape_second_deriv_ptr (map_fe_type, elem, i, 1, qp[p], false);
359 33599894 : 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 1861504 : break;
365 : }
366 :
367 : //------------------------------------------------------------
368 : // 3D
369 23177713 : 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 25288349 : 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 25131435 : if (calculate_dxyz)
417 : {
418 30148509 : std::vector<std::vector<Real>> * comps[3]
419 30148509 : { &this->dphidxi_map, &this->dphideta_map, &this->dphidzeta_map };
420 22599465 : FEInterface::all_shape_derivs(Dim, map_fe_type, elem, qp, comps, false);
421 : }
422 :
423 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
424 25131435 : 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 2110636 : break;
439 : }
440 :
441 : default:
442 : libmesh_error_msg("Invalid Dim = " << Dim);
443 : }
444 47779483 : }
445 :
446 :
447 :
448 165911836 : 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 13954603 : libmesh_assert(elem);
456 13954603 : 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 13954603 : libmesh_ignore(compute_second_derivatives);
463 :
464 165911836 : if (calculate_xyz)
465 3085998 : libmesh_assert_equal_to(phi_map.size(), elem_nodes.size());
466 :
467 165911836 : switch (dim)
468 : {
469 : //--------------------------------------------------------------------
470 : // 0D
471 198780 : case 0:
472 : {
473 19049 : libmesh_assert(elem_nodes[0]);
474 198780 : if (calculate_xyz)
475 0 : xyz[p] = *elem_nodes[0];
476 198780 : if (calculate_dxyz)
477 : {
478 198780 : jac[p] = 1.0;
479 236878 : JxW[p] = qw[p];
480 : }
481 19049 : break;
482 : }
483 :
484 : //--------------------------------------------------------------------
485 : // 1D
486 44981070 : case 1:
487 : {
488 : // Clear the entities that will be summed
489 44981070 : if (calculate_xyz)
490 126767 : xyz[p].zero();
491 44981070 : if (calculate_dxyz)
492 42646910 : dxyzdxi_map[p].zero();
493 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
494 44981070 : if (calculate_d2xyz)
495 : {
496 42506022 : d2xyzdxi2_map[p].zero();
497 : // Inverse map second derivatives
498 45132035 : d2xidxyz2_map[p].assign(6, 0.);
499 : }
500 : #endif
501 :
502 : // compute x, dx, d2x at the quadrature point
503 135176239 : for (auto i : index_range(elem_nodes)) // sum over the nodes
504 : {
505 : // Reference to the point, helps eliminate
506 : // excessive temporaries in the inner loop
507 6135827 : libmesh_assert(elem_nodes[i]);
508 90195169 : const Point & elem_point = *elem_nodes[i];
509 :
510 90195169 : if (calculate_xyz)
511 440855 : xyz[p].add_scaled (elem_point, phi_map[i][p] );
512 90195169 : if (calculate_dxyz)
513 96048227 : dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p]);
514 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
515 90195169 : if (calculate_d2xyz)
516 95558351 : d2xyzdxi2_map[p].add_scaled(elem_point, d2phidxi2_map[i][p]);
517 : #endif
518 : }
519 :
520 : // Compute the jacobian
521 : //
522 : // 1D elements can live in 2D or 3D space.
523 : // The transformation matrix from local->global
524 : // coordinates is
525 : //
526 : // T = | dx/dxi |
527 : // | dy/dxi |
528 : // | dz/dxi |
529 : //
530 : // The generalized determinant of T (from the
531 : // so-called "normal" eqns.) is
532 : // jac = "det(T)" = sqrt(det(T'T))
533 : //
534 : // where T'= transpose of T, so
535 : //
536 : // jac = sqrt( (dx/dxi)^2 + (dy/dxi)^2 + (dz/dxi)^2 )
537 :
538 44981070 : if (calculate_dxyz)
539 : {
540 45284639 : jac[p] = dxyzdxi_map[p].norm();
541 :
542 42646910 : if (jac[p] <= jacobian_tolerance)
543 : {
544 : // Don't call print_info() recursively if we're already
545 : // failing. print_info() calls Elem::volume() which may
546 : // call FE::reinit() and trigger the same failure again.
547 : static bool failing = false;
548 0 : if (!failing)
549 : {
550 0 : failing = true;
551 0 : elem->print_info(libMesh::err);
552 0 : failing = false;
553 0 : if (calculate_xyz)
554 : {
555 0 : libmesh_error_msg("ERROR: negative Jacobian " \
556 : << jac[p] \
557 : << " at point " \
558 : << xyz[p] \
559 : << " in element " \
560 : << elem->id());
561 : }
562 : else
563 : {
564 : // In this case xyz[p] is not defined, so don't
565 : // try to print it out.
566 0 : libmesh_error_msg("ERROR: negative Jacobian " \
567 : << jac[p] \
568 : << " at point index " \
569 : << p \
570 : << " in element " \
571 : << elem->id());
572 : }
573 : }
574 : else
575 : {
576 : // We were already failing when we called this, so just
577 : // stop the current computation and return with
578 : // incomplete results.
579 0 : return;
580 : }
581 : }
582 :
583 : // The inverse Jacobian entries also come from the
584 : // generalized inverse of T (see also the 2D element
585 : // living in 3D code).
586 42646910 : const Real jacm2 = 1./jac[p]/jac[p];
587 45284639 : dxidx_map[p] = jacm2*dxdxi_map(p);
588 : #if LIBMESH_DIM > 1
589 45284639 : dxidy_map[p] = jacm2*dydxi_map(p);
590 : #endif
591 : #if LIBMESH_DIM > 2
592 42646910 : dxidz_map[p] = jacm2*dzdxi_map(p);
593 : #endif
594 :
595 47922368 : JxW[p] = jac[p]*qw[p];
596 : }
597 :
598 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
599 :
600 44981070 : if (calculate_d2xyz)
601 : {
602 : #if LIBMESH_DIM == 1
603 : // Compute inverse map second derivatives for 1D-element-living-in-1D case
604 : this->compute_inverse_map_second_derivs(p);
605 : #elif LIBMESH_DIM == 2
606 : // Compute inverse map second derivatives for 1D-element-living-in-2D case
607 : // See JWP notes for details
608 :
609 : // numer = x_xi*x_{xi xi} + y_xi*y_{xi xi}
610 : Real numer =
611 : dxyzdxi_map[p](0)*d2xyzdxi2_map[p](0) +
612 : dxyzdxi_map[p](1)*d2xyzdxi2_map[p](1);
613 :
614 : // denom = (x_xi)^2 + (y_xi)^2 must be >= 0.0
615 : Real denom =
616 : dxyzdxi_map[p](0)*dxyzdxi_map[p](0) +
617 : dxyzdxi_map[p](1)*dxyzdxi_map[p](1);
618 :
619 : if (denom <= 0.0)
620 : {
621 : // Don't call print_info() recursively if we're already
622 : // failing. print_info() calls Elem::volume() which may
623 : // call FE::reinit() and trigger the same failure again.
624 : static bool failing = false;
625 : if (!failing)
626 : {
627 : failing = true;
628 : elem->print_info(libMesh::err);
629 : failing = false;
630 : libmesh_error_msg("Encountered invalid 1D element!");
631 : }
632 : else
633 : {
634 : // We were already failing when we called this, so just
635 : // stop the current computation and return with
636 : // incomplete results.
637 : return;
638 : }
639 : }
640 :
641 : // xi_{x x}
642 : d2xidxyz2_map[p][0] = -numer * dxidx_map[p]*dxidx_map[p] / denom;
643 :
644 : // xi_{x y}
645 : d2xidxyz2_map[p][1] = -numer * dxidx_map[p]*dxidy_map[p] / denom;
646 :
647 : // xi_{y y}
648 : d2xidxyz2_map[p][3] = -numer * dxidy_map[p]*dxidy_map[p] / denom;
649 :
650 : #elif LIBMESH_DIM == 3
651 : // Compute inverse map second derivatives for 1D-element-living-in-3D case
652 : // See JWP notes for details
653 :
654 : // numer = x_xi*x_{xi xi} + y_xi*y_{xi xi} + z_xi*z_{xi xi}
655 : Real numer =
656 45132035 : dxyzdxi_map[p](0)*d2xyzdxi2_map[p](0) +
657 42506022 : dxyzdxi_map[p](1)*d2xyzdxi2_map[p](1) +
658 42506022 : dxyzdxi_map[p](2)*d2xyzdxi2_map[p](2);
659 :
660 : // denom = (x_xi)^2 + (y_xi)^2 + (z_xi)^2 must be >= 0.0
661 : Real denom =
662 45132035 : dxyzdxi_map[p](0)*dxyzdxi_map[p](0) +
663 42506022 : dxyzdxi_map[p](1)*dxyzdxi_map[p](1) +
664 42506022 : dxyzdxi_map[p](2)*dxyzdxi_map[p](2);
665 :
666 42506022 : if (denom <= 0.0)
667 : {
668 : // Don't call print_info() recursively if we're already
669 : // failing. print_info() calls Elem::volume() which may
670 : // call FE::reinit() and trigger the same failure again.
671 : static bool failing = false;
672 0 : if (!failing)
673 : {
674 0 : failing = true;
675 0 : elem->print_info(libMesh::err);
676 0 : failing = false;
677 0 : libmesh_error_msg("Encountered invalid 1D element!");
678 : }
679 : else
680 : {
681 : // We were already failing when we called this, so just
682 : // stop the current computation and return with
683 : // incomplete results.
684 0 : return;
685 : }
686 : }
687 :
688 : // xi_{x x}
689 45132035 : d2xidxyz2_map[p][0] = -numer * dxidx_map[p]*dxidx_map[p] / denom;
690 :
691 : // xi_{x y}
692 42506022 : d2xidxyz2_map[p][1] = -numer * dxidx_map[p]*dxidy_map[p] / denom;
693 :
694 : // xi_{x z}
695 42506022 : d2xidxyz2_map[p][2] = -numer * dxidx_map[p]*dxidz_map[p] / denom;
696 :
697 : // xi_{y y}
698 42506022 : d2xidxyz2_map[p][3] = -numer * dxidy_map[p]*dxidy_map[p] / denom;
699 :
700 : // xi_{y z}
701 42506022 : d2xidxyz2_map[p][4] = -numer * dxidy_map[p]*dxidz_map[p] / denom;
702 :
703 : // xi_{z z}
704 42506022 : d2xidxyz2_map[p][5] = -numer * dxidz_map[p]*dxidz_map[p] / denom;
705 : #endif //LIBMESH_DIM == 3
706 : } // calculate_d2xyz
707 :
708 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
709 :
710 : // done computing the map
711 3058276 : break;
712 : }
713 :
714 :
715 : //--------------------------------------------------------------------
716 : // 2D
717 87943234 : case 2:
718 : {
719 : //------------------------------------------------------------------
720 : // Compute the (x,y) values at the quadrature points,
721 : // the Jacobian at the quadrature points
722 :
723 87943234 : if (calculate_xyz)
724 17312585 : xyz[p].zero();
725 :
726 87943234 : if (calculate_dxyz)
727 : {
728 80693940 : dxyzdxi_map[p].zero();
729 14994016 : dxyzdeta_map[p].zero();
730 : }
731 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
732 87943234 : if (calculate_d2xyz)
733 : {
734 1357749 : d2xyzdxi2_map[p].zero();
735 230462 : d2xyzdxideta_map[p].zero();
736 230462 : d2xyzdeta2_map[p].zero();
737 : // Inverse map second derivatives
738 1472980 : d2xidxyz2_map[p].assign(6, 0.);
739 1472980 : d2etadxyz2_map[p].assign(6, 0.);
740 : }
741 : #endif
742 :
743 :
744 : // compute (x,y) at the quadrature points, derivatives once
745 521424287 : for (auto i : index_range(elem_nodes)) // sum over the nodes
746 : {
747 : // Reference to the point, helps eliminate
748 : // excessive temporaries in the inner loop
749 39894068 : libmesh_assert(elem_nodes[i]);
750 433481053 : const Point & elem_point = *elem_nodes[i];
751 :
752 433481053 : if (calculate_xyz)
753 131918204 : xyz[p].add_scaled (elem_point, phi_map[i][p] );
754 :
755 433481053 : if (calculate_dxyz)
756 : {
757 444230097 : dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
758 444230097 : dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p]);
759 : }
760 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
761 433481053 : if (calculate_d2xyz)
762 : {
763 10968110 : d2xyzdxi2_map[p].add_scaled (elem_point, d2phidxi2_map[i][p]);
764 10968110 : d2xyzdxideta_map[p].add_scaled (elem_point, d2phidxideta_map[i][p]);
765 10968110 : d2xyzdeta2_map[p].add_scaled (elem_point, d2phideta2_map[i][p]);
766 : }
767 : #endif
768 : }
769 :
770 87943234 : if (calculate_dxyz)
771 : {
772 : // compute the jacobian once
773 7525918 : const Real dx_dxi = dxdxi_map(p),
774 7525918 : dx_deta = dxdeta_map(p),
775 7525918 : dy_dxi = dydxi_map(p),
776 7525918 : dy_deta = dydeta_map(p);
777 :
778 : #if LIBMESH_DIM == 2
779 : // Compute the Jacobian. This assumes the 2D face
780 : // lives in 2D space
781 : //
782 : // Symbolically, the matrix determinant is
783 : //
784 : // | dx/dxi dx/deta |
785 : // jac = | dy/dxi dy/deta |
786 : //
787 : // jac = dx/dxi*dy/deta - dx/deta*dy/dxi
788 : jac[p] = (dx_dxi*dy_deta - dx_deta*dy_dxi);
789 :
790 : if (jac[p] <= jacobian_tolerance)
791 : {
792 : // Don't call print_info() recursively if we're already
793 : // failing. print_info() calls Elem::volume() which may
794 : // call FE::reinit() and trigger the same failure again.
795 : static bool failing = false;
796 : if (!failing)
797 : {
798 : failing = true;
799 : elem->print_info(libMesh::err);
800 : failing = false;
801 : if (calculate_xyz)
802 : {
803 : libmesh_error_msg("ERROR: negative Jacobian " \
804 : << jac[p] \
805 : << " at point " \
806 : << xyz[p] \
807 : << " in element " \
808 : << elem->id());
809 : }
810 : else
811 : {
812 : // In this case xyz[p] is not defined, so don't
813 : // try to print it out.
814 : libmesh_error_msg("ERROR: negative Jacobian " \
815 : << jac[p] \
816 : << " at point index " \
817 : << p \
818 : << " in element " \
819 : << elem->id());
820 : }
821 : }
822 : else
823 : {
824 : // We were already failing when we called this, so just
825 : // stop the current computation and return with
826 : // incomplete results.
827 : return;
828 : }
829 : }
830 :
831 : JxW[p] = jac[p]*qw[p];
832 :
833 : // Compute the shape function derivatives wrt x,y at the
834 : // quadrature points
835 : const Real inv_jac = 1./jac[p];
836 :
837 : dxidx_map[p] = dy_deta*inv_jac; //dxi/dx = (1/J)*dy/deta
838 : dxidy_map[p] = -dx_deta*inv_jac; //dxi/dy = -(1/J)*dx/deta
839 : detadx_map[p] = -dy_dxi* inv_jac; //deta/dx = -(1/J)*dy/dxi
840 : detady_map[p] = dx_dxi* inv_jac; //deta/dy = (1/J)*dx/dxi
841 :
842 : dxidz_map[p] = detadz_map[p] = 0.;
843 :
844 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
845 : if (compute_second_derivatives)
846 : this->compute_inverse_map_second_derivs(p);
847 : #endif
848 : #else // LIBMESH_DIM == 3
849 :
850 7525918 : const Real dz_dxi = dzdxi_map(p),
851 7525918 : dz_deta = dzdeta_map(p);
852 :
853 : // Compute the Jacobian. This assumes a 2D face in
854 : // 3D space.
855 : //
856 : // The transformation matrix T from local to global
857 : // coordinates is
858 : //
859 : // | dx/dxi dx/deta |
860 : // T = | dy/dxi dy/deta |
861 : // | dz/dxi dz/deta |
862 : // note det(T' T) = det(T')det(T) = det(T)det(T)
863 : // so det(T) = std::sqrt(det(T' T))
864 : //
865 : //----------------------------------------------
866 : // Notes:
867 : //
868 : // dX = R dXi -> R'dX = R'R dXi
869 : // (R^-1)dX = dXi [(R'R)^-1 R']dX = dXi
870 : //
871 : // so R^-1 = (R'R)^-1 R'
872 : //
873 : // and R^-1 R = (R'R)^-1 R'R = I.
874 : //
875 110681972 : const Real g11 = (dx_dxi*dx_dxi +
876 80693940 : dy_dxi*dy_dxi +
877 80693940 : dz_dxi*dz_dxi);
878 :
879 110681972 : const Real g12 = (dx_dxi*dx_deta +
880 80693940 : dy_dxi*dy_deta +
881 80693940 : dz_dxi*dz_deta);
882 :
883 7525918 : const Real g21 = g12;
884 :
885 110681972 : const Real g22 = (dx_deta*dx_deta +
886 80693940 : dy_deta*dy_deta +
887 80693940 : dz_deta*dz_deta);
888 :
889 80693940 : const Real det = (g11*g22 - g12*g21);
890 :
891 80693940 : if (det <= jacobian_tolerance)
892 : {
893 : // Don't call print_info() recursively if we're already
894 : // failing. print_info() calls Elem::volume() which may
895 : // call FE::reinit() and trigger the same failure again.
896 : thread_local bool failing = false;
897 0 : if (!failing)
898 : {
899 0 : failing = true;
900 0 : Threads::spin_mutex::scoped_lock lock(_point_inv_err_mutex);
901 0 : elem->print_info(libMesh::err);
902 0 : failing = false;
903 0 : if (calculate_xyz)
904 : {
905 0 : libmesh_error_msg("ERROR: negative Jacobian " \
906 : << det \
907 : << " at point " \
908 : << xyz[p] \
909 : << " in element " \
910 : << elem->id());
911 : }
912 : else
913 : {
914 : // In this case xyz[p] is not defined, so don't
915 : // try to print it out.
916 0 : libmesh_error_msg("ERROR: negative Jacobian " \
917 : << det \
918 : << " at point index " \
919 : << p \
920 : << " in element " \
921 : << elem->id());
922 : }
923 : }
924 : else
925 : {
926 : // We were already failing when we called this, so just
927 : // stop the current computation and return with
928 : // incomplete results.
929 0 : return;
930 : }
931 : }
932 :
933 80693940 : const Real inv_det = 1./det;
934 80693940 : jac[p] = std::sqrt(det);
935 :
936 88162038 : JxW[p] = jac[p]*qw[p];
937 :
938 80693940 : const Real g11inv = g22*inv_det;
939 80693940 : const Real g12inv = -g12*inv_det;
940 7525918 : const Real g21inv = -g21*inv_det;
941 80693940 : const Real g22inv = g11*inv_det;
942 :
943 80693940 : dxidx_map[p] = g11inv*dx_dxi + g12inv*dx_deta;
944 80693940 : dxidy_map[p] = g11inv*dy_dxi + g12inv*dy_deta;
945 80693940 : dxidz_map[p] = g11inv*dz_dxi + g12inv*dz_deta;
946 :
947 80693940 : detadx_map[p] = g21inv*dx_dxi + g22inv*dx_deta;
948 80693940 : detady_map[p] = g21inv*dy_dxi + g22inv*dy_deta;
949 80693940 : detadz_map[p] = g21inv*dz_dxi + g22inv*dz_deta;
950 :
951 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
952 :
953 80693940 : if (calculate_d2xyz)
954 : {
955 : // Compute inverse map second derivative values for
956 : // 2D-element-living-in-3D case. We pursue a least-squares
957 : // solution approach for this "non-square" case, see JWP notes
958 : // for details.
959 :
960 : // A = [ x_{xi xi} x_{eta eta} ]
961 : // [ y_{xi xi} y_{eta eta} ]
962 : // [ z_{xi xi} z_{eta eta} ]
963 1588211 : DenseMatrix<Real> A(3,2);
964 1588211 : A(0,0) = d2xyzdxi2_map[p](0); A(0,1) = d2xyzdeta2_map[p](0);
965 1472980 : A(1,0) = d2xyzdxi2_map[p](1); A(1,1) = d2xyzdeta2_map[p](1);
966 1472980 : A(2,0) = d2xyzdxi2_map[p](2); A(2,1) = d2xyzdeta2_map[p](2);
967 :
968 : // J^T, the transpose of the Jacobian matrix
969 1588211 : DenseMatrix<Real> JT(2,3);
970 1357749 : JT(0,0) = dx_dxi; JT(0,1) = dy_dxi; JT(0,2) = dz_dxi;
971 1472980 : JT(1,0) = dx_deta; JT(1,1) = dy_deta; JT(1,2) = dz_deta;
972 :
973 : // (J^T J)^(-1), this has already been computed for us above...
974 1588211 : DenseMatrix<Real> JTJinv(2,2);
975 1357749 : JTJinv(0,0) = g11inv; JTJinv(0,1) = g12inv;
976 1357749 : JTJinv(1,0) = g21inv; JTJinv(1,1) = g22inv;
977 :
978 : // Some helper variables
979 : RealVectorValue
980 460924 : dxi (dxidx_map[p], dxidy_map[p], dxidz_map[p]),
981 460924 : deta (detadx_map[p], detady_map[p], detadz_map[p]);
982 :
983 : // To be filled in below
984 1357749 : DenseVector<Real> tmp1(2);
985 1357749 : DenseVector<Real> tmp2(3);
986 1357749 : DenseVector<Real> tmp3(2);
987 :
988 : // For (s,t) in {(x,x), (x,y), (x,z), (y,y), (y,z), (z,z)}, compute the
989 : // vector of inverse map second derivatives [xi_{s t}, eta_{s t}]
990 115231 : unsigned ctr=0;
991 5430996 : for (unsigned s=0; s<3; ++s)
992 12219741 : for (unsigned t=s; t<3; ++t)
993 : {
994 : // Construct tmp1 = [xi_s*xi_t, eta_s*eta_t]
995 8146494 : tmp1(0) = dxi(s)*dxi(t);
996 8146494 : tmp1(1) = deta(s)*deta(t);
997 :
998 : // Compute tmp2 = A * tmp1
999 8146494 : A.vector_mult(tmp2, tmp1);
1000 :
1001 : // Compute scalar value "alpha"
1002 8146494 : Real alpha = dxi(s)*deta(t) + deta(s)*dxi(t);
1003 :
1004 : // Compute tmp2 <- tmp2 + alpha * x_{xi eta}
1005 32585976 : for (unsigned i=0; i<3; ++i)
1006 28587798 : tmp2(i) += alpha*d2xyzdxideta_map[p](i);
1007 :
1008 : // Compute tmp3 = J^T * tmp2
1009 8146494 : JT.vector_mult(tmp3, tmp2);
1010 :
1011 : // Compute tmp1 = (J^T J)^(-1) * tmp3. tmp1 is available for us to reuse.
1012 8146494 : JTJinv.vector_mult(tmp1, tmp3);
1013 :
1014 : // Fill in appropriate entries, don't forget to multiply by -1!
1015 8837880 : d2xidxyz2_map[p][ctr] = -tmp1(0);
1016 8837880 : d2etadxyz2_map[p][ctr] = -tmp1(1);
1017 :
1018 : // Increment the counter
1019 8146494 : ctr++;
1020 : }
1021 1127287 : }
1022 :
1023 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1024 :
1025 : #endif // LIBMESH_DIM == 3
1026 : }
1027 : // done computing the map
1028 8143247 : break;
1029 : }
1030 :
1031 :
1032 :
1033 : //--------------------------------------------------------------------
1034 : // 3D
1035 32788752 : case 3:
1036 : {
1037 : //------------------------------------------------------------------
1038 : // Compute the (x,y,z) values at the quadrature points,
1039 : // the Jacobian at the quadrature point
1040 :
1041 : // Clear the entities that will be summed
1042 32788752 : if (calculate_xyz)
1043 16683527 : xyz[p].zero ();
1044 32788752 : if (calculate_dxyz)
1045 : {
1046 28693497 : dxyzdxi_map[p].zero ();
1047 4788902 : dxyzdeta_map[p].zero ();
1048 4788902 : dxyzdzeta_map[p].zero ();
1049 : }
1050 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1051 32788752 : if (calculate_d2xyz)
1052 : {
1053 6135761 : d2xyzdxi2_map[p].zero();
1054 1001996 : d2xyzdxideta_map[p].zero();
1055 1001996 : d2xyzdxidzeta_map[p].zero();
1056 1001996 : d2xyzdeta2_map[p].zero();
1057 1001996 : d2xyzdetadzeta_map[p].zero();
1058 1001996 : d2xyzdzeta2_map[p].zero();
1059 : // Inverse map second derivatives
1060 6636759 : d2xidxyz2_map[p].assign(6, 0.);
1061 6636759 : d2etadxyz2_map[p].assign(6, 0.);
1062 6636759 : d2zetadxyz2_map[p].assign(6, 0.);
1063 : }
1064 : #endif
1065 :
1066 :
1067 : // compute (x,y,z) at the quadrature points,
1068 : // dxdxi, dydxi, dzdxi,
1069 : // dxdeta, dydeta, dzdeta,
1070 : // dxdzeta, dydzeta, dzdzeta all once
1071 503129509 : for (auto i : index_range(elem_nodes)) // sum over the nodes
1072 : {
1073 : // Reference to the point, helps eliminate
1074 : // excessive temporaries in the inner loop
1075 39184488 : libmesh_assert(elem_nodes[i]);
1076 470340757 : const Point & elem_point = *elem_nodes[i];
1077 :
1078 470340757 : if (calculate_xyz)
1079 277453431 : xyz[p].add_scaled (elem_point, phi_map[i][p] );
1080 470340757 : if (calculate_dxyz)
1081 : {
1082 472442733 : dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
1083 472442733 : dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p] );
1084 472442733 : dxyzdzeta_map[p].add_scaled (elem_point, dphidzeta_map[i][p]);
1085 : }
1086 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1087 470340757 : if (calculate_d2xyz)
1088 : {
1089 100131912 : d2xyzdxi2_map[p].add_scaled (elem_point,
1090 23004636 : d2phidxi2_map[i][p]);
1091 100131912 : d2xyzdxideta_map[p].add_scaled (elem_point,
1092 23004636 : d2phidxideta_map[i][p]);
1093 100131912 : d2xyzdxidzeta_map[p].add_scaled (elem_point,
1094 23004636 : d2phidxidzeta_map[i][p]);
1095 100131912 : d2xyzdeta2_map[p].add_scaled (elem_point,
1096 23004636 : d2phideta2_map[i][p]);
1097 100131912 : d2xyzdetadzeta_map[p].add_scaled (elem_point,
1098 23004636 : d2phidetadzeta_map[i][p]);
1099 100131912 : d2xyzdzeta2_map[p].add_scaled (elem_point,
1100 23004636 : d2phidzeta2_map[i][p]);
1101 : }
1102 : #endif
1103 : }
1104 :
1105 32788752 : if (calculate_dxyz)
1106 : {
1107 : // compute the jacobian
1108 : const Real
1109 2394451 : dx_dxi = dxdxi_map(p), dy_dxi = dydxi_map(p), dz_dxi = dzdxi_map(p),
1110 2394451 : dx_deta = dxdeta_map(p), dy_deta = dydeta_map(p), dz_deta = dzdeta_map(p),
1111 2394451 : dx_dzeta = dxdzeta_map(p), dy_dzeta = dydzeta_map(p), dz_dzeta = dzdzeta_map(p);
1112 :
1113 : // Symbolically, the matrix determinant is
1114 : //
1115 : // | dx/dxi dy/dxi dz/dxi |
1116 : // jac = | dx/deta dy/deta dz/deta |
1117 : // | dx/dzeta dy/dzeta dz/dzeta |
1118 : //
1119 : // jac = dx/dxi*(dy/deta*dz/dzeta - dz/deta*dy/dzeta) +
1120 : // dy/dxi*(dz/deta*dx/dzeta - dx/deta*dz/dzeta) +
1121 : // dz/dxi*(dx/deta*dy/dzeta - dy/deta*dx/dzeta)
1122 :
1123 43060203 : jac[p] = (dx_dxi*(dy_deta*dz_dzeta - dz_deta*dy_dzeta) +
1124 33482399 : dy_dxi*(dz_deta*dx_dzeta - dx_deta*dz_dzeta) +
1125 28693497 : dz_dxi*(dx_deta*dy_dzeta - dy_deta*dx_dzeta));
1126 :
1127 28693497 : if (jac[p] <= jacobian_tolerance)
1128 : {
1129 : // Don't call print_info() recursively if we're already
1130 : // failing. print_info() calls Elem::volume() which may
1131 : // call FE::reinit() and trigger the same failure again.
1132 : static bool failing = false;
1133 0 : if (!failing)
1134 : {
1135 0 : failing = true;
1136 0 : elem->print_info(libMesh::err);
1137 0 : failing = false;
1138 0 : if (calculate_xyz)
1139 : {
1140 0 : libmesh_error_msg("ERROR: negative Jacobian " \
1141 : << jac[p] \
1142 : << " at point " \
1143 : << xyz[p] \
1144 : << " in element " \
1145 : << elem->id());
1146 : }
1147 : else
1148 : {
1149 : // In this case xyz[p] is not defined, so don't
1150 : // try to print it out.
1151 0 : libmesh_error_msg("ERROR: negative Jacobian " \
1152 : << jac[p] \
1153 : << " at point index " \
1154 : << p \
1155 : << " in element " \
1156 : << elem->id());
1157 : }
1158 : }
1159 : else
1160 : {
1161 : // We were already failing when we called this, so just
1162 : // stop the current computation and return with
1163 : // incomplete results.
1164 0 : return;
1165 : }
1166 : }
1167 :
1168 31087948 : JxW[p] = jac[p]*qw[p];
1169 :
1170 : // Compute the shape function derivatives wrt x,y at the
1171 : // quadrature points
1172 28693497 : const Real inv_jac = 1./jac[p];
1173 :
1174 28693497 : dxidx_map[p] = (dy_deta*dz_dzeta - dz_deta*dy_dzeta)*inv_jac;
1175 28693497 : dxidy_map[p] = (dz_deta*dx_dzeta - dx_deta*dz_dzeta)*inv_jac;
1176 28693497 : dxidz_map[p] = (dx_deta*dy_dzeta - dy_deta*dx_dzeta)*inv_jac;
1177 :
1178 28693497 : detadx_map[p] = (dz_dxi*dy_dzeta - dy_dxi*dz_dzeta )*inv_jac;
1179 28693497 : detady_map[p] = (dx_dxi*dz_dzeta - dz_dxi*dx_dzeta )*inv_jac;
1180 28693497 : detadz_map[p] = (dy_dxi*dx_dzeta - dx_dxi*dy_dzeta )*inv_jac;
1181 :
1182 28693497 : dzetadx_map[p] = (dy_dxi*dz_deta - dz_dxi*dy_deta )*inv_jac;
1183 28693497 : dzetady_map[p] = (dz_dxi*dx_deta - dx_dxi*dz_deta )*inv_jac;
1184 31087948 : dzetadz_map[p] = (dx_dxi*dy_deta - dy_dxi*dx_deta )*inv_jac;
1185 : }
1186 :
1187 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1188 32788752 : if (compute_second_derivatives)
1189 281008 : this->compute_inverse_map_second_derivs(p);
1190 : #endif
1191 : // done computing the map
1192 2734031 : break;
1193 : }
1194 :
1195 0 : default:
1196 0 : libmesh_error_msg("Invalid dim = " << dim);
1197 : }
1198 : }
1199 :
1200 :
1201 :
1202 161365203 : void FEMap::resize_quadrature_map_vectors(const unsigned int dim, unsigned int n_qp)
1203 : {
1204 : // We're calculating now!
1205 13575032 : this->determine_calculations();
1206 :
1207 : // Resize the vectors to hold data at the quadrature points
1208 161365203 : if (calculate_xyz)
1209 33746293 : xyz.resize(n_qp);
1210 161365203 : if (calculate_dxyz)
1211 : {
1212 147708766 : dxyzdxi_map.resize(n_qp);
1213 147708766 : dxidx_map.resize(n_qp);
1214 147708766 : dxidy_map.resize(n_qp); // 1D element may live in 2D ...
1215 147708766 : dxidz_map.resize(n_qp); // ... or 3D
1216 : }
1217 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1218 161365203 : if (calculate_d2xyz)
1219 : {
1220 49817876 : d2xyzdxi2_map.resize(n_qp);
1221 :
1222 : // Inverse map second derivatives
1223 49817876 : d2xidxyz2_map.resize(n_qp);
1224 245892966 : for (auto i : index_range(d2xidxyz2_map))
1225 208733307 : d2xidxyz2_map[i].assign(6, 0.);
1226 : }
1227 : #endif
1228 161365203 : if (dim > 1)
1229 : {
1230 116187879 : if (calculate_dxyz)
1231 : {
1232 104865602 : dxyzdeta_map.resize(n_qp);
1233 104865602 : detadx_map.resize(n_qp);
1234 104865602 : detady_map.resize(n_qp);
1235 104865602 : detadz_map.resize(n_qp);
1236 : }
1237 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1238 116187879 : if (calculate_d2xyz)
1239 : {
1240 7311854 : d2xyzdxideta_map.resize(n_qp);
1241 7311854 : d2xyzdeta2_map.resize(n_qp);
1242 :
1243 : // Inverse map second derivatives
1244 7311854 : d2etadxyz2_map.resize(n_qp);
1245 33413335 : for (auto i : index_range(d2etadxyz2_map))
1246 28259881 : d2etadxyz2_map[i].assign(6, 0.);
1247 : }
1248 : #endif
1249 116187879 : if (dim > 2)
1250 : {
1251 30051557 : if (calculate_dxyz)
1252 : {
1253 25978574 : dxyzdzeta_map.resize (n_qp);
1254 25978574 : dzetadx_map.resize (n_qp);
1255 25978574 : dzetady_map.resize (n_qp);
1256 25978574 : dzetadz_map.resize (n_qp);
1257 : }
1258 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1259 30051557 : if (calculate_d2xyz)
1260 : {
1261 6113849 : d2xyzdxidzeta_map.resize(n_qp);
1262 6113849 : d2xyzdetadzeta_map.resize(n_qp);
1263 6113849 : d2xyzdzeta2_map.resize(n_qp);
1264 :
1265 : // Inverse map second derivatives
1266 6113849 : d2zetadxyz2_map.resize(n_qp);
1267 23924943 : for (auto i : index_range(d2zetadxyz2_map))
1268 19276489 : d2zetadxyz2_map[i].assign(6, 0.);
1269 : }
1270 : #endif
1271 : }
1272 : }
1273 :
1274 161365203 : if (calculate_dxyz)
1275 : {
1276 147708766 : jac.resize(n_qp);
1277 147708766 : JxW.resize(n_qp);
1278 : }
1279 161365203 : }
1280 :
1281 :
1282 :
1283 160523384 : void FEMap::compute_affine_map(const unsigned int dim,
1284 : const std::vector<Real> & qw,
1285 : const Elem * elem)
1286 : {
1287 : // Start logging the map computation.
1288 27020766 : LOG_SCOPE("compute_affine_map()", "FEMap");
1289 :
1290 13510383 : libmesh_assert(elem);
1291 :
1292 26754086 : const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1293 :
1294 : // Resize the vectors to hold data at the quadrature points
1295 160523384 : this->resize_quadrature_map_vectors(dim, n_qp);
1296 :
1297 : // Determine the nodes contributing to element elem
1298 160523384 : unsigned int n_nodes = elem->n_nodes();
1299 160523384 : _elem_nodes.resize(elem->n_nodes());
1300 1073716886 : for (unsigned int i=0; i<n_nodes; i++)
1301 1068951142 : _elem_nodes[i] = elem->node_ptr(i);
1302 :
1303 : // Compute map at quadrature point 0
1304 160523384 : this->compute_single_point_map(dim, qw, elem, 0, _elem_nodes, /*compute_second_derivatives=*/false);
1305 :
1306 : // Compute xyz at all other quadrature points
1307 160523384 : if (calculate_xyz)
1308 182944057 : for (unsigned int p=1; p<n_qp; p++)
1309 : {
1310 149378338 : xyz[p].zero();
1311 1838201417 : for (auto i : index_range(phi_map)) // sum over the nodes
1312 1986038309 : xyz[p].add_scaled (*_elem_nodes[i], phi_map[i][p]);
1313 : }
1314 :
1315 : // Copy other map data from quadrature point 0
1316 160523384 : if (calculate_dxyz)
1317 628203330 : for (unsigned int p=1; p<n_qp; p++) // for each extra quadrature point
1318 : {
1319 481328211 : dxyzdxi_map[p] = dxyzdxi_map[0];
1320 481328211 : dxidx_map[p] = dxidx_map[0];
1321 481328211 : dxidy_map[p] = dxidy_map[0];
1322 481328211 : dxidz_map[p] = dxidz_map[0];
1323 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1324 : // The map should be affine, so second derivatives are zero
1325 481328211 : if (calculate_d2xyz)
1326 19458636 : d2xyzdxi2_map[p] = 0.;
1327 : #endif
1328 481328211 : if (dim > 1)
1329 : {
1330 353769070 : dxyzdeta_map[p] = dxyzdeta_map[0];
1331 353769070 : detadx_map[p] = detadx_map[0];
1332 353769070 : detady_map[p] = detady_map[0];
1333 353769070 : detadz_map[p] = detadz_map[0];
1334 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1335 353769070 : if (calculate_d2xyz)
1336 : {
1337 3084342 : d2xyzdxideta_map[p] = 0.;
1338 3084342 : d2xyzdeta2_map[p] = 0.;
1339 : }
1340 : #endif
1341 353769070 : if (dim > 2)
1342 : {
1343 96886509 : dxyzdzeta_map[p] = dxyzdzeta_map[0];
1344 96886509 : dzetadx_map[p] = dzetadx_map[0];
1345 96886509 : dzetady_map[p] = dzetady_map[0];
1346 96886509 : dzetadz_map[p] = dzetadz_map[0];
1347 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1348 96886509 : if (calculate_d2xyz)
1349 : {
1350 1928794 : d2xyzdxidzeta_map[p] = 0.;
1351 1928794 : d2xyzdetadzeta_map[p] = 0.;
1352 1928794 : d2xyzdzeta2_map[p] = 0.;
1353 : }
1354 : #endif
1355 : }
1356 : }
1357 481328211 : jac[p] = jac[0];
1358 560753497 : JxW[p] = JxW[0] / qw[0] * qw[p];
1359 : }
1360 160523384 : }
1361 :
1362 :
1363 :
1364 0 : void FEMap::compute_null_map(const unsigned int dim,
1365 : const std::vector<Real> & qw)
1366 : {
1367 : // Start logging the map computation.
1368 0 : LOG_SCOPE("compute_null_map()", "FEMap");
1369 :
1370 0 : const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1371 :
1372 : // Resize the vectors to hold data at the quadrature points
1373 0 : this->resize_quadrature_map_vectors(dim, n_qp);
1374 :
1375 : // Compute "fake" xyz
1376 0 : for (unsigned int p=1; p<n_qp; p++)
1377 : {
1378 0 : if (calculate_xyz)
1379 0 : xyz[p].zero();
1380 :
1381 0 : if (calculate_dxyz)
1382 : {
1383 0 : dxyzdxi_map[p] = 0;
1384 0 : dxidx_map[p] = 0;
1385 0 : dxidy_map[p] = 0;
1386 0 : dxidz_map[p] = 0;
1387 : }
1388 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1389 0 : if (calculate_d2xyz)
1390 : {
1391 0 : d2xyzdxi2_map[p] = 0;
1392 : }
1393 : #endif
1394 0 : if (dim > 1)
1395 : {
1396 0 : if (calculate_dxyz)
1397 : {
1398 0 : dxyzdeta_map[p] = 0;
1399 0 : detadx_map[p] = 0;
1400 0 : detady_map[p] = 0;
1401 0 : detadz_map[p] = 0;
1402 : }
1403 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1404 0 : if (calculate_d2xyz)
1405 : {
1406 0 : d2xyzdxideta_map[p] = 0.;
1407 0 : d2xyzdeta2_map[p] = 0.;
1408 : }
1409 : #endif
1410 0 : if (dim > 2)
1411 : {
1412 0 : if (calculate_dxyz)
1413 : {
1414 0 : dxyzdzeta_map[p] = 0;
1415 0 : dzetadx_map[p] = 0;
1416 0 : dzetady_map[p] = 0;
1417 0 : dzetadz_map[p] = 0;
1418 : }
1419 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1420 0 : if (calculate_d2xyz)
1421 : {
1422 0 : d2xyzdxidzeta_map[p] = 0;
1423 0 : d2xyzdetadzeta_map[p] = 0;
1424 0 : d2xyzdzeta2_map[p] = 0;
1425 : }
1426 : #endif
1427 : }
1428 : }
1429 0 : if (calculate_dxyz)
1430 : {
1431 0 : jac[p] = 1;
1432 0 : JxW[p] = qw[p];
1433 : }
1434 : }
1435 0 : }
1436 :
1437 :
1438 :
1439 161365203 : void FEMap::compute_map(const unsigned int dim,
1440 : const std::vector<Real> & qw,
1441 : const Elem * elem,
1442 : bool calculate_d2phi)
1443 : {
1444 161365203 : if (!elem)
1445 : {
1446 0 : compute_null_map(dim, qw);
1447 13510383 : return;
1448 : }
1449 :
1450 161365203 : if (elem->has_affine_map())
1451 : {
1452 160523384 : compute_affine_map(dim, qw, elem);
1453 160523384 : return;
1454 : }
1455 : #ifndef LIBMESH_ENABLE_SECOND_DERIVATIVES
1456 : libmesh_assert(!calculate_d2phi);
1457 : #endif
1458 :
1459 : // Start logging the map computation.
1460 129298 : LOG_SCOPE("compute_map()", "FEMap");
1461 :
1462 64649 : libmesh_assert(elem);
1463 :
1464 129298 : const unsigned int n_qp = cast_int<unsigned int>(qw.size());
1465 :
1466 : // Resize the vectors to hold data at the quadrature points
1467 841819 : this->resize_quadrature_map_vectors(dim, n_qp);
1468 :
1469 : // Determine the nodes contributing to element elem
1470 841819 : if (elem->type() == TRI3SUBDIVISION)
1471 : {
1472 : // Subdivision surface FE require the 1-ring around elem
1473 256 : libmesh_assert_equal_to (dim, 2);
1474 256 : const Tri3Subdivision * sd_elem = static_cast<const Tri3Subdivision *>(elem);
1475 3072 : MeshTools::Subdivision::find_one_ring(sd_elem, _elem_nodes);
1476 : }
1477 : else
1478 : {
1479 : // All other FE use only the nodes of elem itself
1480 838747 : _elem_nodes.resize(elem->n_nodes(), nullptr);
1481 10172852 : for (auto i : elem->node_index_range())
1482 10762049 : _elem_nodes[i] = elem->node_ptr(i);
1483 : }
1484 :
1485 : // Compute map at all quadrature points
1486 6230271 : for (unsigned int p=0; p!=n_qp; p++)
1487 5388452 : this->compute_single_point_map(dim, qw, elem, p, _elem_nodes, calculate_d2phi);
1488 : }
1489 :
1490 :
1491 :
1492 0 : void FEMap::print_JxW(std::ostream & os) const
1493 : {
1494 0 : for (auto i : index_range(JxW))
1495 0 : os << " [" << i << "]: " << JxW[i] << std::endl;
1496 0 : }
1497 :
1498 :
1499 :
1500 0 : void FEMap::print_xyz(std::ostream & os) const
1501 : {
1502 0 : for (auto i : index_range(xyz))
1503 0 : os << " [" << i << "]: " << xyz[i];
1504 0 : }
1505 :
1506 :
1507 :
1508 281008 : void FEMap::compute_inverse_map_second_derivs(unsigned p)
1509 : {
1510 : #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
1511 : // Only certain second derivatives are valid depending on the
1512 : // dimension...
1513 35742 : std::set<unsigned> valid_indices;
1514 :
1515 : // Construct J^{-1}, A, and B matrices (see JWP's notes for details)
1516 : // for cases in which the element dimension matches LIBMESH_DIM.
1517 : #if LIBMESH_DIM==1
1518 : RealTensor
1519 : Jinv(dxidx_map[p], 0., 0.,
1520 : 0., 0., 0.,
1521 : 0., 0., 0.),
1522 :
1523 : A(d2xyzdxi2_map[p](0), 0., 0.,
1524 : 0., 0., 0.,
1525 : 0., 0., 0.),
1526 :
1527 : B(0., 0., 0.,
1528 : 0., 0., 0.,
1529 : 0., 0., 0.);
1530 :
1531 : RealVectorValue
1532 : dxi (dxidx_map[p], 0., 0.),
1533 : deta (0., 0., 0.),
1534 : dzeta(0., 0., 0.);
1535 :
1536 : // In 1D, we have only the xx second derivative
1537 : valid_indices.insert(0);
1538 :
1539 : #elif LIBMESH_DIM==2
1540 : RealTensor
1541 : Jinv(dxidx_map[p], dxidy_map[p], 0.,
1542 : detadx_map[p], detady_map[p], 0.,
1543 : 0., 0., 0.),
1544 :
1545 : A(d2xyzdxi2_map[p](0), d2xyzdeta2_map[p](0), 0.,
1546 : d2xyzdxi2_map[p](1), d2xyzdeta2_map[p](1), 0.,
1547 : 0., 0., 0.),
1548 :
1549 : B(d2xyzdxideta_map[p](0), 0., 0.,
1550 : d2xyzdxideta_map[p](1), 0., 0.,
1551 : 0., 0., 0.);
1552 :
1553 : RealVectorValue
1554 : dxi (dxidx_map[p], dxidy_map[p], 0.),
1555 : deta (detadx_map[p], detady_map[p], 0.),
1556 : dzeta(0., 0., 0.);
1557 :
1558 : // In 2D, we have xx, xy, and yy second derivatives
1559 : const unsigned tmp[3] = {0,1,3};
1560 : valid_indices.insert(tmp, tmp+3);
1561 :
1562 : #elif LIBMESH_DIM==3
1563 : RealTensor
1564 71484 : Jinv(dxidx_map[p], dxidy_map[p], dxidz_map[p],
1565 71484 : detadx_map[p], detady_map[p], detadz_map[p],
1566 352492 : dzetadx_map[p], dzetady_map[p], dzetadz_map[p]),
1567 :
1568 17871 : A(d2xyzdxi2_map[p](0), d2xyzdeta2_map[p](0), d2xyzdzeta2_map[p](0),
1569 17871 : d2xyzdxi2_map[p](1), d2xyzdeta2_map[p](1), d2xyzdzeta2_map[p](1),
1570 107226 : d2xyzdxi2_map[p](2), d2xyzdeta2_map[p](2), d2xyzdzeta2_map[p](2)),
1571 :
1572 17871 : B(d2xyzdxideta_map[p](0), d2xyzdxidzeta_map[p](0), d2xyzdetadzeta_map[p](0),
1573 17871 : d2xyzdxideta_map[p](1), d2xyzdxidzeta_map[p](1), d2xyzdetadzeta_map[p](1),
1574 107226 : d2xyzdxideta_map[p](2), d2xyzdxidzeta_map[p](2), d2xyzdetadzeta_map[p](2));
1575 :
1576 : RealVectorValue
1577 17871 : dxi (dxidx_map[p], dxidy_map[p], dxidz_map[p]),
1578 17871 : deta (detadx_map[p], detady_map[p], detadz_map[p]),
1579 17871 : dzeta(dzetadx_map[p], dzetady_map[p], dzetadz_map[p]);
1580 :
1581 : // In 3D, we have xx, xy, xz, yy, yz, and zz second derivatives
1582 281008 : const unsigned tmp[6] = {0,1,2,3,4,5};
1583 17871 : valid_indices.insert(tmp, tmp+6);
1584 :
1585 : #endif
1586 :
1587 : // For (s,t) in {(x,x), (x,y), (x,z), (y,y), (y,z), (z,z)}, compute the
1588 : // vector of inverse map second derivatives [xi_{s t}, eta_{s t}, zeta_{s t}]
1589 263137 : unsigned ctr=0;
1590 1124032 : for (unsigned s=0; s<3; ++s)
1591 2529072 : for (unsigned t=s; t<3; ++t)
1592 : {
1593 107226 : if (valid_indices.count(ctr))
1594 : {
1595 : RealVectorValue
1596 1686048 : v1(dxi(s)*dxi(t),
1597 1686048 : deta(s)*deta(t),
1598 1793274 : dzeta(s)*dzeta(t)),
1599 :
1600 1686048 : v2(dxi(s)*deta(t) + deta(s)*dxi(t),
1601 1686048 : dxi(s)*dzeta(t) + dzeta(s)*dxi(t),
1602 1793274 : deta(s)*dzeta(t) + dzeta(s)*deta(t));
1603 :
1604 : // Compute the inverse map second derivatives
1605 1686048 : RealVectorValue v3 = -Jinv*(A*v1 + B*v2);
1606 :
1607 : // Store them in the appropriate locations in the class data structures
1608 1793274 : d2xidxyz2_map[p][ctr] = v3(0);
1609 :
1610 : if (LIBMESH_DIM > 1)
1611 1793274 : d2etadxyz2_map[p][ctr] = v3(1);
1612 :
1613 : if (LIBMESH_DIM > 2)
1614 1900500 : d2zetadxyz2_map[p][ctr] = v3(2);
1615 : }
1616 :
1617 : // Increment the counter
1618 1686048 : ctr++;
1619 : }
1620 : #else
1621 : // to avoid compiler warnings:
1622 : libmesh_ignore(p);
1623 : #endif // LIBMESH_ENABLE_SECOND_DERIVATIVES
1624 281008 : }
1625 :
1626 :
1627 :
1628 106911640 : Point FEMap::inverse_map (const unsigned int dim,
1629 : const Elem * elem,
1630 : const Point & physical_point,
1631 : const Real tolerance,
1632 : const bool secure,
1633 : const bool
1634 : extra_checks
1635 : )
1636 : {
1637 7606423 : libmesh_assert(elem);
1638 7606423 : libmesh_assert_greater_equal (tolerance, 0.);
1639 :
1640 7606423 : libmesh_ignore(extra_checks);
1641 :
1642 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1643 :
1644 : // TODO: possibly use the extra_checks parameter in InfFEMap::inverse_map() as well.
1645 :
1646 20908861 : if (elem->infinite())
1647 714 : return InfFEMap::inverse_map(dim, elem, physical_point, tolerance,
1648 711 : secure);
1649 : #endif
1650 :
1651 : // Start logging the map inversion.
1652 15212132 : LOG_SCOPE("inverse_map()", "FEMap");
1653 :
1654 : // How much did the point on the reference
1655 : // element change by in this Newton step?
1656 7606066 : Real inverse_map_error = 0.;
1657 :
1658 : // The point on the reference element. This is
1659 : // the "initial guess" for Newton's method. The
1660 : // centroid seems like a good idea, but computing
1661 : // it is a little more intensive than, say taking
1662 : // the zero point.
1663 : //
1664 : // Convergence should be insensitive of this choice
1665 : // for "good" elements.
1666 7606066 : Point p; // the zero point. No computation required
1667 :
1668 : // The number of iterations in the map inversion process.
1669 7606066 : unsigned int cnt = 0;
1670 :
1671 : // The number of iterations after which we give up and declare
1672 : // divergence
1673 7606066 : const unsigned int max_cnt = 10;
1674 :
1675 : // The distance (in master element space) beyond which we give up
1676 : // and declare divergence. This is no longer used...
1677 : // Real max_step_length = 4.;
1678 :
1679 :
1680 :
1681 : // Newton iteration loop.
1682 7395909 : do
1683 : {
1684 : // Where our current iterate \p p maps to.
1685 212472739 : const Point physical_guess = map(dim, elem, p);
1686 :
1687 : // How far our current iterate is from the actual point.
1688 15001975 : const Point delta = physical_point - physical_guess;
1689 :
1690 : // Increment in current iterate \p p, will be computed.
1691 15001975 : Point dp;
1692 :
1693 :
1694 : // The form of the map and how we invert it depends
1695 : // on the dimension that we are in.
1696 212472739 : switch (dim)
1697 : {
1698 : // ------------------------------------------------------------------
1699 : // 0D map inversion is trivial
1700 1070 : case 0:
1701 : {
1702 1070 : break;
1703 : }
1704 :
1705 : // ------------------------------------------------------------------
1706 : // 1D map inversion
1707 : //
1708 : // Here we find the point on a 1D reference element that maps to
1709 : // the point \p physical_point in the domain. This is a bit tricky
1710 : // since we do not want to assume that the point \p physical_point
1711 : // is also in a 1D domain. In particular, this method might get
1712 : // called on the edge of a 3D element, in which case
1713 : // \p physical_point actually lives in 3D.
1714 8694626 : case 1:
1715 : {
1716 8694626 : const Point dxi = map_deriv (dim, elem, 0, p);
1717 :
1718 : // Newton's method in this case looks like
1719 : //
1720 : // {X} - {X_n} = [J]*dp
1721 : //
1722 : // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x1 matrix
1723 : // d(x,y,z)/dxi, and we seek dp, a scalar. Since the above
1724 : // system is either overdetermined or rank-deficient, we will
1725 : // solve the normal equations for this system
1726 : //
1727 : // [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
1728 : //
1729 : // which involves the trivial inversion of the scalar
1730 : // G = [J]^T [J]
1731 889150 : const Real G = dxi*dxi;
1732 :
1733 889150 : if (secure)
1734 146265 : libmesh_assert_greater (G, 0.);
1735 :
1736 8694626 : const Real Ginv = 1./G;
1737 :
1738 889150 : const Real dxidelta = dxi*delta;
1739 :
1740 8694626 : dp(0) = Ginv*dxidelta;
1741 :
1742 : // No master elements have radius > 4, but sometimes we
1743 : // can take a step that big while still converging
1744 : // if (secure)
1745 : // libmesh_assert_less (dp.size(), max_step_length);
1746 :
1747 889150 : break;
1748 : }
1749 :
1750 :
1751 :
1752 : // ------------------------------------------------------------------
1753 : // 2D map inversion
1754 : //
1755 : // Here we find the point on a 2D reference element that maps to
1756 : // the point \p physical_point in the domain. This is a bit tricky
1757 : // since we do not want to assume that the point \p physical_point
1758 : // is also in a 2D domain. In particular, this method might get
1759 : // called on the face of a 3D element, in which case
1760 : // \p physical_point actually lives in 3D.
1761 99921930 : case 2:
1762 : {
1763 99921930 : const Point dxi = map_deriv (dim, elem, 0, p);
1764 99921930 : const Point deta = map_deriv (dim, elem, 1, p);
1765 :
1766 : // Newton's method in this case looks like
1767 : //
1768 : // {X} - {X_n} = [J]*{dp}
1769 : //
1770 : // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x2 matrix
1771 : // d(x,y,z)/d(xi,eta), and we seek {dp}, a 2x1 vector. Since
1772 : // the above system is either over-determined or rank-deficient,
1773 : // we will solve the normal equations for this system
1774 : //
1775 : // [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
1776 : //
1777 : // which involves the inversion of the 2x2 matrix
1778 : // [G] = [J]^T [J]
1779 : const Real
1780 8449359 : G11 = dxi*dxi, G12 = dxi*deta,
1781 8449359 : G21 = dxi*deta, G22 = deta*deta;
1782 :
1783 :
1784 99921930 : const Real det = (G11*G22 - G12*G21);
1785 :
1786 8449359 : if (secure)
1787 5303021 : libmesh_assert_not_equal_to (det, 0.);
1788 :
1789 99921930 : const Real inv_det = 1./det;
1790 :
1791 : const Real
1792 99921930 : Ginv11 = G22*inv_det,
1793 99921930 : Ginv12 = -G12*inv_det,
1794 :
1795 8449359 : Ginv21 = -G21*inv_det,
1796 99921930 : Ginv22 = G11*inv_det;
1797 :
1798 :
1799 8449359 : const Real dxidelta = dxi*delta;
1800 8449359 : const Real detadelta = deta*delta;
1801 :
1802 99921930 : dp(0) = (Ginv11*dxidelta + Ginv12*detadelta);
1803 99921930 : dp(1) = (Ginv21*dxidelta + Ginv22*detadelta);
1804 :
1805 : // No master elements have radius > 4, but sometimes we
1806 : // can take a step that big while still converging
1807 : // if (secure)
1808 : // libmesh_assert_less (dp.size(), max_step_length);
1809 :
1810 8449359 : break;
1811 : }
1812 :
1813 :
1814 :
1815 : // ------------------------------------------------------------------
1816 : // 3D map inversion
1817 : //
1818 : // Here we find the point in a 3D reference element that maps to
1819 : // the point \p physical_point in a 3D domain. Nothing special
1820 : // has to happen here, since (unless the map is singular because
1821 : // you have a BAD element) the map will be invertible and we can
1822 : // apply Newton's method directly.
1823 103838701 : case 3:
1824 : {
1825 103838701 : const Point dxi = map_deriv (dim, elem, 0, p);
1826 103838701 : const Point deta = map_deriv (dim, elem, 1, p);
1827 103838701 : const Point dzeta = map_deriv (dim, elem, 2, p);
1828 :
1829 : // Newton's method in this case looks like
1830 : //
1831 : // {X} = {X_n} + [J]*{dp}
1832 : //
1833 : // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x3 matrix
1834 : // d(x,y,z)/d(xi,eta,zeta), and we seek {dp}, a 3x1 vector.
1835 : // Since the above system is nonsingular for invertible maps
1836 : // we will solve
1837 : //
1838 : // {dp} = [J]^-1 ({X} - {X_n})
1839 : //
1840 : // which involves the inversion of the 3x3 matrix [J]
1841 : libmesh_try
1842 : {
1843 103868307 : RealTensorValue(dxi(0), deta(0), dzeta(0),
1844 : dxi(1), deta(1), dzeta(1),
1845 103838701 : dxi(2), deta(2), dzeta(2)).solve(delta, dp);
1846 : }
1847 198138 : libmesh_catch (ConvergenceFailure &)
1848 : {
1849 : // If we encountered a singular Jacobian, but are at
1850 : // a singular node, return said singular node.
1851 : const unsigned int local_singular_node =
1852 187738 : elem->local_singular_node(physical_point, tolerance);
1853 187738 : if (local_singular_node != invalid_uint)
1854 : {
1855 240 : libmesh_assert_less(local_singular_node, elem->n_nodes());
1856 5136 : return elem->master_point(local_singular_node);
1857 : }
1858 :
1859 : // We encountered a singular Jacobian. The value of
1860 : // dp is zero, since it was never changed during the
1861 : // call to RealTensorValue::solve(). We don't want to
1862 : // continue iterating until max_cnt since there is no
1863 : // update to the Newton iterate, and we don't want to
1864 : // print the inverse_map_error value since it will
1865 : // confusingly be 0. Therefore, in the secure case we
1866 : // need to throw an error message while in the !secure
1867 : // case we can just return a far away point.
1868 182602 : if (secure)
1869 : {
1870 0 : libMesh::err << "ERROR: Newton scheme encountered a singular Jacobian in element: "
1871 0 : << elem->id()
1872 0 : << std::endl;
1873 :
1874 0 : elem->print_info(libMesh::err);
1875 :
1876 0 : libmesh_error_msg("Exiting...");
1877 : }
1878 : else
1879 : {
1880 730408 : for (unsigned int i=0; i != dim; ++i)
1881 547806 : p(i) = 1e6;
1882 182602 : return p;
1883 : }
1884 177338 : }
1885 :
1886 : // No master elements have radius > 4, but sometimes we
1887 : // can take a step that big while still converging
1888 : // if (secure)
1889 : // libmesh_assert_less (dp.size(), max_step_length);
1890 :
1891 103650963 : break;
1892 : }
1893 :
1894 :
1895 : // Some other dimension?
1896 0 : default:
1897 0 : libmesh_error_msg("Invalid dim = " << dim);
1898 : } // end switch(Dim), dp now computed
1899 :
1900 :
1901 :
1902 : // ||P_n+1 - P_n||
1903 197441129 : inverse_map_error = dp.norm();
1904 :
1905 : // P_n+1 = P_n + dp
1906 14996775 : p.add (dp);
1907 :
1908 : // Increment the iteration count.
1909 212285001 : cnt++;
1910 :
1911 : // Watch for divergence of Newton's
1912 : // method. Here's how it goes:
1913 : // (1) For good elements, we expect convergence in 10
1914 : // iterations, with no too-large steps.
1915 : // - If called with (secure == true) and we have not yet converged
1916 : // print out a warning message.
1917 : // - If called with (secure == true) and we have not converged in
1918 : // 20 iterations abort
1919 : // (2) This method may be called in cases when the target point is not
1920 : // inside the element and we have no business expecting convergence.
1921 : // For these cases if we have not converged in 10 iterations forget
1922 : // about it.
1923 212285001 : if (cnt > max_cnt)
1924 : {
1925 : // Warn about divergence when secure is true - this
1926 : // shouldn't happen
1927 77576 : if (secure)
1928 : {
1929 : // Print every time in devel/dbg modes
1930 : #ifndef NDEBUG
1931 0 : libmesh_here();
1932 0 : libMesh::err << "WARNING: Newton scheme has not converged in "
1933 0 : << cnt << " iterations:" << std::endl
1934 0 : << " physical_point="
1935 0 : << physical_point
1936 0 : << " physical_guess="
1937 0 : << physical_guess
1938 0 : << " dp="
1939 0 : << dp
1940 0 : << " p="
1941 0 : << p
1942 0 : << " error=" << inverse_map_error
1943 0 : << " in element " << elem->id()
1944 0 : << std::endl;
1945 :
1946 0 : elem->print_info(libMesh::err);
1947 : #else
1948 : // In optimized mode, just print once that an inverse_map() call
1949 : // had trouble converging its Newton iteration.
1950 0 : libmesh_do_once(libMesh::err << "WARNING: At least one element took more than "
1951 : << max_cnt
1952 : << " iterations to converge in inverse_map()...\n"
1953 : << "Rerun in devel/dbg mode for more details."
1954 : << std::endl;);
1955 :
1956 : #endif // NDEBUG
1957 :
1958 0 : if (cnt > 2*max_cnt)
1959 : {
1960 0 : libMesh::err << "ERROR: Newton scheme FAILED to converge in "
1961 0 : << cnt
1962 0 : << " iterations in element "
1963 0 : << elem->id()
1964 0 : << " for physical point = "
1965 0 : << physical_point
1966 0 : << std::endl;
1967 :
1968 0 : elem->print_info(libMesh::err);
1969 :
1970 0 : libmesh_error_msg("Exiting...");
1971 : }
1972 : }
1973 : // Return a far off point when secure is false - this
1974 : // should only happen when we're trying to map a point
1975 : // that's outside the element
1976 : else
1977 : {
1978 308737 : for (unsigned int i=0; i != dim; ++i)
1979 231161 : p(i) = 1e6;
1980 :
1981 77576 : return p;
1982 : }
1983 : }
1984 : }
1985 212207425 : while (inverse_map_error > tolerance);
1986 :
1987 :
1988 :
1989 : // If we are in debug mode and the user requested it, do two extra sanity checks.
1990 : #ifdef DEBUG
1991 :
1992 7599161 : if (extra_checks)
1993 : {
1994 : // Make sure the point \p p on the reference element actually
1995 : // does map to the point \p physical_point within a tolerance.
1996 :
1997 3598908 : const Point check = map (dim, elem, p);
1998 3598908 : const Point diff = physical_point - check;
1999 :
2000 3598908 : if (diff.norm() > tolerance)
2001 : {
2002 0 : libmesh_here();
2003 0 : libMesh::err << "WARNING: diff is "
2004 0 : << diff.norm()
2005 0 : << std::endl
2006 0 : << " point="
2007 0 : << physical_point;
2008 0 : libMesh::err << " local=" << check;
2009 0 : libMesh::err << " lref= " << p;
2010 :
2011 0 : elem->print_info(libMesh::err);
2012 : }
2013 :
2014 : // Make sure the point \p p on the reference element actually
2015 : // is
2016 :
2017 3598908 : if (!elem->on_reference_element(p, 2*tolerance))
2018 : {
2019 0 : libmesh_here();
2020 0 : libMesh::err << "WARNING: inverse_map of physical point "
2021 0 : << physical_point
2022 0 : << " is not on element." << '\n';
2023 0 : elem->print_info(libMesh::err);
2024 : }
2025 : }
2026 :
2027 : #endif
2028 :
2029 106644901 : return p;
2030 : }
2031 :
2032 :
2033 :
2034 2081755 : void FEMap::inverse_map (const unsigned int dim,
2035 : const Elem * elem,
2036 : const std::vector<Point> & physical_points,
2037 : std::vector<Point> & reference_points,
2038 : const Real tolerance,
2039 : const bool secure,
2040 : const bool extra_checks)
2041 : {
2042 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2043 699771 : if (elem->infinite())
2044 : {
2045 : // TODO: possibly use the extra_checks parameter in InfFEMap::inverse_map() as well.
2046 0 : libmesh_ignore(extra_checks);
2047 :
2048 0 : return InfFEMap::inverse_map(dim, elem, physical_points, reference_points, tolerance, secure);
2049 : }
2050 : #endif
2051 :
2052 : // The number of points to find the
2053 : // inverse map of
2054 468547 : const std::size_t n_points = physical_points.size();
2055 :
2056 : // Resize the vector to hold the points
2057 : // on the reference element
2058 2081755 : reference_points.resize(n_points);
2059 :
2060 : // Find the coordinates on the reference
2061 : // element of each point in physical space
2062 5213646 : for (std::size_t p=0; p<n_points; p++)
2063 3465432 : reference_points[p] =
2064 3465432 : inverse_map (dim, elem, physical_points[p], tolerance, secure, extra_checks);
2065 1381984 : }
2066 :
2067 :
2068 :
2069 242822912 : Point FEMap::map (const unsigned int dim,
2070 : const Elem * elem,
2071 : const Point & reference_point)
2072 : {
2073 20712656 : libmesh_assert(elem);
2074 20712656 : libmesh_ignore(dim);
2075 :
2076 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
2077 50322524 : if (elem->infinite())
2078 80 : return InfFEMap::map(dim, elem, reference_point);
2079 : #endif
2080 :
2081 20712624 : Point p;
2082 :
2083 242822832 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
2084 242822832 : const FEType fe_type (elem->default_order(), mapping_family);
2085 :
2086 : // Do not consider the Elem::p_level(), if any, when computing the
2087 : // number of shape functions.
2088 242822832 : const unsigned int n_sf = FEInterface::n_shape_functions(fe_type, /*extra_order=*/0, elem);
2089 :
2090 : FEInterface::shape_ptr shape_ptr =
2091 242822832 : FEInterface::shape_function(fe_type, elem);
2092 :
2093 : // Lagrange basis functions are used for mapping
2094 2718497885 : for (unsigned int i=0; i<n_sf; i++)
2095 2475675053 : p.add_scaled (elem->point(i),
2096 2659346590 : shape_ptr(fe_type, elem, i, reference_point, false));
2097 :
2098 242822832 : return p;
2099 : }
2100 :
2101 :
2102 :
2103 565121476 : Point FEMap::map_deriv (const unsigned int dim,
2104 : const Elem * elem,
2105 : const unsigned int j,
2106 : const Point & reference_point)
2107 : {
2108 38894836 : libmesh_assert(elem);
2109 38894836 : libmesh_ignore(dim);
2110 :
2111 108241134 : if (elem->infinite())
2112 0 : libmesh_not_implemented();
2113 :
2114 38894836 : Point p;
2115 :
2116 565121476 : const FEFamily mapping_family = FEMap::map_fe_type(*elem);
2117 565121476 : const FEType fe_type (elem->default_order(), mapping_family);
2118 :
2119 : // Do not consider the Elem::p_level(), if any, when computing the
2120 : // number of shape functions.
2121 : const unsigned int n_sf =
2122 565121476 : FEInterface::n_shape_functions(fe_type, /*total_order=*/0, elem);
2123 :
2124 : FEInterface::shape_deriv_ptr shape_deriv_ptr =
2125 565121476 : FEInterface::shape_deriv_function(fe_type, elem);
2126 :
2127 : // Lagrange basis functions are used for mapping
2128 6928475125 : for (unsigned int i=0; i<n_sf; i++)
2129 6363353649 : p.add_scaled (elem->point(i),
2130 6757685217 : shape_deriv_ptr(fe_type, elem, i, j, reference_point,
2131 : /*add_p_level=*/false));
2132 :
2133 604016312 : return p;
2134 : }
2135 :
2136 :
2137 :
2138 : // Explicit instantiation of FEMap member functions
2139 : template LIBMESH_EXPORT void FEMap::init_reference_to_physical_map<0>( const std::vector<Point> &, const Elem *);
2140 : template LIBMESH_EXPORT void FEMap::init_reference_to_physical_map<1>( const std::vector<Point> &, const Elem *);
2141 : template LIBMESH_EXPORT void FEMap::init_reference_to_physical_map<2>( const std::vector<Point> &, const Elem *);
2142 : template LIBMESH_EXPORT void FEMap::init_reference_to_physical_map<3>( const std::vector<Point> &, const Elem *);
2143 :
2144 : // subdivision elements are implemented only for 2D meshes & reimplement
2145 : // the inverse_maps method separately
2146 : INSTANTIATE_SUBDIVISION_MAPS;
2147 :
2148 : } // namespace libMesh
|
