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