inf_fe_map.c

一个用来实现偏微分方程中网格的计算库

	    
	    // compute only the first two coordinates.
	    dp(0) = (Ginv11*dxidelta + Ginv12*detadelta);
	    dp(1) = (Ginv21*dxidelta + Ginv22*detadelta);

	    break;
	  }


	  
	  // Some other dimension?
	default:
	  libmesh_error();
	} // end switch(Dim), dp now computed


      
      // determine the error in computing the local coordinates
      // in the base: ||P_n+1 - P_n||
      inverse_map_error = dp.size();

      
      // P_n+1 = P_n + dp
      p.add (dp);

      
      // Increment the iteration count.
      cnt++;

      
      // Watch for divergence of Newton's
      // method.
      if (cnt > 10)
	{
	  if (secure || !secure)
	    {
	      here();
	      {
		std::cerr << "WARNING: Newton scheme has not converged in "
			  << cnt << " iterations:\n" 
			  << "   physical_point="
			  << physical_point
			  << "   dp="
			  << dp
			  << "   p="
			  << p
			  << "   error=" << inverse_map_error
			  << std::endl;
	      }
	    }

	  if (cnt > 20)
	    {
	      std::cerr << "ERROR: Newton scheme FAILED to converge in "
			<< cnt << " iterations!" << std::endl;
	      libmesh_error();
	    }

	  // else
	  //  {
	  //    break;
	  //  }
	}
    }
  while (inverse_map_error > tolerance);


  
  // 4.
  //
  // Now that we have the local coordinates in the base,
  // compute the interpolated radial distance a(s,t) \p a_interpolated
  if (interpolated)
    switch (Dim)
      {
      case 1:
	{
	  Real a_interpolated = Point( inf_elem->point(0) 
				       - inf_elem->point(n_base_mapping_sf) ).size();
	  
	  p(0) = 1. - 2*a_interpolated/physical_point(0);
	  
#ifdef DEBUG
	  // the radial distance should always be >= -1.
	  
	  if (p(0)+1 < tolerance)
	    {
	      here();
	      std::cerr << "WARNING: radial distance p(0) is "
			<< p(0)
			<< std::endl;
	    }
#endif
	  
	  break;
	}
	
	
      case 2:
	{
	  Real a_interpolated = 0.;
	  
	  // the distance between the origin and the physical point
	  const Real fp_o_dist = Point(o-physical_point).size(); 
	  
	  for (unsigned int i=0; i<n_base_mapping_sf; i++)
	    {
	      // the radial distance of the i-th base mapping point
	      const Real dist_i = Point( inf_elem->point(i) 
					 - inf_elem->point(i+n_base_mapping_sf) ).size();
	      // weight with the corresponding shape function
	      a_interpolated += dist_i * FE<1,LAGRANGE>::shape(base_mapping_elem_type,
							       base_mapping_order,
							       i,
							       p);
	    }
	  
	  p(1) = 1. - 2*a_interpolated/fp_o_dist;
	  
#ifdef DEBUG
	  // the radial distance should always be >= -1.
	  
	  // if (p(1)+1 < tolerance)
	  //  {
	  //    here();
	  //    std::cerr << "WARNING: radial distance p(1) is "
	  //	      << p(1)
	  //	      << std::endl;
	  //  }
#endif
	  
	  break;
	}
	
	
      case 3:
	{
	  Real a_interpolated = 0.;
	  
	  
	  // the distance between the origin and the physical point
	  const Real fp_o_dist = Point(o-physical_point).size(); 
	  
	  for (unsigned int i=0; i<n_base_mapping_sf; i++)
	    {
	      // the radial distance of the i-th base mapping point
	      const Real dist_i = Point( inf_elem->point(i) 
					 - inf_elem->point(i+n_base_mapping_sf) ).size();
	      
	      // weight with the corresponding shape function
	      a_interpolated += dist_i * FE<2,LAGRANGE>::shape(base_mapping_elem_type,
							       base_mapping_order,
							       i,
							       p);
	      
	    }
	  
	  p(2) = 1. - 2*a_interpolated/fp_o_dist;
	  
#ifdef DEBUG
	  
	  
	  // the radial distance should always be >= -1.
	  
	  // if (p(2)+1 < tolerance)
	  //  {
	  // here();
	  // std::cerr << "WARNING: radial distance p(2) is "
	  //	      << p(2)
	  //	      << std::endl;
	  //  }
#endif
	  
	  break;
	}
	
      default:
	libmesh_error();
      } // end switch(Dim), p fully computed, including radial part

  // if we do not want the interpolated distance, then
  // use newton iteration to get the actual distance
  else
    {
      // distance from the physical point to the ifem origin
      const Real fp_o_dist = Point(o-physical_point).size(); 

      // the distance from the intersection on the
      // base to the origin
      const Real a_dist = intersection.size();
      
      // element coordinate in radial direction
      // here our first guess is 0.
      Real v = 0.;
      
      // the order of the radial mapping
      const Order radial_mapping_order (Radial::mapping_order());    
      
      unsigned int cnt = 0;
      inverse_map_error = 0.;
      
      // Newton iteration in 1-D
      do
	{
	  // the mapping in radial direction
	  // note that we only have two mapping functions in
	  // radial direction
	  const Real r = a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval (v, radial_mapping_order, 0)
	    + 2. * a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval (v, radial_mapping_order, 1);
	  
	  const Real dr = a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (v, radial_mapping_order, 0)
	    + 2. * a_dist * InfFE<Dim,INFINITE_MAP,T_map>::eval_deriv (v, radial_mapping_order, 1);
	  
	  const Real G = dr*dr;
	  const Real Ginv = 1./G;
	  
	  const Real delta = fp_o_dist - r;
	  const Real drdelta = dr*delta;
	  
	  Real dp = Ginv*drdelta;
	  
	  // update the radial coordinate
	  v += dp;
	  
	  // note that v should be smaller than 1,
	  // since radial mapping function tends to infinity
	  if (v >= 1.)
	    v = .9999;
	  
	  inverse_map_error = std::fabs(dp);
	  
	  // increment iteration count
	  cnt ++;
	  if (cnt > 20)
	    {
	      
	      std::cerr << "ERROR: 1D Newton scheme FAILED to converge"
			<< std::endl;
	      
	      libmesh_error();
	    }
	  
	  
	}
      while (inverse_map_error > tolerance);

      switch (Dim)
	{
	case 1:
	  {
	    p(0) = v;
	    break;
	  }
	case 2:
	  {
	    p(1) = v;
	    break;
	  }
	case 3:
	  {
	    p(2) = v;
	    break;
	  }
	default:
	  libmesh_error();
	}
    }
  
  // If we are in debug mode do a sanity check.  Make sure
  // the point \p p on the reference element actually does
  // map to the point \p physical_point within a tolerance. 
#ifdef DEBUG
  /*	
	const Point check = InfFE<Dim,T_radial,T_map>::map (inf_elem, p);
	const Point diff  = physical_point - check;
  
	if (diff.size() > tolerance)
	{
	here();
	std::cerr << "WARNING:  diff is "
	<< diff.size()
	<< std::endl;
	}
  */
#endif



  
  // Stop logging the map inversion.
  STOP_LOG("inverse_map()", "InfFE");
  
  return p;
}

template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::inverse_map (const Elem* elem,
					     const std::vector<Point>& physical_points,
					     std::vector<Point>&       reference_points,
					     const Real tolerance,
					     const bool secure)
{
  // The number of points to find the
  // inverse map of
  const unsigned int 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] =
      InfFE<Dim,T_radial,T_map>::inverse_map (elem, physical_points[p], tolerance, secure, false);
}




//--------------------------------------------------------------
// Explicit instantiations using the macro from inf_fe_macro.h
//INSTANTIATE_INF_FE(1,CARTESIAN);

//INSTANTIATE_INF_FE(2,CARTESIAN);

//INSTANTIATE_INF_FE(3,CARTESIAN);

INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,Point,inverse_map(const Elem*,const Point&,const Real,const bool,const bool));
INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,Point,inverse_map(const Elem*,const Point&,const Real,const bool,const bool));
INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,Point,inverse_map(const Elem*,const Point&,const Real,const bool,const bool));

INSTANTIATE_INF_FE_MBRF(1,CARTESIAN,void,inverse_map(const Elem*,const std::vector<Point>&,std::vector<Point>&,const Real, const bool));
INSTANTIATE_INF_FE_MBRF(2,CARTESIAN,void,inverse_map(const Elem*,const std::vector<Point>&,std::vector<Point>&,const Real, const bool));
INSTANTIATE_INF_FE_MBRF(3,CARTESIAN,void,inverse_map(const Elem*,const std::vector<Point>&,std::vector<Point>&,const Real, const bool));



#endif //ifdef ENABLE_INFINITE_ELEMENTS