libMesh::InfFEMap Class Reference

Class that encapsulates mapping (i.e. More...

#include <inf_fe_map.h>

Static Public Member Functions
static Point	map (const unsigned int dim, const Elem *inf_elem, const Point &reference_point)
static Point	inverse_map (const unsigned int dim, const Elem *elem, const Point &p, const Real tolerance=TOLERANCE, const bool secure=true)
static void	inverse_map (const unsigned int dim, const Elem *elem, const std::vector< Point > &physical_points, std::vector< Point > &reference_points, const Real tolerance=TOLERANCE, const bool secure=true)
	Takes a number points in physical space (in the physical_points vector) and finds their location on the reference element for the input element elem. More...
static Real	eval (Real v, Order o, unsigned int i)
static Real	eval_deriv (Real v, Order o, unsigned int i)

Detailed Description

Class that encapsulates mapping (i.e.

from physical space to reference space and vice-versa) operations on infinite elements.

Definition at line 47 of file inf_fe_map.h.

Member Function Documentation

eval()

Real libMesh::InfFEMap::eval ( Real  v, Order  o, unsigned int  i  )

static

Returns

The value of the \( i^{th} \) polynomial evaluated at v. This method provides the approximation in radial direction for the overall shape functions, which is defined in InfFE::shape(). This method is allowed to be static, since it is independent of dimension and base_family. It is templated, though, w.r.t. to radial FEFamily.

The value of the \( i^{th} \) mapping shape function in radial direction evaluated at v when T_radial == INFINITE_MAP. Currently, only one specific mapping shape is used. Namely the one by Marques JMMC, Owen DRJ: Infinite elements in quasi-static materially nonlinear problems, Computers and Structures, 1984.

Definition at line 27 of file inf_fe_map_eval.C.

References libMesh::libmesh_assert().

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::eval(), and map().

{
  libmesh_assert (-1.-1.e-5 <= v && v < 1.);

  switch (i)
    {
    case 0:
      return -2.*v/(1.-v);
    case 1:
      return (1.+v)/(1.-v);

    default:
      libmesh_error_msg("bad index i = " << i);
    }
}

eval_deriv()

Real libMesh::InfFEMap::eval_deriv ( Real  v, Order  o, unsigned int  i  )

static

Returns

The value of the first derivative of the \( i^{th} \) polynomial at coordinate v. See eval for details.

Definition at line 45 of file inf_fe_map_eval.C.

References libMesh::libmesh_assert().

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::eval_deriv().

{
  libmesh_assert (-1.-1.e-5 <= v && v < 1.);

  switch (i)
    {
    case 0:
      return -2./((1.-v)*(1.-v));
    case 1:
      return 2./((1.-v)*(1.-v));

    default:
      libmesh_error_msg("bad index i = " << i);
    }
}

inverse_map() [1/2]

Point libMesh::InfFEMap::inverse_map ( const unsigned int  dim, const Elem *  elem, const Point &  p, const Real  tolerance = TOLERANCE, const bool  secure = true  )

static

Returns

The location (on the reference element) of the point p located in physical space. First, the location in the base face is computed. This requires inverting the (possibly nonlinear) transformation map in the base, so it is not trivial. The optional parameter tolerance defines how close is "good enough." The map inversion iteration computes the sequence \( \{ p_n \} \), and the iteration is terminated when \( \|p - p_n\| < \mbox{\texttt{tolerance}} \). Once the base face point is determined, the radial local coordinate is directly evaluated. If interpolated is true, the interpolated distance from the base element to the infinite element origin is used for the map in radial direction.

Definition at line 96 of file inf_fe_map.C.

References libMesh::TypeVector< T >::add_scaled(), libMesh::InfFEBase::build_elem(), dim, libMesh::FEMap::inverse_map(), libMesh::libmesh_assert(), libMesh::libmesh_isinf(), map(), libMesh::FE< Dim, T >::n_shape_functions(), libMesh::TypeVector< T >::norm(), libMesh::Elem::origin(), libMesh::Real, libMesh::FE< Dim, T >::shape(), libMesh::FE< Dim, T >::shape_deriv(), libMesh::TypeTensor< T >::solve(), and libMesh::TOLERANCE.

Referenced by libMesh::InfQuad4::contains_point(), libMesh::InfPrism::contains_point(), libMesh::InfHex::contains_point(), inverse_map(), libMesh::FEMap::inverse_map(), and libMesh::InfFE< Dim, T_radial, T_map >::inverse_map().

{
  libmesh_assert(inf_elem);
  libmesh_assert_greater_equal (tolerance, 0.);
  libmesh_assert(dim > 0);

  // Start logging the map inversion.
  LOG_SCOPE("inverse_map()", "InfFEMap");

  // The strategy is:
  // compute the intersection of the line
  // physical_point - origin with the base element,
  // find its internal coordinatels using FEMap::inverse_map():
  // The radial part can then be computed directly later on.

  // 1.)
  // build a base element to do the map inversion in the base face
  std::unique_ptr<const Elem> base_elem = InfFEBase::build_elem (inf_elem);

  // The point on the reference element (which we are looking for).
  // start with an invalid guess:
  Point p;
  p(dim-1)=-2.;

  // 2.)
  // Now find the intersection of a plane represented by the base
  // element nodes and the line given by the origin of the infinite
  // element and the physical point.
  Point intersection;

  // the origin of the infinite element
  const Point o = inf_elem->origin();

  switch (dim)
    {
      // unnecessary for 1D
    case 1:
      break;

    case 2:
      libmesh_error_msg("ERROR: InfFE::inverse_map is not yet implemented in 2d");

    case 3:
      {
        const Point xi ( base_elem->point(1) - base_elem->point(0));
        const Point eta( base_elem->point(2) - base_elem->point(0));
        const Point zeta( physical_point - o);

        // normal vector of the base elements plane
        Point n(xi.cross(eta));
        Real c_factor = (base_elem->point(0) - o)*n/(zeta*n) - 1.;

        // Check whether the above system is ill-posed.  It should
        // only happen when \p physical_point is not in \p
        // inf_elem. In this case, any point that is not in the
        // element is a valid answer.
        if (libmesh_isinf(c_factor))
          return p;

        // Compute the intersection with
        // {intersection} = {physical_point} + c*({physical_point}-{o}).
        intersection.add_scaled(physical_point,1.+c_factor);
        intersection.add_scaled(o,-c_factor);

        // For non-planar elements, the intersecting point is obtained via Newton-iteration
        if (!base_elem->has_affine_map())
          {
            unsigned int iter_max = 20;

            // the number of shape functions needed for the base_elem
            unsigned int n_sf = FE<2,LAGRANGE>::n_shape_functions(base_elem->type(),base_elem->default_order());

            // guess base element coordinates: p=xi,eta,0
            // in first iteration, never run with 'secure' to avoid false warnings.
            Point ref_point= FEMap::inverse_map(dim-1, base_elem.get(), intersection,
                                                tolerance, false);

            // Newton iteration
            for(unsigned int it=0; it<iter_max; ++it)
              {
                // Get the shape function and derivative values at the reference coordinate
                // phi.size() == dphi.size()
                Point dxyz_dxi;
                Point dxyz_deta;

                Point intersection_guess;
                for(unsigned int i=0; i<n_sf; ++i)
                  {

                    intersection_guess += base_elem->node_ref(i) * FE<2,LAGRANGE>::shape(base_elem->type(),
                                                                                         base_elem->default_order(),
                                                                                         i,
                                                                                         ref_point);

                    dxyz_dxi += base_elem->node_ref(i) * FE<2,LAGRANGE>::shape_deriv(base_elem->type(),
                                                                                     base_elem->default_order(),
                                                                                     i,
                                                                                     0, // d()/dxi
                                                                                     ref_point);

                    dxyz_deta+= base_elem->node_ref(i) * FE<2,LAGRANGE>::shape_deriv(base_elem->type(),
                                                                                     base_elem->default_order(),
                                                                                     i,
                                                                                     1, // d()/deta
                                                                                     ref_point);
                  } // for i

                TypeVector<Real> F(physical_point + c_factor*(physical_point-o) - intersection_guess);

                TypeTensor<Real> J;
                J(0,0) = (physical_point-o)(0);
                J(0,1) = -dxyz_dxi(0);
                J(0,2) = -dxyz_deta(0);
                J(1,0) = (physical_point-o)(1);
                J(1,1) = -dxyz_dxi(1);
                J(1,2) = -dxyz_deta(1);
                J(2,0) = (physical_point-o)(2);
                J(2,1) = -dxyz_dxi(2);
                J(2,2) = -dxyz_deta(2);

                // delta will be the newton step
                TypeVector<Real> delta;
                J.solve(F,delta);

                // check for convergence
                Real tol = std::min( TOLERANCE*0.01, TOLERANCE*base_elem->hmax() );
                if ( delta.norm() < tol )
                  {
                    // newton solver converged, now make sure it converged to a point on the base_elem
                    if (base_elem->contains_point(intersection_guess,TOLERANCE*0.1))
                      {
                        intersection = intersection_guess;
                      }
                    //else
                    //  {
                    //     err<<" Error: inverse_map failed!"<<std::endl;
                    //  }
                    break; // break out of 'for it'
                  }
                else
                  {
                    c_factor     -= delta(0);
                    ref_point(0) -= delta(1);
                    ref_point(1) -= delta(2);
                  }

              }

          }
        break;
      }

    default:
      libmesh_error_msg("Invalid dim = " << dim);
    }

  // 3.)
  // Now we have the intersection-point (projection of physical point onto base-element).
  // Lets compute its internal coordinates (being p(0) and p(1)):
  p= FEMap::inverse_map(dim-1, base_elem.get(), intersection,
                        tolerance, secure, secure);

  // 4.
  // Now that we have the local coordinates in the base,
  // we compute the radial distance with Newton iteration.

  // distance from the physical point to the ifem origin
  const Real fp_o_dist = Point(o-physical_point).norm();

  // the distance from the intersection on the
  // base to the origin
  const Real a_dist = Point(o-intersection).norm();

  // element coordinate in radial direction

  // fp_o_dist is at infinity.
  if (libmesh_isinf(fp_o_dist))
    {
      p(dim-1)=1;
      return p;
    }

  // when we are somewhere in this element:
  Real v = 0;

  // For now we're sticking with T_map == CARTESIAN
  // if (T_map == CARTESIAN)
  v = 1.-2.*a_dist/fp_o_dist;
  // else
  //   libmesh_not_implemented();


  p(dim-1)=v;
#ifdef DEBUG
  // first check whether we are in the reference-element:
  if ((-1.-1.e-5 < v && v < 1.) || secure)
    {
      const Point check = map (dim, inf_elem, p);
      const Point diff  = physical_point - check;

      if (diff.norm() > tolerance)
        libmesh_warning("WARNING:  diff is " << diff.norm());
    }
#endif

  return p;
}

inverse_map() [2/2]

void libMesh::InfFEMap::inverse_map ( const unsigned int  dim, const Elem *  elem, const std::vector< Point > &  physical_points, std::vector< Point > &  reference_points, const Real  tolerance = TOLERANCE, const bool  secure = true  )

static

Takes a number points in physical space (in the physical_points vector) and finds their location on the reference element for the input element elem.

The values on the reference element are returned in the vector reference_points

Definition at line 310 of file inf_fe_map.C.

References dim, and inverse_map().

{
  // The number of points to find the
  // inverse map of
  const std::size_t n_points = physical_points.size();

  // Resize the vector to hold the points
  // on the reference element
  reference_points.resize(n_points);

  // Find the coordinates on the reference
  // element of each point in physical space
  for (unsigned int p=0; p<n_points; p++)
    reference_points[p] =
      inverse_map (dim, elem, physical_points[p], tolerance, secure);
}

map()

Point libMesh::InfFEMap::map ( const unsigned int  dim, const Elem *  inf_elem, const Point &  reference_point  )

static

Returns

The location (in physical space) of the point p located on the reference element.

Definition at line 40 of file inf_fe_map.C.

References libMesh::TypeVector< T >::add_scaled(), libMesh::InfFEBase::build_elem(), dim, eval(), libMesh::libmesh_assert(), libMesh::FEMap::map(), libMesh::InfFERadial::mapping_order(), libMesh::Elem::origin(), libMesh::Elem::point(), and libMesh::Real.

Referenced by libMesh::InfFE< Dim, T_radial, T_map >::compute_shape_functions(), inverse_map(), libMesh::FEMap::map(), and libMesh::InfFE< Dim, T_radial, T_map >::map().

{
  libmesh_assert(inf_elem);
  libmesh_assert_not_equal_to (dim, 0);

  std::unique_ptr<const Elem> base_elem = InfFEBase::build_elem (inf_elem);

  const Real         v                    (reference_point(dim-1));

  // map in the base face
  Point base_point;
  switch (dim)
    {
    case 1:
      base_point = inf_elem->point(0);
      break;
    case 2:
    case 3:
      base_point = FEMap::map (dim-1, base_elem.get(), reference_point);
      break;
    default:
#ifdef DEBUG
      libmesh_error_msg("Unknown dim = " << dim);
#endif
      break;
    }

  // This is the same as the algorithm used below,
  // but is more explicit in the actual form in the end.
  //
  // NOTE: the form used below can be implemented to yield
  // more general/flexible mappings, but the current form is
  // used e.g. for \p inverse_map() and \p reinit() explicitly.
  return (base_point-inf_elem->origin())*2/(1-v)+inf_elem->origin();

#if 0 // Leave this documented, but w/o unreachable-code warnings
  const Order radial_mapping_order (InfFERadial::mapping_order());

  // map in the outer node face not necessary. Simply
  // compute the outer_point = base_point + (base_point-origin)
  const Point outer_point (base_point * 2. - inf_elem->origin());

  Point p;

  // there are only two mapping shapes in radial direction
  p.add_scaled (base_point,  eval (v, radial_mapping_order, 0));
  p.add_scaled (outer_point, eval (v, radial_mapping_order, 1));

  return p;
#endif
}

---

The documentation for this class was generated from the following files:

• include/fe/inf_fe_map.h
• src/fe/inf_fe_map.C
• src/fe/inf_fe_map_eval.C