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