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