libMesh
inf_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 // Local includes
21 #include "libmesh/libmesh_config.h"
22 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
23 #include "libmesh/inf_fe.h"
24 #include "libmesh/fe.h"
25 #include "libmesh/elem.h"
26 #include "libmesh/inf_fe_macro.h"
27 #include "libmesh/libmesh_logging.h"
28 
29 namespace libMesh
30 {
31 
32 
33 
34 // ------------------------------------------------------------
35 // InfFE static class members concerned with coordinate
36 // mapping
37 
38 
39 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
41  const Point & reference_point)
42 {
43  libmesh_assert(inf_elem);
44  libmesh_assert_not_equal_to (Dim, 0);
45 
46  UniquePtr<Elem> base_elem (Base::build_elem (inf_elem));
47 
48  const Order radial_mapping_order (Radial::mapping_order());
49  const Real v (reference_point(Dim-1));
50 
51  // map in the base face
52  Point base_point;
53  switch (Dim)
54  {
55  case 1:
56  base_point = inf_elem->point(0);
57  break;
58  case 2:
59  base_point = FE<1,LAGRANGE>::map (base_elem.get(), reference_point);
60  break;
61  case 3:
62  base_point = FE<2,LAGRANGE>::map (base_elem.get(), reference_point);
63  break;
64 #ifdef DEBUG
65  default:
66  libmesh_error_msg("Unknown Dim = " << Dim);
67 #endif
68  }
69 
70 
71  // map in the outer node face not necessary. Simply
72  // compute the outer_point = base_point + (base_point-origin)
73  const Point outer_point (base_point * 2. - inf_elem->origin());
74 
75  Point p;
76 
77  // there are only two mapping shapes in radial direction
78  p.add_scaled (base_point, InfFE<Dim,INFINITE_MAP,T_map>::eval (v, radial_mapping_order, 0));
79  p.add_scaled (outer_point, InfFE<Dim,INFINITE_MAP,T_map>::eval (v, radial_mapping_order, 1));
80 
81  return p;
82 }
83 
84 
85 
86 
87 
88 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
90  const Point & physical_point,
91  const Real tolerance,
92  const bool secure)
93 {
94  libmesh_assert(inf_elem);
95  libmesh_assert_greater_equal (tolerance, 0.);
96 
97  // Start logging the map inversion.
98  LOG_SCOPE("inverse_map()", "InfFE");
99 
100  // 1.)
101  // build a base element to do the map inversion in the base face
102  UniquePtr<Elem> base_elem (Base::build_elem (inf_elem));
103 
104  // 2.)
105  // just like in FE<Dim-1,LAGRANGE>::inverse_map(): compute
106  // the local coordinates, but only in the base element.
107  // The radial part can then be computed directly later on.
108 
109  // How much did the point on the reference
110  // element change by in this Newton step?
111  Real inverse_map_error = 0.;
112 
113  // The point on the reference element. This is
114  // the "initial guess" for Newton's method. The
115  // centroid seems like a good idea, but computing
116  // it is a little more intensive than, say taking
117  // the zero point.
118  //
119  // Convergence should be insensitive of this choice
120  // for "good" elements.
121  Point p; // the zero point. No computation required
122 
123  // Now find the intersection of a plane represented by the base
124  // element nodes and the line given by the origin of the infinite
125  // element and the physical point.
126  Point intersection;
127 
128  // the origin of the infinite element
129  const Point o = inf_elem->origin();
130 
131  switch (Dim)
132  {
133  // unnecessary for 1D
134  case 1:
135  break;
136 
137  case 2:
138  libmesh_error_msg("ERROR: InfFE::inverse_map is not yet implemented in 2d");
139 
140  case 3:
141  {
142  // references to the nodal points of the base element
143  const Point & p0 = base_elem->point(0);
144  const Point & p1 = base_elem->point(1);
145  const Point & p2 = base_elem->point(2);
146 
147  // a reference to the physical point
148  const Point & fp = physical_point;
149 
150  // The intersection of the plane and the line is given by
151  // can be computed solving a linear 3x3 system
152  // a*({p1}-{p0})+b*({p2}-{p0})-c*({fp}-{o})={fp}-{p0}.
153  const Real c_factor = -(p1(0)*fp(1)*p0(2)-p1(0)*fp(2)*p0(1)
154  +fp(0)*p1(2)*p0(1)-p0(0)*fp(1)*p1(2)
155  +p0(0)*fp(2)*p1(1)+p2(0)*fp(2)*p0(1)
156  -p2(0)*fp(1)*p0(2)-fp(0)*p2(2)*p0(1)
157  +fp(0)*p0(2)*p2(1)+p0(0)*fp(1)*p2(2)
158  -p0(0)*fp(2)*p2(1)-fp(0)*p0(2)*p1(1)
159  +p0(2)*p2(0)*p1(1)-p0(1)*p2(0)*p1(2)
160  -fp(0)*p1(2)*p2(1)+p2(1)*p0(0)*p1(2)
161  -p2(0)*fp(2)*p1(1)-p1(0)*fp(1)*p2(2)
162  +p2(2)*p1(0)*p0(1)+p1(0)*fp(2)*p2(1)
163  -p0(2)*p1(0)*p2(1)-p2(2)*p0(0)*p1(1)
164  +fp(0)*p2(2)*p1(1)+p2(0)*fp(1)*p1(2))/
165  (fp(0)*p1(2)*p0(1)-p1(0)*fp(2)*p0(1)
166  +p1(0)*fp(1)*p0(2)-p1(0)*o(1)*p0(2)
167  +o(0)*p2(2)*p0(1)-p0(0)*fp(2)*p2(1)
168  +p1(0)*o(1)*p2(2)+fp(0)*p0(2)*p2(1)
169  -fp(0)*p1(2)*p2(1)-p0(0)*o(1)*p2(2)
170  +p0(0)*fp(1)*p2(2)-o(0)*p0(2)*p2(1)
171  +o(0)*p1(2)*p2(1)-p2(0)*fp(2)*p1(1)
172  +fp(0)*p2(2)*p1(1)-p2(0)*fp(1)*p0(2)
173  -o(2)*p0(0)*p1(1)-fp(0)*p0(2)*p1(1)
174  +p0(0)*o(1)*p1(2)+p0(0)*fp(2)*p1(1)
175  -p0(0)*fp(1)*p1(2)-o(0)*p1(2)*p0(1)
176  -p2(0)*o(1)*p1(2)-o(2)*p2(0)*p0(1)
177  -o(2)*p1(0)*p2(1)+o(2)*p0(0)*p2(1)
178  -fp(0)*p2(2)*p0(1)+o(2)*p2(0)*p1(1)
179  +p2(0)*o(1)*p0(2)+p2(0)*fp(1)*p1(2)
180  +p2(0)*fp(2)*p0(1)-p1(0)*fp(1)*p2(2)
181  +p1(0)*fp(2)*p2(1)-o(0)*p2(2)*p1(1)
182  +o(2)*p1(0)*p0(1)+o(0)*p0(2)*p1(1));
183 
184 
185  // Check whether the above system is ill-posed. It should
186  // only happen when \p physical_point is not in \p
187  // inf_elem. In this case, any point that is not in the
188  // element is a valid answer.
189  if (libmesh_isinf(c_factor))
190  return o;
191 
192  // Compute the intersection with
193  // {intersection} = {fp} + c*({fp}-{o}).
194  intersection.add_scaled(fp,1.+c_factor);
195  intersection.add_scaled(o,-c_factor);
196 
197  break;
198  }
199 
200  default:
201  libmesh_error_msg("Invalid dim = " << Dim);
202  }
203 
204  // The number of iterations in the map inversion process.
205  unsigned int cnt = 0;
206 
207  // Newton iteration loop.
208  do
209  {
210  // Increment in current iterate \p p, will be computed.
211  // Automatically initialized to all zero. Note that
212  // in 3D, actually only the first two entries are
213  // filled by the inverse map, and in 2D only the first.
214  Point dp;
215 
216  // The form of the map and how we invert it depends
217  // on the dimension that we are in.
218  switch (Dim)
219  {
220  // 1D infinite element - no map inversion necessary
221  case 1:
222  break;
223 
224  // 2D infinite element - 1D map inversion
225  //
226  // In this iteration scheme only search for the local coordinate
227  // in xi direction. Once xi is determined, the radial coordinate eta is
228  // uniquely determined, and there is no need to iterate in that direction.
229  case 2:
230  {
231  // Where our current iterate \p p maps to.
232  const Point physical_guess = FE<1,LAGRANGE>::map (base_elem.get(), p);
233 
234  // How far our current iterate is from the actual point.
235  const Point delta = physical_point - physical_guess;
236 
237  const Point dxi = FE<1,LAGRANGE>::map_xi (base_elem.get(), p);
238 
239  // For details on Newton's method see fe_map.C
240  const Real G = dxi*dxi;
241 
242  if (secure)
243  libmesh_assert_greater (G, 0.);
244 
245  const Real Ginv = 1./G;
246 
247  const Real dxidelta = dxi*delta;
248 
249  // compute only the first coordinate
250  dp(0) = Ginv*dxidelta;
251 
252  break;
253  }
254 
255  // 3D infinite element - 2D map inversion
256  //
257  // In this iteration scheme only search for the local coordinates
258  // in xi and eta direction. Once xi, eta are determined, the radial
259  // coordinate zeta may directly computed.
260  case 3:
261  {
262  // Where our current iterate \p p maps to.
263  const Point physical_guess = FE<2,LAGRANGE>::map (base_elem.get(), p);
264 
265  // How far our current iterate is from the actual point.
266  // const Point delta = physical_point - physical_guess;
267  const Point delta = intersection - physical_guess;
268 
269  const Point dxi = FE<2,LAGRANGE>::map_xi (base_elem.get(), p);
270  const Point deta = FE<2,LAGRANGE>::map_eta (base_elem.get(), p);
271 
272  // For details on Newton's method see fe_map.C
273  const Real
274  G11 = dxi*dxi, G12 = dxi*deta,
275  G21 = dxi*deta, G22 = deta*deta;
276 
277  const Real det = (G11*G22 - G12*G21);
278 
279  if (secure)
280  {
281  libmesh_assert_greater (det, 0.);
282  libmesh_assert_greater (std::abs(det), 1.e-10);
283  }
284 
285  const Real inv_det = 1./det;
286 
287  const Real
288  Ginv11 = G22*inv_det,
289  Ginv12 = -G12*inv_det,
290 
291  Ginv21 = -G21*inv_det,
292  Ginv22 = G11*inv_det;
293 
294  const Real dxidelta = dxi*delta;
295  const Real detadelta = deta*delta;
296 
297  // compute only the first two coordinates.
298  dp(0) = (Ginv11*dxidelta + Ginv12*detadelta);
299  dp(1) = (Ginv21*dxidelta + Ginv22*detadelta);
300 
301  break;
302  }
303 
304  // Some other dimension?
305  default:
306  libmesh_error_msg("Unknown Dim = " << Dim);
307  } // end switch(Dim), dp now computed
308 
309  // determine the error in computing the local coordinates
310  // in the base: ||P_n+1 - P_n||
311  inverse_map_error = dp.norm();
312 
313  // P_n+1 = P_n + dp
314  p.add (dp);
315 
316  // Increment the iteration count.
317  cnt++;
318 
319  // Watch for divergence of Newton's
320  // method.
321  if (cnt > 10)
322  {
323  if (secure || !secure)
324  {
325  libmesh_here();
326  {
327  libMesh::err << "WARNING: Newton scheme has not converged in "
328  << cnt << " iterations:\n"
329  << " physical_point="
330  << physical_point
331  << " dp="
332  << dp
333  << " p="
334  << p
335  << " error=" << inverse_map_error
336  << std::endl;
337  }
338  }
339 
340  if (cnt > 20)
341  libmesh_error_msg("ERROR: Newton scheme FAILED to converge in " << cnt << " iterations!");
342  }
343  }
344  while (inverse_map_error > tolerance);
345 
346  // 4.
347  //
348  // Now that we have the local coordinates in the base,
349  // we compute the radial distance with Newton iteration.
350 
351  // distance from the physical point to the ifem origin
352  const Real fp_o_dist = Point(o-physical_point).norm();
353 
354  // the distance from the intersection on the
355  // base to the origin
356  const Real a_dist = Point(o-intersection).norm();
357 
358  // element coordinate in radial direction
359  // here our first guess is 0.
360  Real v = 0.;
361 
362  // We know the analytic answer for CARTESIAN,
363  // other schemes do not yet have a better guess.
364  if (T_map == CARTESIAN)
365  v = 1.-2.*a_dist/fp_o_dist;
366 
367  // the order of the radial mapping
368  const Order radial_mapping_order (Radial::mapping_order());
369 
370  unsigned int cnt2 = 0;
371  inverse_map_error = 0.;
372 
373  // Newton iteration in 1-D
374  do
375  {
376  if (v < -1.)
377  {
378  // in this case, physical_point is not in this element.
379  // We therefore give back the best approximation:
380  // TODO: is there a strict argument for this?
381  p(Dim-1) = -1;
382  return p;
383  }
384 
385  // the mapping in radial direction
386  // note that we only have two mapping functions in
387  // radial direction
388  const Real r = a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval (v, radial_mapping_order, 0)
389  + 2. * a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval (v, radial_mapping_order, 1);
390 
391  const Real dr = a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (v, radial_mapping_order, 0)
392  + 2. * a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (v, radial_mapping_order, 1);
393 
394  const Real G = dr*dr;
395  const Real Ginv = 1./G;
396 
397  const Real delta = fp_o_dist - r;
398  const Real drdelta = dr*delta;
399 
400  Real dp = Ginv*drdelta;
401 
402  // update the radial coordinate
403  v += dp;
404 
405  // note that v should be smaller than 1,
406  // since radial mapping function tends to infinity
407  if (v >= 1.)
408  v = .9999;
409 
410  inverse_map_error = std::abs(dp);
411 
412  // increment iteration count
413  cnt2 ++;
414  if (cnt2 > 20)
415  libmesh_error_msg("ERROR: 1D Newton scheme FAILED to converge");
416  }
417  while (inverse_map_error > tolerance);
418 
419  // Set the last coordinate of the Point to v.
420  p(Dim-1) = v;
421 
422  // If we are in debug mode do a sanity check. Make sure
423  // the point \p p on the reference element actually does
424  // map to the point \p physical_point within a tolerance.
425 #ifdef DEBUG
426  const Point check = InfFE<Dim,T_radial,T_map>::map (inf_elem, p);
427  const Point diff = physical_point - check;
428 
429  if (diff.norm() > tolerance)
430  libmesh_warning("WARNING: diff is " << diff.norm());
431 #endif
432 
433  return p;
434 }
435 
436 
437 
438 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
440  const std::vector<Point> & physical_points,
441  std::vector<Point> & reference_points,
442  const Real tolerance,
443  const bool secure)
444 {
445  // The number of points to find the
446  // inverse map of
447  const std::size_t n_points = physical_points.size();
448 
449  // Resize the vector to hold the points
450  // on the reference element
451  reference_points.resize(n_points);
452 
453  // Find the coordinates on the reference
454  // element of each point in physical space
455  for (unsigned int p=0; p<n_points; p++)
456  reference_points[p] =
457  InfFE<Dim,T_radial,T_map>::inverse_map (elem, physical_points[p], tolerance, secure);
458 }
459 
460 
461 
462 
463 //--------------------------------------------------------------
464 // Explicit instantiations using the macro from inf_fe_macro.h
465 //INSTANTIATE_INF_FE(1,CARTESIAN);
466 
467 //INSTANTIATE_INF_FE(2,CARTESIAN);
468 
469 //INSTANTIATE_INF_FE(3,CARTESIAN);
470 
471 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Point, map(const Elem *, const Point &));
472 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Point, map(const Elem *, const Point &));
473 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Point, map(const Elem *, const Point &));
474 
475 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Point, inverse_map(const Elem *, const Point &, const Real, const bool));
476 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, Point, inverse_map(const Elem *, const Point &, const Real, const bool));
477 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, Point, inverse_map(const Elem *, const Point &, const Real, const bool));
478 
479 INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, void, inverse_map(const Elem *, const std::vector<Point> &, std::vector<Point> &, const Real, const bool));
480 INSTANTIATE_INF_FE_MBRF(2, CARTESIAN, void, inverse_map(const Elem *, const std::vector<Point> &, std::vector<Point> &, const Real, const bool));
481 INSTANTIATE_INF_FE_MBRF(3, CARTESIAN, void, inverse_map(const Elem *, const std::vector<Point> &, std::vector<Point> &, const Real, const bool));
482 
483 
484 } // namespace libMesh
485 
486 
487 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
OStreamProxy err
double abs(double a)
A specific instantiation of the FEBase class.
Definition: fe.h:40
static Real eval_deriv(Real v, Order o_radial, unsigned int i)
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
This is the base class from which all geometric element types are derived.
Definition: elem.h:89
Real norm() const
Definition: type_vector.h:909
INSTANTIATE_INF_FE_MBRF(1, CARTESIAN, Elem *, Base::build_elem(const Elem *))
The libMesh namespace provides an interface to certain functionality in the library.
bool libmesh_isinf(T x)
void add(const TypeVector< T2 > &)
Add to this vector without creating a temporary.
Definition: type_vector.h:600
libmesh_assert(j)
std::unique_ptr< T > UniquePtr
Definition: auto_ptr.h:46
static Point map(const Elem *elem, const Point &reference_point)
Definition: fe_map.C:1975
static Real eval(Real v, Order o_radial, unsigned int i)
virtual Point origin() const
Definition: elem.h:1480
static Point inverse_map(const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
Definition: inf_fe_map.C:89
static Point map_eta(const Elem *elem, const Point &reference_point)
Definition: fe_map.C:2030
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
const Point & point(const unsigned int i) const
Definition: elem.h:1809
static Point map(const Elem *inf_elem, const Point &reference_point)
Definition: inf_fe_map.C:40
Order
defines an enum for polynomial orders.
Definition: enum_order.h:32
Real size() const
Definition: type_vector.h:898
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:38