fe_map.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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"
29 #include "libmesh/fe_macro.h"
30 #include "libmesh/fe_map.h"
31 #include "libmesh/fe_xyz_map.h"
32 
33 namespace libMesh
34 {
35 
36  // Constructor (empty)
38 
39 
40 
42 {
43  switch( fe_type.family )
44  {
45  case XYZ:
46  {
47  AutoPtr<FEMap> ap( new FEXYZMap );
48  return ap;
49  }
50  default:
51  {
52  AutoPtr<FEMap> ap( new FEMap );
53  return ap;
54  }
55  }
56 
57  //Shouldn't ever get here
58  libmesh_error();
59  return AutoPtr<FEMap>();
60 }
61 
62 
63 
64 template<unsigned int Dim>
65 void FEMap::init_reference_to_physical_map( const std::vector<Point>& qp,
66  const Elem* elem)
67 {
68  // Start logging the reference->physical map initialization
69  START_LOG("init_reference_to_physical_map()", "FEMap");
70 
71  // The number of quadrature points.
72  const std::size_t n_qp = qp.size();
73 
74  // The element type and order to use in
75  // the map
76  const Order mapping_order (elem->default_order());
77  const ElemType mapping_elem_type (elem->type());
78 
79  // Number of shape functions used to construt the map
80  // (Lagrange shape functions are used for mapping)
81  const unsigned int n_mapping_shape_functions =
82  FE<Dim,LAGRANGE>::n_shape_functions (mapping_elem_type,
83  mapping_order);
84 
85  this->phi_map.resize (n_mapping_shape_functions);
86  if (Dim > 0)
87  {
88  this->dphidxi_map.resize (n_mapping_shape_functions);
89 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
90  this->d2phidxi2_map.resize (n_mapping_shape_functions);
91 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
92  }
93 
94  if (Dim > 1)
95  {
96  this->dphideta_map.resize (n_mapping_shape_functions);
97 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
98  this->d2phidxideta_map.resize (n_mapping_shape_functions);
99  this->d2phideta2_map.resize (n_mapping_shape_functions);
100 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
101  }
102 
103  if (Dim > 2)
104  {
105  this->dphidzeta_map.resize (n_mapping_shape_functions);
106 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
107  this->d2phidxidzeta_map.resize (n_mapping_shape_functions);
108  this->d2phidetadzeta_map.resize (n_mapping_shape_functions);
109  this->d2phidzeta2_map.resize (n_mapping_shape_functions);
110 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
111  }
112 
113 
114  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
115  {
116  this->phi_map[i].resize (n_qp);
117  if (Dim > 0)
118  {
119  this->dphidxi_map[i].resize (n_qp);
120 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
121  this->d2phidxi2_map[i].resize (n_qp);
122  if (Dim > 1)
123  {
124  this->d2phidxideta_map[i].resize (n_qp);
125  this->d2phideta2_map[i].resize (n_qp);
126  }
127  if (Dim > 2)
128  {
129  this->d2phidxidzeta_map[i].resize (n_qp);
130  this->d2phidetadzeta_map[i].resize (n_qp);
131  this->d2phidzeta2_map[i].resize (n_qp);
132  }
133 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
134 
135  if (Dim > 1)
136  this->dphideta_map[i].resize (n_qp);
137 
138  if (Dim > 2)
139  this->dphidzeta_map[i].resize (n_qp);
140  }
141  }
142 
143  // Optimize for the *linear* geometric elements case:
144  bool is_linear = elem->is_linear();
145 
146  switch (Dim)
147  {
148 
149  //------------------------------------------------------------
150  // 0D
151  case 0:
152  {
153  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
154  for (std::size_t p=0; p<n_qp; p++)
155  this->phi_map[i][p] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
156 
157  break;
158  }
159 
160  //------------------------------------------------------------
161  // 1D
162  case 1:
163  {
164  // Compute the value of the mapping shape function i at quadrature point p
165  // (Lagrange shape functions are used for mapping)
166  if (is_linear)
167  {
168  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
169  {
170  this->phi_map[i][0] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[0]);
171  this->dphidxi_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
172 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
173  this->d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
174 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
175  for (std::size_t p=1; p<n_qp; p++)
176  {
177  this->phi_map[i][p] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
178  this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
179 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
180  this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
181 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
182  }
183  }
184  }
185  else
186  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
187  for (std::size_t p=0; p<n_qp; p++)
188  {
189  this->phi_map[i][p] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
190  this->dphidxi_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
191 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
192  this->d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
193 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
194  }
195 
196  break;
197  }
198  //------------------------------------------------------------
199  // 2D
200  case 2:
201  {
202  // Compute the value of the mapping shape function i at quadrature point p
203  // (Lagrange shape functions are used for mapping)
204  if (is_linear)
205  {
206  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
207  {
208  this->phi_map[i][0] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[0]);
209  this->dphidxi_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
210  this->dphideta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
211 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
212  this->d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
213  this->d2phidxideta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
214  this->d2phideta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
215 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
216  for (std::size_t p=1; p<n_qp; p++)
217  {
218  this->phi_map[i][p] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
219  this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
220  this->dphideta_map[i][p] = this->dphideta_map[i][0];
221 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
222  this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
223  this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
224  this->d2phideta2_map[i][p] = this->d2phideta2_map[i][0];
225 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
226  }
227  }
228  }
229  else
230  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
231  for (std::size_t p=0; p<n_qp; p++)
232  {
233  this->phi_map[i][p] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
234  this->dphidxi_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
235  this->dphideta_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
236 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
237  this->d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
238  this->d2phidxideta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
239  this->d2phideta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
240 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
241  }
242 
243  break;
244  }
245 
246  //------------------------------------------------------------
247  // 3D
248  case 3:
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  this->phi_map[i][0] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[0]);
257  this->dphidxi_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
258  this->dphideta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
259  this->dphidzeta_map[i][0] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
260 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
261  this->d2phidxi2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[0]);
262  this->d2phidxideta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[0]);
263  this->d2phideta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[0]);
264  this->d2phidxidzeta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 3, qp[0]);
265  this->d2phidetadzeta_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 4, qp[0]);
266  this->d2phidzeta2_map[i][0] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 5, qp[0]);
267 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
268  for (std::size_t p=1; p<n_qp; p++)
269  {
270  this->phi_map[i][p] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
271  this->dphidxi_map[i][p] = this->dphidxi_map[i][0];
272  this->dphideta_map[i][p] = this->dphideta_map[i][0];
273  this->dphidzeta_map[i][p] = this->dphidzeta_map[i][0];
274 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
275  this->d2phidxi2_map[i][p] = this->d2phidxi2_map[i][0];
276  this->d2phidxideta_map[i][p] = this->d2phidxideta_map[i][0];
277  this->d2phideta2_map[i][p] = this->d2phideta2_map[i][0];
278  this->d2phidxidzeta_map[i][p] = this->d2phidxidzeta_map[i][0];
279  this->d2phidetadzeta_map[i][p] = this->d2phidetadzeta_map[i][0];
280  this->d2phidzeta2_map[i][p] = this->d2phidzeta2_map[i][0];
281 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
282  }
283  }
284  }
285  else
286  for (unsigned int i=0; i<n_mapping_shape_functions; i++)
287  for (std::size_t p=0; p<n_qp; p++)
288  {
289  this->phi_map[i][p] = FE<Dim,LAGRANGE>::shape (mapping_elem_type, mapping_order, i, qp[p]);
290  this->dphidxi_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
291  this->dphideta_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
292  this->dphidzeta_map[i][p] = FE<Dim,LAGRANGE>::shape_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
293 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
294  this->d2phidxi2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 0, qp[p]);
295  this->d2phidxideta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 1, qp[p]);
296  this->d2phideta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 2, qp[p]);
297  this->d2phidxidzeta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 3, qp[p]);
298  this->d2phidetadzeta_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 4, qp[p]);
299  this->d2phidzeta2_map[i][p] = FE<Dim,LAGRANGE>::shape_second_deriv (mapping_elem_type, mapping_order, i, 5, qp[p]);
300 #endif // ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
301  }
302 
303  break;
304  }
305 
306  default:
307  libmesh_error();
308  }
309 
310  // Stop logging the reference->physical map initialization
311  STOP_LOG("init_reference_to_physical_map()", "FEMap");
312  return;
313 }
314 
315 
316 
317 void FEMap::compute_single_point_map(const unsigned int dim,
318  const std::vector<Real>& qw,
319  const Elem* elem,
320  unsigned int p)
321 {
322  libmesh_assert(elem);
323 
324  switch (dim)
325  {
326  //--------------------------------------------------------------------
327  // 0D
328  case 0:
329  {
330  xyz[p] = elem->point(0);
331  jac[p] = 1.0;
332  JxW[p] = qw[p];
333  break;
334  }
335 
336  //--------------------------------------------------------------------
337  // 1D
338  case 1:
339  {
340  // Clear the entities that will be summed
341  xyz[p].zero();
342  dxyzdxi_map[p].zero();
343 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
344  d2xyzdxi2_map[p].zero();
345 #endif
346 
347  // compute x, dx, d2x at the quadrature point
348  for (unsigned int i=0; i<phi_map.size(); i++) // sum over the nodes
349  {
350  // Reference to the point, helps eliminate
351  // exessive temporaries in the inner loop
352  const Point& elem_point = elem->point(i);
353 
354  xyz[p].add_scaled (elem_point, phi_map[i][p] );
355  dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p]);
356 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
357  d2xyzdxi2_map[p].add_scaled(elem_point, d2phidxi2_map[i][p]);
358 #endif
359  }
360 
361  // Compute the jacobian
362  //
363  // 1D elements can live in 2D or 3D space.
364  // The transformation matrix from local->global
365  // coordinates is
366  //
367  // T = | dx/dxi |
368  // | dy/dxi |
369  // | dz/dxi |
370  //
371  // The generalized determinant of T (from the
372  // so-called "normal" eqns.) is
373  // jac = "det(T)" = sqrt(det(T'T))
374  //
375  // where T'= transpose of T, so
376  //
377  // jac = sqrt( (dx/dxi)^2 + (dy/dxi)^2 + (dz/dxi)^2 )
378  jac[p] = dxyzdxi_map[p].size();
379 
380  if (jac[p] <= 0.)
381  {
382  libMesh::err << "ERROR: negative Jacobian: "
383  << jac[p]
384  << " in element "
385  << elem->id()
386  << std::endl;
387  libmesh_error();
388  }
389 
390  // The inverse Jacobian entries also come from the
391  // generalized inverse of T (see also the 2D element
392  // living in 3D code).
393  const Real jacm2 = 1./jac[p]/jac[p];
394  dxidx_map[p] = jacm2*dxdxi_map(p);
395 #if LIBMESH_DIM > 1
396  dxidy_map[p] = jacm2*dydxi_map(p);
397 #endif
398 #if LIBMESH_DIM > 2
399  dxidz_map[p] = jacm2*dzdxi_map(p);
400 #endif
401 
402  JxW[p] = jac[p]*qw[p];
403 
404  // done computing the map
405  break;
406  }
407 
408 
409  //--------------------------------------------------------------------
410  // 2D
411  case 2:
412  {
413  //------------------------------------------------------------------
414  // Compute the (x,y) values at the quadrature points,
415  // the Jacobian at the quadrature points
416 
417  xyz[p].zero();
418 
419  dxyzdxi_map[p].zero();
420  dxyzdeta_map[p].zero();
421 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
422  d2xyzdxi2_map[p].zero();
423  d2xyzdxideta_map[p].zero();
424  d2xyzdeta2_map[p].zero();
425 #endif
426 
427 
428  // compute (x,y) at the quadrature points, derivatives once
429  for (unsigned int i=0; i<phi_map.size(); i++) // sum over the nodes
430  {
431  // Reference to the point, helps eliminate
432  // exessive temporaries in the inner loop
433  const Point& elem_point = elem->point(i);
434 
435  xyz[p].add_scaled (elem_point, phi_map[i][p] );
436 
437  dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
438  dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p]);
439 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
440  d2xyzdxi2_map[p].add_scaled (elem_point, d2phidxi2_map[i][p]);
441  d2xyzdxideta_map[p].add_scaled (elem_point, d2phidxideta_map[i][p]);
442  d2xyzdeta2_map[p].add_scaled (elem_point, d2phideta2_map[i][p]);
443 #endif
444  }
445 
446  // compute the jacobian once
447  const Real dx_dxi = dxdxi_map(p),
448  dx_deta = dxdeta_map(p),
449  dy_dxi = dydxi_map(p),
450  dy_deta = dydeta_map(p);
451 
452 #if LIBMESH_DIM == 2
453  // Compute the Jacobian. This assumes the 2D face
454  // lives in 2D space
455  //
456  // Symbolically, the matrix determinant is
457  //
458  // | dx/dxi dx/deta |
459  // jac = | dy/dxi dy/deta |
460  //
461  // jac = dx/dxi*dy/deta - dx/deta*dy/dxi
462  jac[p] = (dx_dxi*dy_deta - dx_deta*dy_dxi);
463 
464  if (jac[p] <= 0.)
465  {
466  libMesh::err << "ERROR: negative Jacobian: "
467  << jac[p]
468  << " in element "
469  << elem->id()
470  << std::endl;
471  libmesh_error();
472  }
473 
474  JxW[p] = jac[p]*qw[p];
475 
476  // Compute the shape function derivatives wrt x,y at the
477  // quadrature points
478  const Real inv_jac = 1./jac[p];
479 
480  dxidx_map[p] = dy_deta*inv_jac; //dxi/dx = (1/J)*dy/deta
481  dxidy_map[p] = -dx_deta*inv_jac; //dxi/dy = -(1/J)*dx/deta
482  detadx_map[p] = -dy_dxi* inv_jac; //deta/dx = -(1/J)*dy/dxi
483  detady_map[p] = dx_dxi* inv_jac; //deta/dy = (1/J)*dx/dxi
484 
485  dxidz_map[p] = detadz_map[p] = 0.;
486 #else
487 
488  const Real dz_dxi = dzdxi_map(p),
489  dz_deta = dzdeta_map(p);
490 
491  // Compute the Jacobian. This assumes a 2D face in
492  // 3D space.
493  //
494  // The transformation matrix T from local to global
495  // coordinates is
496  //
497  // | dx/dxi dx/deta |
498  // T = | dy/dxi dy/deta |
499  // | dz/dxi dz/deta |
500  // note det(T' T) = det(T')det(T) = det(T)det(T)
501  // so det(T) = std::sqrt(det(T' T))
502  //
503  //----------------------------------------------
504  // Notes:
505  //
506  // dX = R dXi -> R'dX = R'R dXi
507  // (R^-1)dX = dXi [(R'R)^-1 R']dX = dXi
508  //
509  // so R^-1 = (R'R)^-1 R'
510  //
511  // and R^-1 R = (R'R)^-1 R'R = I.
512  //
513  const Real g11 = (dx_dxi*dx_dxi +
514  dy_dxi*dy_dxi +
515  dz_dxi*dz_dxi);
516 
517  const Real g12 = (dx_dxi*dx_deta +
518  dy_dxi*dy_deta +
519  dz_dxi*dz_deta);
520 
521  const Real g21 = g12;
522 
523  const Real g22 = (dx_deta*dx_deta +
524  dy_deta*dy_deta +
525  dz_deta*dz_deta);
526 
527  const Real det = (g11*g22 - g12*g21);
528 
529  if (det <= 0.)
530  {
531  libMesh::err << "ERROR: negative Jacobian! "
532  << " in element "
533  << elem->id()
534  << std::endl;
535  libmesh_error();
536  }
537 
538  const Real inv_det = 1./det;
539  jac[p] = std::sqrt(det);
540 
541  JxW[p] = jac[p]*qw[p];
542 
543  const Real g11inv = g22*inv_det;
544  const Real g12inv = -g12*inv_det;
545  const Real g21inv = -g21*inv_det;
546  const Real g22inv = g11*inv_det;
547 
548  dxidx_map[p] = g11inv*dx_dxi + g12inv*dx_deta;
549  dxidy_map[p] = g11inv*dy_dxi + g12inv*dy_deta;
550  dxidz_map[p] = g11inv*dz_dxi + g12inv*dz_deta;
551 
552  detadx_map[p] = g21inv*dx_dxi + g22inv*dx_deta;
553  detady_map[p] = g21inv*dy_dxi + g22inv*dy_deta;
554  detadz_map[p] = g21inv*dz_dxi + g22inv*dz_deta;
555 
556 #endif
557  // done computing the map
558  break;
559  }
560 
561 
562 
563  //--------------------------------------------------------------------
564  // 3D
565  case 3:
566  {
567  //------------------------------------------------------------------
568  // Compute the (x,y,z) values at the quadrature points,
569  // the Jacobian at the quadrature point
570 
571  // Clear the entities that will be summed
572  xyz[p].zero ();
573  dxyzdxi_map[p].zero ();
574  dxyzdeta_map[p].zero ();
575  dxyzdzeta_map[p].zero ();
576 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
577  d2xyzdxi2_map[p].zero();
578  d2xyzdxideta_map[p].zero();
579  d2xyzdxidzeta_map[p].zero();
580  d2xyzdeta2_map[p].zero();
581  d2xyzdetadzeta_map[p].zero();
582  d2xyzdzeta2_map[p].zero();
583 #endif
584 
585 
586  // compute (x,y,z) at the quadrature points,
587  // dxdxi, dydxi, dzdxi,
588  // dxdeta, dydeta, dzdeta,
589  // dxdzeta, dydzeta, dzdzeta all once
590  for (unsigned int i=0; i<phi_map.size(); i++) // sum over the nodes
591  {
592  // Reference to the point, helps eliminate
593  // exessive temporaries in the inner loop
594  const Point& elem_point = elem->point(i);
595 
596  xyz[p].add_scaled (elem_point, phi_map[i][p] );
597  dxyzdxi_map[p].add_scaled (elem_point, dphidxi_map[i][p] );
598  dxyzdeta_map[p].add_scaled (elem_point, dphideta_map[i][p] );
599  dxyzdzeta_map[p].add_scaled (elem_point, dphidzeta_map[i][p]);
600 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
601  d2xyzdxi2_map[p].add_scaled (elem_point,
602  d2phidxi2_map[i][p]);
603  d2xyzdxideta_map[p].add_scaled (elem_point,
604  d2phidxideta_map[i][p]);
605  d2xyzdxidzeta_map[p].add_scaled (elem_point,
606  d2phidxidzeta_map[i][p]);
607  d2xyzdeta2_map[p].add_scaled (elem_point,
608  d2phideta2_map[i][p]);
609  d2xyzdetadzeta_map[p].add_scaled (elem_point,
610  d2phidetadzeta_map[i][p]);
611  d2xyzdzeta2_map[p].add_scaled (elem_point,
612  d2phidzeta2_map[i][p]);
613 #endif
614  }
615 
616  // compute the jacobian
617  const Real
618  dx_dxi = dxdxi_map(p), dy_dxi = dydxi_map(p), dz_dxi = dzdxi_map(p),
619  dx_deta = dxdeta_map(p), dy_deta = dydeta_map(p), dz_deta = dzdeta_map(p),
620  dx_dzeta = dxdzeta_map(p), dy_dzeta = dydzeta_map(p), dz_dzeta = dzdzeta_map(p);
621 
622  // Symbolically, the matrix determinant is
623  //
624  // | dx/dxi dy/dxi dz/dxi |
625  // jac = | dx/deta dy/deta dz/deta |
626  // | dx/dzeta dy/dzeta dz/dzeta |
627  //
628  // jac = dx/dxi*(dy/deta*dz/dzeta - dz/deta*dy/dzeta) +
629  // dy/dxi*(dz/deta*dx/dzeta - dx/deta*dz/dzeta) +
630  // dz/dxi*(dx/deta*dy/dzeta - dy/deta*dx/dzeta)
631 
632  jac[p] = (dx_dxi*(dy_deta*dz_dzeta - dz_deta*dy_dzeta) +
633  dy_dxi*(dz_deta*dx_dzeta - dx_deta*dz_dzeta) +
634  dz_dxi*(dx_deta*dy_dzeta - dy_deta*dx_dzeta));
635 
636  if (jac[p] <= 0.)
637  {
638  libMesh::err << "ERROR: negative Jacobian: "
639  << jac[p]
640  << " in element "
641  << elem->id()
642  << std::endl;
643  libmesh_error();
644  }
645 
646  JxW[p] = jac[p]*qw[p];
647 
648  // Compute the shape function derivatives wrt x,y at the
649  // quadrature points
650  const Real inv_jac = 1./jac[p];
651 
652  dxidx_map[p] = (dy_deta*dz_dzeta - dz_deta*dy_dzeta)*inv_jac;
653  dxidy_map[p] = (dz_deta*dx_dzeta - dx_deta*dz_dzeta)*inv_jac;
654  dxidz_map[p] = (dx_deta*dy_dzeta - dy_deta*dx_dzeta)*inv_jac;
655 
656  detadx_map[p] = (dz_dxi*dy_dzeta - dy_dxi*dz_dzeta )*inv_jac;
657  detady_map[p] = (dx_dxi*dz_dzeta - dz_dxi*dx_dzeta )*inv_jac;
658  detadz_map[p] = (dy_dxi*dx_dzeta - dx_dxi*dy_dzeta )*inv_jac;
659 
660  dzetadx_map[p] = (dy_dxi*dz_deta - dz_dxi*dy_deta )*inv_jac;
661  dzetady_map[p] = (dz_dxi*dx_deta - dx_dxi*dz_deta )*inv_jac;
662  dzetadz_map[p] = (dx_dxi*dy_deta - dy_dxi*dx_deta )*inv_jac;
663 
664  // done computing the map
665  break;
666  }
667 
668  default:
669  libmesh_error();
670  }
671 }
672 
673 
674 
675 void FEMap::resize_quadrature_map_vectors(const unsigned int dim, unsigned int n_qp)
676 {
677  // Resize the vectors to hold data at the quadrature points
678  xyz.resize(n_qp);
679  dxyzdxi_map.resize(n_qp);
680  dxidx_map.resize(n_qp);
681  dxidy_map.resize(n_qp); // 1D element may live in 2D ...
682  dxidz_map.resize(n_qp); // ... or 3D
683 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
684  d2xyzdxi2_map.resize(n_qp);
685 #endif
686  if (dim > 1)
687  {
688  dxyzdeta_map.resize(n_qp);
689  detadx_map.resize(n_qp);
690  detady_map.resize(n_qp);
691  detadz_map.resize(n_qp);
692 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
693  d2xyzdxideta_map.resize(n_qp);
694  d2xyzdeta2_map.resize(n_qp);
695 #endif
696  if (dim > 2)
697  {
698  dxyzdzeta_map.resize (n_qp);
699  dzetadx_map.resize (n_qp);
700  dzetady_map.resize (n_qp);
701  dzetadz_map.resize (n_qp);
702 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
703  d2xyzdxidzeta_map.resize(n_qp);
704  d2xyzdetadzeta_map.resize(n_qp);
705  d2xyzdzeta2_map.resize(n_qp);
706 #endif
707  }
708  }
709 
710  jac.resize(n_qp);
711  JxW.resize(n_qp);
712 }
713 
714 
715 
716 void FEMap::compute_affine_map( const unsigned int dim,
717  const std::vector<Real>& qw,
718  const Elem* elem )
719 {
720  // Start logging the map computation.
721  START_LOG("compute_affine_map()", "FEMap");
722 
723  libmesh_assert(elem);
724 
725  const unsigned int n_qp = libmesh_cast_int<unsigned int>(qw.size());
726 
727  // Resize the vectors to hold data at the quadrature points
728  this->resize_quadrature_map_vectors(dim, n_qp);
729 
730  // Compute map at quadrature point 0
731  this->compute_single_point_map(dim, qw, elem, 0);
732 
733  // Compute xyz at all other quadrature points
734  for (unsigned int p=1; p<n_qp; p++)
735  {
736  xyz[p].zero();
737  for (unsigned int i=0; i<phi_map.size(); i++) // sum over the nodes
738  xyz[p].add_scaled (elem->point(i), phi_map[i][p] );
739  }
740 
741  // Copy other map data from quadrature point 0
742  for (unsigned int p=1; p<n_qp; p++) // for each extra quadrature point
743  {
744  dxyzdxi_map[p] = dxyzdxi_map[0];
745  dxidx_map[p] = dxidx_map[0];
746  dxidy_map[p] = dxidy_map[0];
747  dxidz_map[p] = dxidz_map[0];
748 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
749  // The map should be affine, so second derivatives are zero
750  d2xyzdxi2_map[p] = 0.;
751 #endif
752  if (dim > 1)
753  {
754  dxyzdeta_map[p] = dxyzdeta_map[0];
755  detadx_map[p] = detadx_map[0];
756  detady_map[p] = detady_map[0];
757  detadz_map[p] = detadz_map[0];
758 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
759  d2xyzdxideta_map[p] = 0.;
760  d2xyzdeta2_map[p] = 0.;
761 #endif
762  if (dim > 2)
763  {
765  dzetadx_map[p] = dzetadx_map[0];
766  dzetady_map[p] = dzetady_map[0];
767  dzetadz_map[p] = dzetadz_map[0];
768 #ifdef LIBMESH_ENABLE_SECOND_DERIVATIVES
769  d2xyzdxidzeta_map[p] = 0.;
770  d2xyzdetadzeta_map[p] = 0.;
771  d2xyzdzeta2_map[p] = 0.;
772 #endif
773  }
774  }
775  jac[p] = jac[0];
776  JxW[p] = JxW[0] / qw[0] * qw[p];
777  }
778 
779  STOP_LOG("compute_affine_map()", "FEMap");
780 }
781 
782 
783 
784 void FEMap::compute_map(const unsigned int dim,
785  const std::vector<Real>& qw,
786  const Elem* elem)
787 {
788  if (elem->has_affine_map())
789  {
790  compute_affine_map(dim, qw, elem);
791  return;
792  }
793 
794  // Start logging the map computation.
795  START_LOG("compute_map()", "FEMap");
796 
797  libmesh_assert(elem);
798 
799  const unsigned int n_qp = libmesh_cast_int<unsigned int>(qw.size());
800 
801  // Resize the vectors to hold data at the quadrature points
802  this->resize_quadrature_map_vectors(dim, n_qp);
803 
804  // Compute map at all quadrature points
805  for (unsigned int p=0; p!=n_qp; p++)
806  this->compute_single_point_map(dim, qw, elem, p);
807 
808  // Stop logging the map computation.
809  STOP_LOG("compute_map()", "FEMap");
810 }
811 
812 
813 
814 void FEMap::print_JxW(std::ostream& os) const
815 {
816  for (unsigned int i=0; i<JxW.size(); ++i)
817  os << " [" << i << "]: " << JxW[i] << std::endl;
818 }
819 
820 
821 
822 void FEMap::print_xyz(std::ostream& os) const
823 {
824  for (unsigned int i=0; i<xyz.size(); ++i)
825  os << " [" << i << "]: " << xyz[i];
826 }
827 
828 
829 
830 // TODO: PB: We should consider moving this to the FEMap class
831 template <unsigned int Dim, FEFamily T>
833  const Point& physical_point,
834  const Real tolerance,
835  const bool secure)
836 {
837  libmesh_assert(elem);
838  libmesh_assert_greater_equal (tolerance, 0.);
839 
840 
841  // Start logging the map inversion.
842  START_LOG("inverse_map()", "FE");
843 
844  // How much did the point on the reference
845  // element change by in this Newton step?
846  Real inverse_map_error = 0.;
847 
848  // The point on the reference element. This is
849  // the "initial guess" for Newton's method. The
850  // centroid seems like a good idea, but computing
851  // it is a little more intensive than, say taking
852  // the zero point.
853  //
854  // Convergence should be insensitive of this choice
855  // for "good" elements.
856  Point p; // the zero point. No computation required
857 
858  // The number of iterations in the map inversion process.
859  unsigned int cnt = 0;
860 
861  // The number of iterations after which we give up and declare
862  // divergence
863  const unsigned int max_cnt = 10;
864 
865  // The distance (in master element space) beyond which we give up
866  // and declare divergence. This is no longer used...
867  // Real max_step_length = 4.;
868 
869 
870 
871  // Newton iteration loop.
872  do
873  {
874  // Where our current iterate \p p maps to.
875  const Point physical_guess = FE<Dim,T>::map (elem, p);
876 
877  // How far our current iterate is from the actual point.
878  const Point delta = physical_point - physical_guess;
879 
880  // Increment in current iterate \p p, will be computed.
881  Point dp;
882 
883 
884  // The form of the map and how we invert it depends
885  // on the dimension that we are in.
886  switch (Dim)
887  {
888  // ------------------------------------------------------------------
889  // 0D map inversion is trivial
890  case 0:
891  {
892  break;
893  }
894 
895  // ------------------------------------------------------------------
896  // 1D map inversion
897  //
898  // Here we find the point on a 1D reference element that maps to
899  // the point \p physical_point in the domain. This is a bit tricky
900  // since we do not want to assume that the point \p physical_point
901  // is also in a 1D domain. In particular, this method might get
902  // called on the edge of a 3D element, in which case
903  // \p physical_point actually lives in 3D.
904  case 1:
905  {
906  const Point dxi = FE<Dim,T>::map_xi (elem, p);
907 
908  // Newton's method in this case looks like
909  //
910  // {X} - {X_n} = [J]*dp
911  //
912  // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x1 matrix
913  // d(x,y,z)/dxi, and we seek dp, a scalar. Since the above
914  // system is either overdetermined or rank-deficient, we will
915  // solve the normal equations for this system
916  //
917  // [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
918  //
919  // which involves the trivial inversion of the scalar
920  // G = [J]^T [J]
921  const Real G = dxi*dxi;
922 
923  if (secure)
924  libmesh_assert_greater (G, 0.);
925 
926  const Real Ginv = 1./G;
927 
928  const Real dxidelta = dxi*delta;
929 
930  dp(0) = Ginv*dxidelta;
931 
932  // No master elements have radius > 4, but sometimes we
933  // can take a step that big while still converging
934  // if (secure)
935  // libmesh_assert_less (dp.size(), max_step_length);
936 
937  break;
938  }
939 
940 
941 
942  // ------------------------------------------------------------------
943  // 2D map inversion
944  //
945  // Here we find the point on a 2D reference element that maps to
946  // the point \p physical_point in the domain. This is a bit tricky
947  // since we do not want to assume that the point \p physical_point
948  // is also in a 2D domain. In particular, this method might get
949  // called on the face of a 3D element, in which case
950  // \p physical_point actually lives in 3D.
951  case 2:
952  {
953  const Point dxi = FE<Dim,T>::map_xi (elem, p);
954  const Point deta = FE<Dim,T>::map_eta (elem, p);
955 
956  // Newton's method in this case looks like
957  //
958  // {X} - {X_n} = [J]*{dp}
959  //
960  // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x2 matrix
961  // d(x,y,z)/d(xi,eta), and we seek {dp}, a 2x1 vector. Since
962  // the above system is either overdermined or rank-deficient,
963  // we will solve the normal equations for this system
964  //
965  // [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
966  //
967  // which involves the inversion of the 2x2 matrix
968  // [G] = [J]^T [J]
969  const Real
970  G11 = dxi*dxi, G12 = dxi*deta,
971  G21 = dxi*deta, G22 = deta*deta;
972 
973 
974  const Real det = (G11*G22 - G12*G21);
975 
976  if (secure)
977  libmesh_assert_not_equal_to (det, 0.);
978 
979  const Real inv_det = 1./det;
980 
981  const Real
982  Ginv11 = G22*inv_det,
983  Ginv12 = -G12*inv_det,
984 
985  Ginv21 = -G21*inv_det,
986  Ginv22 = G11*inv_det;
987 
988 
989  const Real dxidelta = dxi*delta;
990  const Real detadelta = deta*delta;
991 
992  dp(0) = (Ginv11*dxidelta + Ginv12*detadelta);
993  dp(1) = (Ginv21*dxidelta + Ginv22*detadelta);
994 
995  // No master elements have radius > 4, but sometimes we
996  // can take a step that big while still converging
997  // if (secure)
998  // libmesh_assert_less (dp.size(), max_step_length);
999 
1000  break;
1001  }
1002 
1003 
1004 
1005  // ------------------------------------------------------------------
1006  // 3D map inversion
1007  //
1008  // Here we find the point in a 3D reference element that maps to
1009  // the point \p physical_point in a 3D domain. Nothing special
1010  // has to happen here, since (unless the map is singular because
1011  // you have a BAD element) the map will be invertable and we can
1012  // apply Newton's method directly.
1013  case 3:
1014  {
1015  const Point dxi = FE<Dim,T>::map_xi (elem, p);
1016  const Point deta = FE<Dim,T>::map_eta (elem, p);
1017  const Point dzeta = FE<Dim,T>::map_zeta (elem, p);
1018 
1019  // Newton's method in this case looks like
1020  //
1021  // {X} = {X_n} + [J]*{dp}
1022  //
1023  // Where {X}, {X_n} are 3x1 vectors, [J] is a 3x3 matrix
1024  // d(x,y,z)/d(xi,eta,zeta), and we seek {dp}, a 3x1 vector.
1025  // Since the above system is nonsingular for invertable maps
1026  // we will solve
1027  //
1028  // {dp} = [J]^-1 ({X} - {X_n})
1029  //
1030  // which involves the inversion of the 3x3 matrix [J]
1031  const Real
1032  J11 = dxi(0), J12 = deta(0), J13 = dzeta(0),
1033  J21 = dxi(1), J22 = deta(1), J23 = dzeta(1),
1034  J31 = dxi(2), J32 = deta(2), J33 = dzeta(2);
1035 
1036  const Real det = (J11*(J22*J33 - J23*J32) +
1037  J12*(J23*J31 - J21*J33) +
1038  J13*(J21*J32 - J22*J31));
1039 
1040  if (secure)
1041  libmesh_assert_not_equal_to (det, 0.);
1042 
1043  const Real inv_det = 1./det;
1044 
1045  const Real
1046  Jinv11 = (J22*J33 - J23*J32)*inv_det,
1047  Jinv12 = -(J12*J33 - J13*J32)*inv_det,
1048  Jinv13 = (J12*J23 - J13*J22)*inv_det,
1049 
1050  Jinv21 = -(J21*J33 - J23*J31)*inv_det,
1051  Jinv22 = (J11*J33 - J13*J31)*inv_det,
1052  Jinv23 = -(J11*J23 - J13*J21)*inv_det,
1053 
1054  Jinv31 = (J21*J32 - J22*J31)*inv_det,
1055  Jinv32 = -(J11*J32 - J12*J31)*inv_det,
1056  Jinv33 = (J11*J22 - J12*J21)*inv_det;
1057 
1058  dp(0) = (Jinv11*delta(0) +
1059  Jinv12*delta(1) +
1060  Jinv13*delta(2));
1061 
1062  dp(1) = (Jinv21*delta(0) +
1063  Jinv22*delta(1) +
1064  Jinv23*delta(2));
1065 
1066  dp(2) = (Jinv31*delta(0) +
1067  Jinv32*delta(1) +
1068  Jinv33*delta(2));
1069 
1070  // No master elements have radius > 4, but sometimes we
1071  // can take a step that big while still converging
1072  // if (secure)
1073  // libmesh_assert_less (dp.size(), max_step_length);
1074 
1075  break;
1076  }
1077 
1078 
1079  // Some other dimension?
1080  default:
1081  libmesh_error();
1082  } // end switch(Dim), dp now computed
1083 
1084 
1085 
1086  // ||P_n+1 - P_n||
1087  inverse_map_error = dp.size();
1088 
1089  // P_n+1 = P_n + dp
1090  p.add (dp);
1091 
1092  // Increment the iteration count.
1093  cnt++;
1094 
1095  // Watch for divergence of Newton's
1096  // method. Here's how it goes:
1097  // (1) For good elements, we expect convergence in 10
1098  // iterations, with no too-large steps.
1099  // - If called with (secure == true) and we have not yet converged
1100  // print out a warning message.
1101  // - If called with (secure == true) and we have not converged in
1102  // 20 iterations abort
1103  // (2) This method may be called in cases when the target point is not
1104  // inside the element and we have no business expecting convergence.
1105  // For these cases if we have not converged in 10 iterations forget
1106  // about it.
1107  if (cnt > max_cnt)
1108  {
1109  // Warn about divergence when secure is true - this
1110  // shouldn't happen
1111  if (secure)
1112  {
1113  // Print every time in devel/dbg modes
1114 #ifndef NDEBUG
1115  libmesh_here();
1116  libMesh::err << "WARNING: Newton scheme has not converged in "
1117  << cnt << " iterations:" << std::endl
1118  << " physical_point="
1119  << physical_point
1120  << " physical_guess="
1121  << physical_guess
1122  << " dp="
1123  << dp
1124  << " p="
1125  << p
1126  << " error=" << inverse_map_error
1127  << " in element " << elem->id()
1128  << std::endl;
1129 
1130  elem->print_info(libMesh::err);
1131 #else
1132  // In optimized mode, just print once that an inverse_map() call
1133  // had trouble converging its Newton iteration.
1134  libmesh_do_once(libMesh::err << "WARNING: At least one element took more than "
1135  << max_cnt
1136  << " iterations to converge in inverse_map()...\n"
1137  << "Rerun in devel/dbg mode for more details."
1138  << std::endl;);
1139 
1140 #endif // NDEBUG
1141 
1142  if (cnt > 2*max_cnt)
1143  {
1144  libMesh::err << "ERROR: Newton scheme FAILED to converge in "
1145  << cnt
1146  << " iterations!"
1147  << " in element " << elem->id()
1148  << std::endl;
1149 
1150  elem->print_info(libMesh::err);
1151 
1152  libmesh_error();
1153  }
1154  }
1155  // Return a far off point when secure is false - this
1156  // should only happen when we're trying to map a point
1157  // that's outside the element
1158  else
1159  {
1160  for (unsigned int i=0; i != Dim; ++i)
1161  p(i) = 1e6;
1162 
1163  STOP_LOG("inverse_map()", "FE");
1164  return p;
1165  }
1166  }
1167  }
1168  while (inverse_map_error > tolerance);
1169 
1170 
1171 
1172  // If we are in debug mode do two sanity checks.
1173 #ifdef DEBUG
1174 
1175  if (secure)
1176  {
1177  // Make sure the point \p p on the reference element actually
1178  // does map to the point \p physical_point within a tolerance.
1179 
1180  const Point check = FE<Dim,T>::map (elem, p);
1181  const Point diff = physical_point - check;
1182 
1183  if (diff.size() > tolerance)
1184  {
1185  libmesh_here();
1186  libMesh::err << "WARNING: diff is "
1187  << diff.size()
1188  << std::endl
1189  << " point="
1190  << physical_point;
1191  libMesh::err << " local=" << check;
1192  libMesh::err << " lref= " << p;
1193 
1194  elem->print_info(libMesh::err);
1195  }
1196 
1197  // Make sure the point \p p on the reference element actually
1198  // is
1199 
1200  if (!FEAbstract::on_reference_element(p, elem->type(), 2*tolerance))
1201  {
1202  libmesh_here();
1203  libMesh::err << "WARNING: inverse_map of physical point "
1204  << physical_point
1205  << "is not on element." << '\n';
1206  elem->print_info(libMesh::err);
1207  }
1208  }
1209 
1210 #endif
1211 
1212 
1213 
1214  // Stop logging the map inversion.
1215  STOP_LOG("inverse_map()", "FE");
1216 
1217  return p;
1218 }
1219 
1220 
1221 
1222 // TODO: PB: We should consider moving this to the FEMap class
1223 template <unsigned int Dim, FEFamily T>
1224 void FE<Dim,T>::inverse_map (const Elem* elem,
1225  const std::vector<Point>& physical_points,
1226  std::vector<Point>& reference_points,
1227  const Real tolerance,
1228  const bool secure)
1229 {
1230  // The number of points to find the
1231  // inverse map of
1232  const std::size_t n_points = physical_points.size();
1233 
1234  // Resize the vector to hold the points
1235  // on the reference element
1236  reference_points.resize(n_points);
1237 
1238  // Find the coordinates on the reference
1239  // element of each point in physical space
1240  for (std::size_t p=0; p<n_points; p++)
1241  reference_points[p] =
1242  FE<Dim,T>::inverse_map (elem, physical_points[p], tolerance, secure);
1243 }
1244 
1245 
1246 
1247 // TODO: PB: We should consider moving this to the FEMap class
1248 template <unsigned int Dim, FEFamily T>
1250  const Point& reference_point)
1251 {
1252  libmesh_assert(elem);
1253 
1254  Point p;
1255 
1256  const ElemType type = elem->type();
1257  const Order order = elem->default_order();
1258  const unsigned int n_sf = FE<Dim,LAGRANGE>::n_shape_functions(type, order);
1259 
1260  // Lagrange basis functions are used for mapping
1261  for (unsigned int i=0; i<n_sf; i++)
1262  p.add_scaled (elem->point(i),
1264  order,
1265  i,
1266  reference_point)
1267  );
1268 
1269  return p;
1270 }
1271 
1272 
1273 
1274 // TODO: PB: We should consider moving this to the FEMap class
1275 template <unsigned int Dim, FEFamily T>
1277  const Point& reference_point)
1278 {
1279  libmesh_assert(elem);
1280 
1281  Point p;
1282 
1283  const ElemType type = elem->type();
1284  const Order order = elem->default_order();
1285  const unsigned int n_sf = FE<Dim,LAGRANGE>::n_shape_functions(type, order);
1286 
1287  // Lagrange basis functions are used for mapping
1288  for (unsigned int i=0; i<n_sf; i++)
1289  p.add_scaled (elem->point(i),
1291  order,
1292  i,
1293  0,
1294  reference_point)
1295  );
1296 
1297  return p;
1298 }
1299 
1300 
1301 
1302 // TODO: PB: We should consider moving this to the FEMap class
1303 template <unsigned int Dim, FEFamily T>
1305  const Point& reference_point)
1306 {
1307  libmesh_assert(elem);
1308 
1309  Point p;
1310 
1311  const ElemType type = elem->type();
1312  const Order order = elem->default_order();
1313  const unsigned int n_sf = FE<Dim,LAGRANGE>::n_shape_functions(type, order);
1314 
1315  // Lagrange basis functions are used for mapping
1316  for (unsigned int i=0; i<n_sf; i++)
1317  p.add_scaled (elem->point(i),
1319  order,
1320  i,
1321  1,
1322  reference_point)
1323  );
1324 
1325  return p;
1326 }
1327 
1328 
1329 
1330 // TODO: PB: We should consider moving this to the FEMap class
1331 template <unsigned int Dim, FEFamily T>
1333  const Point& reference_point)
1334 {
1335  libmesh_assert(elem);
1336 
1337  Point p;
1338 
1339  const ElemType type = elem->type();
1340  const Order order = elem->default_order();
1341  const unsigned int n_sf = FE<Dim,LAGRANGE>::n_shape_functions(type, order);
1342 
1343  // Lagrange basis functions are used for mapping
1344  for (unsigned int i=0; i<n_sf; i++)
1345  p.add_scaled (elem->point(i),
1347  order,
1348  i,
1349  2,
1350  reference_point)
1351  );
1352 
1353  return p;
1354 }
1355 
1356 
1357 
1358 // Explicit instantiation of FEMap member functions
1359 template void FEMap::init_reference_to_physical_map<0>( const std::vector<Point>&, const Elem*);
1360 template void FEMap::init_reference_to_physical_map<1>( const std::vector<Point>&, const Elem*);
1361 template void FEMap::init_reference_to_physical_map<2>( const std::vector<Point>&, const Elem*);
1362 template void FEMap::init_reference_to_physical_map<3>( const std::vector<Point>&, const Elem*);
1363 
1364 //--------------------------------------------------------------
1365 // Explicit instantiations using the macro from fe_macro.h
1370 
1371 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:05 UTC

Hosted By:
SourceForge.net Logo