inf_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 // 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"
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  AutoPtr<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();
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  const bool interpolated)
94 {
95  libmesh_assert(inf_elem);
96  libmesh_assert_greater_equal (tolerance, 0.);
97 
98 
99  // Start logging the map inversion.
100  START_LOG("inverse_map()", "InfFE");
101 
102 
103  // 1.)
104  // build a base element to do the map inversion in the base face
105  AutoPtr<Elem> base_elem (Base::build_elem (inf_elem));
106 
107  const Order base_mapping_order (base_elem->default_order());
108  const ElemType base_mapping_elem_type (base_elem->type());
109  const unsigned int n_base_mapping_sf = Base::n_base_mapping_sf (base_mapping_elem_type,
110  base_mapping_order);
111 
112  const ElemType inf_elem_type = inf_elem->type();
113  if (inf_elem_type != INFHEX8 &&
114  inf_elem_type != INFPRISM6)
115  {
116  libMesh::err << "ERROR: InfFE::inverse_map is currently implemented only for\n"
117  << " infinite elments of type InfHex8 and InfPrism6." << std::endl;
118  libmesh_error();
119  }
120 
121 
122  // 2.)
123  // just like in FE<Dim-1,LAGRANGE>::inverse_map(): compute
124  // the local coordinates, but only in the base element.
125  // The radial part can then be computed directly later on.
126 
127 
128  // How much did the point on the reference
129  // element change by in this Newton step?
130  Real inverse_map_error = 0.;
131 
132 
133  // The point on the reference element. This is
134  // the "initial guess" for Newton's method. The
135  // centroid seems like a good idea, but computing
136  // it is a little more intensive than, say taking
137  // the zero point.
138  //
139  // Convergence should be insensitive of this choice
140  // for "good" elements.
141  Point p; // the zero point. No computation required
142 
143 
144 
145  // Now find the intersection of a plane represented by the base
146  // element nodes and the line given by the origin of the infinite
147  // element and the physical point.
148  Point intersection;
149 
150  // the origin of the infinite lement
151  const Point o = inf_elem->origin();
152 
153  switch (Dim)
154  {
155 
156  // unnecessary for 1D
157  case 1:
158  {
159  break;
160  }
161 
162  case 2:
163  {
164  libMesh::err << "ERROR: InfFE::inverse_map is not yet implemented"
165  << " in 2d" << std::endl;
166  libmesh_error();
167  break;
168  }
169 
170 
171  case 3:
172  {
173  // references to the nodal points of the base element
174  const Point& p0 = base_elem->point(0);
175  const Point& p1 = base_elem->point(1);
176  const Point& p2 = base_elem->point(2);
177 
178  // a reference to the physical point
179  const Point& fp = physical_point;
180 
181  // The intersection of the plane and the line is given by
182  // can be computed solving a linear 3x3 system
183  // a*({p1}-{p0})+b*({p2}-{p0})-c*({fp}-{o})={fp}-{p0}.
184  const Real c_factor = -(p1(0)*fp(1)*p0(2)-p1(0)*fp(2)*p0(1)
185  +fp(0)*p1(2)*p0(1)-p0(0)*fp(1)*p1(2)
186  +p0(0)*fp(2)*p1(1)+p2(0)*fp(2)*p0(1)
187  -p2(0)*fp(1)*p0(2)-fp(0)*p2(2)*p0(1)
188  +fp(0)*p0(2)*p2(1)+p0(0)*fp(1)*p2(2)
189  -p0(0)*fp(2)*p2(1)-fp(0)*p0(2)*p1(1)
190  +p0(2)*p2(0)*p1(1)-p0(1)*p2(0)*p1(2)
191  -fp(0)*p1(2)*p2(1)+p2(1)*p0(0)*p1(2)
192  -p2(0)*fp(2)*p1(1)-p1(0)*fp(1)*p2(2)
193  +p2(2)*p1(0)*p0(1)+p1(0)*fp(2)*p2(1)
194  -p0(2)*p1(0)*p2(1)-p2(2)*p0(0)*p1(1)
195  +fp(0)*p2(2)*p1(1)+p2(0)*fp(1)*p1(2))/
196  (fp(0)*p1(2)*p0(1)-p1(0)*fp(2)*p0(1)
197  +p1(0)*fp(1)*p0(2)-p1(0)*o(1)*p0(2)
198  +o(0)*p2(2)*p0(1)-p0(0)*fp(2)*p2(1)
199  +p1(0)*o(1)*p2(2)+fp(0)*p0(2)*p2(1)
200  -fp(0)*p1(2)*p2(1)-p0(0)*o(1)*p2(2)
201  +p0(0)*fp(1)*p2(2)-o(0)*p0(2)*p2(1)
202  +o(0)*p1(2)*p2(1)-p2(0)*fp(2)*p1(1)
203  +fp(0)*p2(2)*p1(1)-p2(0)*fp(1)*p0(2)
204  -o(2)*p0(0)*p1(1)-fp(0)*p0(2)*p1(1)
205  +p0(0)*o(1)*p1(2)+p0(0)*fp(2)*p1(1)
206  -p0(0)*fp(1)*p1(2)-o(0)*p1(2)*p0(1)
207  -p2(0)*o(1)*p1(2)-o(2)*p2(0)*p0(1)
208  -o(2)*p1(0)*p2(1)+o(2)*p0(0)*p2(1)
209  -fp(0)*p2(2)*p0(1)+o(2)*p2(0)*p1(1)
210  +p2(0)*o(1)*p0(2)+p2(0)*fp(1)*p1(2)
211  +p2(0)*fp(2)*p0(1)-p1(0)*fp(1)*p2(2)
212  +p1(0)*fp(2)*p2(1)-o(0)*p2(2)*p1(1)
213  +o(2)*p1(0)*p0(1)+o(0)*p0(2)*p1(1));
214 
215 
216  // Compute the intersection with
217  // {intersection} = {fp} + c*({fp}-{o}).
218  intersection.add_scaled(fp,1.+c_factor);
219  intersection.add_scaled(o,-c_factor);
220 
221  break;
222  }
223  }
224 
228  unsigned int cnt = 0;
229 
230 
234  do
235  {
236 
237  // Increment in current iterate \p p, will be computed.
238  // Automatically initialized to all zero. Note that
239  // in 3D, actually only the first two entries are
240  // filled by the inverse map, and in 2D only the first.
241  Point dp;
242 
243  // The form of the map and how we invert it depends
244  // on the dimension that we are in.
245  switch (Dim)
246  {
247 
248  //------------------------------------------------------------------
249  // 1D infinite element - no map inversion necessary
250  case 1:
251  {
252  break;
253  }
254 
255 
256  //------------------------------------------------------------------
257  // 2D infinite element - 1D map inversion
258  //
259  // In this iteration scheme only search for the local coordinate
260  // in xi direction. Once xi is determined, the radial coordinate eta is
261  // uniquely determined, and there is no need to iterate in that direction.
262  case 2:
263  {
264 
265  // Where our current iterate \p p maps to.
266  const Point physical_guess = FE<1,LAGRANGE>::map (base_elem.get(), p);
267 
268 
269  // How far our current iterate is from the actual point.
270  const Point delta = physical_point - physical_guess;
271 
272 
273  const Point dxi = FE<1,LAGRANGE>::map_xi (base_elem.get(), p);
274 
275 
276  // For details on Newton's method see fe_map.C
277  const Real G = dxi*dxi;
278 
279  if (secure)
280  libmesh_assert_greater (G, 0.);
281 
282  const Real Ginv = 1./G;
283 
284  const Real dxidelta = dxi*delta;
285 
286  // compute only the first coordinate
287  dp(0) = Ginv*dxidelta;
288 
289  break;
290  }
291 
292 
293 
294  //------------------------------------------------------------------
295  // 3D infinite element - 2D map inversion
296  //
297  // In this iteration scheme only search for the local coordinates
298  // in xi and eta direction. Once xi, eta are determined, the radial
299  // coordinate zeta may directly computed.
300  case 3:
301  {
302 
303 
304  // Where our current iterate \p p maps to.
305  const Point physical_guess = FE<2,LAGRANGE>::map (base_elem.get(), p);
306 
307  // How far our current iterate is from the actual point.
308  // const Point delta = physical_point - physical_guess;
309  const Point delta = intersection - physical_guess;
310 
311  const Point dxi = FE<2,LAGRANGE>::map_xi (base_elem.get(), p);
312  const Point deta = FE<2,LAGRANGE>::map_eta (base_elem.get(), p);
313 
314 
315  // For details on Newton's method see fe_map.C
316  const Real
317  G11 = dxi*dxi, G12 = dxi*deta,
318  G21 = dxi*deta, G22 = deta*deta;
319 
320 
321  const Real det = (G11*G22 - G12*G21);
322 
323  if (secure)
324  {
325  libmesh_assert_greater (det, 0.);
326  libmesh_assert_greater (std::abs(det), 1.e-10);
327  }
328 
329  const Real inv_det = 1./det;
330 
331  const Real
332  Ginv11 = G22*inv_det,
333  Ginv12 = -G12*inv_det,
334 
335  Ginv21 = -G21*inv_det,
336  Ginv22 = G11*inv_det;
337 
338 
339  const Real dxidelta = dxi*delta;
340  const Real detadelta = deta*delta;
341 
342  // compute only the first two coordinates.
343  dp(0) = (Ginv11*dxidelta + Ginv12*detadelta);
344  dp(1) = (Ginv21*dxidelta + Ginv22*detadelta);
345 
346  break;
347  }
348 
349 
350 
351  // Some other dimension?
352  default:
353  libmesh_error();
354  } // end switch(Dim), dp now computed
355 
356 
357 
358  // determine the error in computing the local coordinates
359  // in the base: ||P_n+1 - P_n||
360  inverse_map_error = dp.size();
361 
362 
363  // P_n+1 = P_n + dp
364  p.add (dp);
365 
366 
367  // Increment the iteration count.
368  cnt++;
369 
370 
371  // Watch for divergence of Newton's
372  // method.
373  if (cnt > 10)
374  {
375  if (secure || !secure)
376  {
377  libmesh_here();
378  {
379  libMesh::err << "WARNING: Newton scheme has not converged in "
380  << cnt << " iterations:\n"
381  << " physical_point="
382  << physical_point
383  << " dp="
384  << dp
385  << " p="
386  << p
387  << " error=" << inverse_map_error
388  << std::endl;
389  }
390  }
391 
392  if (cnt > 20)
393  {
394  libMesh::err << "ERROR: Newton scheme FAILED to converge in "
395  << cnt << " iterations!" << std::endl;
396  libmesh_error();
397  }
398 
399  // else
400  // {
401  // break;
402  // }
403  }
404  }
405  while (inverse_map_error > tolerance);
406 
407 
408 
409  // 4.
410  //
411  // Now that we have the local coordinates in the base,
412  // compute the interpolated radial distance a(s,t) \p a_interpolated
413  if (interpolated)
414  switch (Dim)
415  {
416  case 1:
417  {
418  Real a_interpolated = Point( inf_elem->point(0)
419  - inf_elem->point(n_base_mapping_sf) ).size();
420 
421  p(0) = 1. - 2*a_interpolated/physical_point(0);
422 
423 #ifdef DEBUG
424  // the radial distance should always be >= -1.
425 
426  if (p(0)+1 < tolerance)
427  {
428  libmesh_here();
429  libMesh::err << "WARNING: radial distance p(0) is "
430  << p(0)
431  << std::endl;
432  }
433 #endif
434 
435  break;
436  }
437 
438 
439  case 2:
440  {
441  Real a_interpolated = 0.;
442 
443  // the distance between the origin and the physical point
444  const Real fp_o_dist = Point(o-physical_point).size();
445 
446  for (unsigned int i=0; i<n_base_mapping_sf; i++)
447  {
448  // the radial distance of the i-th base mapping point
449  const Real dist_i = Point( inf_elem->point(i)
450  - inf_elem->point(i+n_base_mapping_sf) ).size();
451  // weight with the corresponding shape function
452  a_interpolated += dist_i * FE<1,LAGRANGE>::shape(base_mapping_elem_type,
453  base_mapping_order,
454  i,
455  p);
456  }
457 
458  p(1) = 1. - 2*a_interpolated/fp_o_dist;
459 
460 #ifdef DEBUG
461  // the radial distance should always be >= -1.
462 
463  // if (p(1)+1 < tolerance)
464  // {
465  // libmesh_here();
466  // libMesh::err << "WARNING: radial distance p(1) is "
467  // << p(1)
468  // << std::endl;
469  // }
470 #endif
471 
472  break;
473  }
474 
475 
476  case 3:
477  {
478  Real a_interpolated = 0.;
479 
480 
481  // the distance between the origin and the physical point
482  const Real fp_o_dist = Point(o-physical_point).size();
483 
484  for (unsigned int i=0; i<n_base_mapping_sf; i++)
485  {
486  // the radial distance of the i-th base mapping point
487  const Real dist_i = Point( inf_elem->point(i)
488  - inf_elem->point(i+n_base_mapping_sf) ).size();
489 
490  // weight with the corresponding shape function
491  a_interpolated += dist_i * FE<2,LAGRANGE>::shape(base_mapping_elem_type,
492  base_mapping_order,
493  i,
494  p);
495 
496  }
497 
498  p(2) = 1. - 2*a_interpolated/fp_o_dist;
499 
500 #ifdef DEBUG
501 
502 
503  // the radial distance should always be >= -1.
504 
505  // if (p(2)+1 < tolerance)
506  // {
507  // libmesh_here();
508  // libMesh::err << "WARNING: radial distance p(2) is "
509  // << p(2)
510  // << std::endl;
511  // }
512 #endif
513 
514  break;
515  }
516 
517  default:
518  libmesh_error();
519  } // end switch(Dim), p fully computed, including radial part
520 
521  // if we do not want the interpolated distance, then
522  // use newton iteration to get the actual distance
523  else
524  {
525  // distance from the physical point to the ifem origin
526  const Real fp_o_dist = Point(o-physical_point).size();
527 
528  // the distance from the intersection on the
529  // base to the origin
530  const Real a_dist = intersection.size();
531 
532  // element coordinate in radial direction
533  // here our first guess is 0.
534  Real v = 0.;
535 
536  // the order of the radial mapping
537  const Order radial_mapping_order (Radial::mapping_order());
538 
539  unsigned int cnt2 = 0;
540  inverse_map_error = 0.;
541 
542  // Newton iteration in 1-D
543  do
544  {
545  // the mapping in radial direction
546  // note that we only have two mapping functions in
547  // radial direction
548  const Real r = a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval (v, radial_mapping_order, 0)
549  + 2. * a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval (v, radial_mapping_order, 1);
550 
551  const Real dr = a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (v, radial_mapping_order, 0)
552  + 2. * a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (v, radial_mapping_order, 1);
553 
554  const Real G = dr*dr;
555  const Real Ginv = 1./G;
556 
557  const Real delta = fp_o_dist - r;
558  const Real drdelta = dr*delta;
559 
560  Real dp = Ginv*drdelta;
561 
562  // update the radial coordinate
563  v += dp;
564 
565  // note that v should be smaller than 1,
566  // since radial mapping function tends to infinity
567  if (v >= 1.)
568  v = .9999;
569 
570  inverse_map_error = std::fabs(dp);
571 
572  // increment iteration count
573  cnt2 ++;
574  if (cnt2 > 20)
575  {
576 
577  libMesh::err << "ERROR: 1D Newton scheme FAILED to converge"
578  << std::endl;
579 
580  libmesh_error();
581  }
582 
583 
584  }
585  while (inverse_map_error > tolerance);
586 
587  switch (Dim)
588  {
589  case 1:
590  {
591  p(0) = v;
592  break;
593  }
594  case 2:
595  {
596  p(1) = v;
597  break;
598  }
599  case 3:
600  {
601  p(2) = v;
602  break;
603  }
604  default:
605  libmesh_error();
606  }
607  }
608 
609  // If we are in debug mode do a sanity check. Make sure
610  // the point \p p on the reference element actually does
611  // map to the point \p physical_point within a tolerance.
612 #ifdef DEBUG
613  /*
614  const Point check = InfFE<Dim,T_radial,T_map>::map (inf_elem, p);
615  const Point diff = physical_point - check;
616 
617  if (diff.size() > tolerance)
618  {
619  libmesh_here();
620  libMesh::err << "WARNING: diff is "
621  << diff.size()
622  << std::endl;
623  }
624  */
625 #endif
626 
627 
628 
629 
630  // Stop logging the map inversion.
631  STOP_LOG("inverse_map()", "InfFE");
632 
633  return p;
634 }
635 
636 template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
638  const std::vector<Point>& physical_points,
639  std::vector<Point>& reference_points,
640  const Real tolerance,
641  const bool secure)
642 {
643  // The number of points to find the
644  // inverse map of
645  const std::size_t n_points = physical_points.size();
646 
647  // Resize the vector to hold the points
648  // on the reference element
649  reference_points.resize(n_points);
650 
651  // Find the coordinates on the reference
652  // element of each point in physical space
653  for (unsigned int p=0; p<n_points; p++)
654  reference_points[p] =
655  InfFE<Dim,T_radial,T_map>::inverse_map (elem, physical_points[p], tolerance, secure, false);
656 }
657 
658 
659 
660 
661 //--------------------------------------------------------------
662 // Explicit instantiations using the macro from inf_fe_macro.h
663 //INSTANTIATE_INF_FE(1,CARTESIAN);
664 
665 //INSTANTIATE_INF_FE(2,CARTESIAN);
666 
667 //INSTANTIATE_INF_FE(3,CARTESIAN);
668 
669 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,Point,inverse_map(const Elem*,const Point&,const Real,const bool,const bool));
670 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,Point,inverse_map(const Elem*,const Point&,const Real,const bool,const bool));
671 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,Point,inverse_map(const Elem*,const Point&,const Real,const bool,const bool));
672 
673 INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,void,inverse_map(const Elem*,const std::vector<Point>&,std::vector<Point>&,const Real, const bool));
674 INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,void,inverse_map(const Elem*,const std::vector<Point>&,std::vector<Point>&,const Real, const bool));
675 INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,void,inverse_map(const Elem*,const std::vector<Point>&,std::vector<Point>&,const Real, const bool));
676 
677 
678 } // namespace libMesh
679 
680 
681 #endif //ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS

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

Hosted By:
SourceForge.net Logo