fe_map.c

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

  // Resize the vectors to hold data at the quadrature points
  this->resize_map_vectors(n_qp);

  // Compute map at all quadrature points
  for (unsigned int p=0; p!=n_qp; p++)
    this->compute_single_point_map(qw, elem, p);
  
  // Stop logging the map computation.
  STOP_LOG("compute_map()", "FE");  
}




template <unsigned int Dim, FEFamily T>
Point FE<Dim,T>::map (const Elem* elem,
		      const Point& reference_point)
{
  libmesh_assert (elem != NULL);
    
  Point p;

  const ElemType type     = elem->type();
  const Order order       = elem->default_order();
  const unsigned int n_sf = FE<Dim,LAGRANGE>::n_shape_functions(type, order);

  // Lagrange basis functions are used for mapping
  for (unsigned int i=0; i<n_sf; i++)
    p.add_scaled (elem->point(i),
		  FE<Dim,LAGRANGE>::shape(type,
					  order,
					  i,
					  reference_point)
		  );

  return p;
}



template <unsigned int Dim, FEFamily T>
Point FE<Dim,T>::map_xi (const Elem* elem,
			 const Point& reference_point)
{
  libmesh_assert (elem != NULL);
    
  Point p;

  const ElemType type     = elem->type();
  const Order order       = elem->default_order();
  const unsigned int n_sf = FE<Dim,LAGRANGE>::n_shape_functions(type, order);

  // Lagrange basis functions are used for mapping    
  for (unsigned int i=0; i<n_sf; i++)
    p.add_scaled (elem->point(i),
		  FE<Dim,LAGRANGE>::shape_deriv(type,
						order,
						i,
						0,
						reference_point)
		  );
    
  return p;
}



template <unsigned int Dim, FEFamily T>
Point FE<Dim,T>::map_eta (const Elem* elem,
			  const Point& reference_point)
{
  libmesh_assert (elem != NULL);
    
  Point p;

  const ElemType type     = elem->type();
  const Order order       = elem->default_order();
  const unsigned int n_sf = FE<Dim,LAGRANGE>::n_shape_functions(type, order);
  
  // Lagrange basis functions are used for mapping
  for (unsigned int i=0; i<n_sf; i++)
    p.add_scaled (elem->point(i),
		  FE<Dim,LAGRANGE>::shape_deriv(type,
						order,
						i,
						1,
						reference_point)
		  );
    
  return p;
}



template <unsigned int Dim, FEFamily T>
Point FE<Dim,T>::map_zeta (const Elem* elem,
			   const Point& reference_point)
{
  libmesh_assert (elem != NULL);
    
  Point p;

  const ElemType type     = elem->type();
  const Order order       = elem->default_order();
  const unsigned int n_sf = FE<Dim,LAGRANGE>::n_shape_functions(type, order);

  // Lagrange basis functions are used for mapping
  for (unsigned int i=0; i<n_sf; i++)
    p.add_scaled (elem->point(i),
		  FE<Dim,LAGRANGE>::shape_deriv(type,
						order,
						i,
						2,
						reference_point)
		  );
    
  return p;
}



template <unsigned int Dim, FEFamily T>
Point FE<Dim,T>::inverse_map (const Elem* elem,
			      const Point& physical_point,
			      const Real tolerance,
			      const bool secure)
{
  libmesh_assert (elem != NULL);
  libmesh_assert (tolerance >= 0.);

  
  // Start logging the map inversion.  
  START_LOG("inverse_map()", "FE");
  
  // How much did the point on the reference
  // element change by in this Newton step?
  Real inverse_map_error = 0.;
  
  //  The point on the reference element.  This is
  //  the "initial guess" for Newton's method.  The
  //  centroid seems like a good idea, but computing
  //  it is a little more intensive than, say taking
  //  the zero point.  
  // 
  //  Convergence should be insensitive of this choice
  //  for "good" elements.
  Point p; // the zero point.  No computation required

  //  The number of iterations in the map inversion process.
  unsigned int cnt = 0;




  //  Newton iteration loop.
  do
    {
      //  Where our current iterate \p p maps to.
      const Point physical_guess = FE<Dim,T>::map (elem, p);

      //  How far our current iterate is from the actual point.
      const Point delta = physical_point - physical_guess;

      //  Increment in current iterate \p p, will be computed.
      Point dp;


      //  The form of the map and how we invert it depends
      //  on the dimension that we are in.
      switch (Dim)
	{
	  
	  // ------------------------------------------------------------------
	  //  1D map inversion
	  // 
	  //  Here we find the point on a 1D reference element that maps to
	  //  the point \p physical_point in the domain.  This is a bit tricky
	  //  since we do not want to assume that the point \p physical_point
	  //  is also in a 1D domain.  In particular, this method might get
	  //  called on the edge of a 3D element, in which case
	  //  \p physical_point actually lives in 3D.
	case 1:
	  {
	    const Point dxi = FE<Dim,T>::map_xi (elem, p);
	    
	    //  Newton's method in this case looks like
	    // 
	    //  {X} - {X_n} = [J]*dp
	    // 
	    //  Where {X}, {X_n} are 3x1 vectors, [J] is a 3x1 matrix
	    //  d(x,y,z)/dxi, and we seek dp, a scalar.  Since the above
	    //  system is either overdetermined or rank-deficient, we will
	    //  solve the normal equations for this system
	    // 
	    //  [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
	    // 
	    //  which involves the trivial inversion of the scalar
	    //  G = [J]^T [J]
	    const Real G = dxi*dxi;
	    
	    if (secure)
	      libmesh_assert (G > 0.);
	    
	    const Real Ginv = 1./G;
	    
	    const Real  dxidelta = dxi*delta;
	    
	    dp(0) = Ginv*dxidelta;

            // Assume that no master elements have radius > 4
	    if (secure)
	      libmesh_assert (dp.size() < 4);

	    break;
	  }



	  // ------------------------------------------------------------------
	  //  2D map inversion
	  // 
	  //  Here we find the point on a 2D reference element that maps to
	  //  the point \p physical_point in the domain.  This is a bit tricky
	  //  since we do not want to assume that the point \p physical_point
	  //  is also in a 2D domain.  In particular, this method might get
	  //  called on the face of a 3D element, in which case
	  //  \p physical_point actually lives in 3D.
	case 2:
	  {
	    const Point dxi  = FE<Dim,T>::map_xi  (elem, p);
	    const Point deta = FE<Dim,T>::map_eta (elem, p);
	    
	    //  Newton's method in this case looks like
	    // 
	    //  {X} - {X_n} = [J]*{dp}
	    // 
	    //  Where {X}, {X_n} are 3x1 vectors, [J] is a 3x2 matrix
	    //  d(x,y,z)/d(xi,eta), and we seek {dp}, a 2x1 vector.  Since
	    //  the above system is either overdermined or rank-deficient,
	    //  we will solve the normal equations for this system
	    // 
	    //  [J]^T ({X} - {X_n}) = [J]^T [J] {dp}
	    // 
	    //  which involves the inversion of the 2x2 matrix
	    //  [G] = [J]^T [J]
	    const Real
	      G11 = dxi*dxi,  G12 = dxi*deta,
	      G21 = dxi*deta, G22 = deta*deta;
	    
	    
	    const Real det = (G11*G22 - G12*G21);
	    
	    if (secure)
	      libmesh_assert (det != 0.);

	    const Real inv_det = 1./det;
	    
	    const Real
	      Ginv11 =  G22*inv_det,
	      Ginv12 = -G12*inv_det,
	      
	      Ginv21 = -G21*inv_det,
	      Ginv22 =  G11*inv_det;
	    
	    
	    const Real  dxidelta  = dxi*delta;
	    const Real  detadelta = deta*delta;
	    
	    dp(0) = (Ginv11*dxidelta + Ginv12*detadelta);
	    dp(1) = (Ginv21*dxidelta + Ginv22*detadelta);

            // Assume that no master elements have radius > 4
	    if (secure)
	      libmesh_assert (dp.size() < 4);

	    break;
	  }
	  

	  
	  // ------------------------------------------------------------------
	  //  3D map inversion
	  // 
	  //  Here we find the point in a 3D reference element that maps to
	  //  the point \p physical_point in a 3D domain. Nothing special
	  //  has to happen here, since (unless the map is singular because
	  //  you have a BAD element) the map will be invertable and we can
	  //  apply Newton's method directly.
	case 3:
	  {
       	    const Point dxi   = FE<Dim,T>::map_xi   (elem, p);
	    const Point deta  = FE<Dim,T>::map_eta  (elem, p);
	    const Point dzeta = FE<Dim,T>::map_zeta (elem, p);
	    
	    //  Newton's method in this case looks like
	    // 
	    //  {X} = {X_n} + [J]*{dp}
	    // 
	    //  Where {X}, {X_n} are 3x1 vectors, [J] is a 3x3 matrix
	    //  d(x,y,z)/d(xi,eta,zeta), and we seek {dp}, a 3x1 vector.
	    //  Since the above system is nonsingular for invertable maps
	    //  we will solve 
	    // 
	    //  {dp} = [J]^-1 ({X} - {X_n})
	    // 
	    //  which involves the inversion of the 3x3 matrix [J]
	    const Real
	      J11 = dxi(0), J12 = deta(0), J13 = dzeta(0),
	      J21 = dxi(1), J22 = deta(1), J23 = dzeta(1),
	      J31 = dxi(2), J32 = deta(2), J33 = dzeta(2);
	    
	    const Real det = (J11*(J22*J33 - J23*J32) +
			      J12*(J23*J31 - J21*J33) +
			      J13*(J21*J32 - J22*J31));
	    
	    if (secure)
	      libmesh_assert (det != 0.);

	    const Real inv_det = 1./det;
	    
	    const Real
	      Jinv11 =  (J22*J33 - J23*J32)*inv_det,
	      Jinv12 = -(J12*J33 - J13*J32)*inv_det,
	      Jinv13 =  (J12*J23 - J13*J22)*inv_det,
	      
	      Jinv21 = -(J21*J33 - J23*J31)*inv_det,
	      Jinv22 =  (J11*J33 - J13*J31)*inv_det,
	      Jinv23 = -(J11*J23 - J13*J21)*inv_det,
	      
	      Jinv31 =  (J21*J32 - J22*J31)*inv_det,
	      Jinv32 = -(J11*J32 - J12*J31)*inv_det,
	      Jinv33 =  (J11*J22 - J12*J21)*inv_det;
	    
	    dp(0) = (Jinv11*delta(0) +
		     Jinv12*delta(1) +
		     Jinv13*delta(2));
	    
	    dp(1) = (Jinv21*delta(0) +
		     Jinv22*delta(1) +
		     Jinv23*delta(2));
	    
	    dp(2) = (Jinv31*delta(0) +
		     Jinv32*delta(1) +
		     Jinv33*delta(2));

            // Assume that no master elements have radius > 4
	    if (secure)
	      libmesh_assert (dp.size() < 4);

	    break;
	  }


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



      //  ||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.  Here's how it goes:
      //  (1) For good elements, we expect convergence in 10 iterations.
      //      - If called with (secure == true) and we have not yet converged
      //        print out a warning message.
      //      - If called with (secure == true) and we have not converged in
      //        20 iterations abort
      //  (2) This method may be called in cases when the target point is not
      //      inside the element and we have no business expecting convergence.
      //      For these cases if we have not converged in 10 iterations forget
      //      about it.
      if (cnt > 10)
	{
	  //  Warn about divergence when secure is true - this
	  //  shouldn't happen
	  if (secure)
	    {
	      here();
	      std::cerr << "WARNING: Newton scheme has not converged in "
			<< cnt << " iterations:" << std::endl
			<< "   physical_point="
			<< physical_point
			<< "   physical_guess="
			<< physical_guess
			<< "   dp="
			<< dp
			<< "   p="
			<< p
			<< "   error=" << inverse_map_error
                        << "   in element " << elem->id()
			<< std::endl;

	      if (cnt > 20)
		{
		  std::cerr << "ERROR: Newton scheme FAILED to converge in "
			    << cnt
			    << " iterations!"
                            << " in element " << elem->id()
			    << std::endl;

		  libmesh_error();
		}
	    }
	  //  Return a far off point when secure is false - this
	  //  should only happen when we're trying to map a point
	  //  that's outside the element
	  else
	    {
	      for (unsigned int i=0; i != Dim; ++i)
		p(i) = 1e6;
	      
	      STOP_LOG("inverse_map()", "FE");
	      return p;
	    }
	}
    }
  while (inverse_map_error > tolerance);



  //  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
  
  if (secure)
    {  
      const Point check = FE<Dim,T>::map (elem, p);
      const Point diff  = physical_point - check;
      
      if (diff.size() > tolerance)
	{
	  here();
	  std::cerr << "WARNING:  diff is "
		    << diff.size()
		    << std::endl
		    << " point="
		    << physical_point;
	  std::cerr << " local=" << check;
	  std::cerr << " lref= " << p;
	}
    }
  
#endif


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



template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::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] =
      FE<Dim,T>::inverse_map (elem, physical_points[p], tolerance, secure);
}



//--------------------------------------------------------------
// Explicit instantiations using the macro from fe_macro.h
INSTANTIATE_IMAP(1);
INSTANTIATE_IMAP(2);
INSTANTIATE_IMAP(3);

#if defined(__INTEL_COMPILER)
INSTANTIATE_MAP(1);
INSTANTIATE_MAP(2);
INSTANTIATE_MAP(3);
#endif