inf_fe.c

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

  // when computing the phase, we need the base approximations
  // therefore, initialize the phase here, but evaluate it 
  // in combine_base_radial().
  //
  // the weight, though, is only needed at the radial quadrature points, n_radial_qp.
  // but for a uniform interface to the protected data fields
  // the weight data field (which are accessible from the outside) are expanded to n_total_qp.
  weight.resize      (n_total_qp);
  dweightdv.resize   (n_total_qp);
  dweight.resize     (n_total_qp);

  dphase.resize      (n_total_qp);
  dphasedxi.resize   (n_total_qp);
  dphasedeta.resize  (n_total_qp);
  dphasedzeta.resize (n_total_qp);

  // this vector contains the integration weights for the combined quadrature rule
  _total_qrule_weights.resize(n_total_qp);


  // -----------------------------------------------------------------
  // InfFE's data fields phi, dphi, dphidx, phi_map etc hold the _total_
  // shape and mapping functions, respectively
  {
    phi.resize     (n_total_approx_shape_functions);
    dphi.resize    (n_total_approx_shape_functions);
    dphidx.resize  (n_total_approx_shape_functions);
    dphidy.resize  (n_total_approx_shape_functions);
    dphidz.resize  (n_total_approx_shape_functions);
    dphidxi.resize (n_total_approx_shape_functions);
#ifdef ENABLE_SECOND_DERIVATIVES
    static bool warning_given = false;
    if (!warning_given)
      std::cerr << "Second derivatives for Infinite elements"
		<< " are not yet implemented!"
		<< std::endl;

    d2phi.resize     (n_total_approx_shape_functions);
    d2phidx2.resize  (n_total_approx_shape_functions);
    d2phidxdy.resize (n_total_approx_shape_functions);
    d2phidxdz.resize (n_total_approx_shape_functions);
    d2phidy2.resize  (n_total_approx_shape_functions);
    d2phidydz.resize (n_total_approx_shape_functions);
    d2phidz2.resize  (n_total_approx_shape_functions);
    d2phidxi2.resize (n_total_approx_shape_functions);

    if (Dim > 1)
      {
        d2phidxideta.resize   (n_total_approx_shape_functions);
        d2phideta2.resize     (n_total_approx_shape_functions);
      }

    if (Dim > 2)
      {
        d2phidetadzeta.resize (n_total_approx_shape_functions);
        d2phidxidzeta.resize  (n_total_approx_shape_functions);
        d2phidzeta2.resize    (n_total_approx_shape_functions);
      }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
    
    if (Dim > 1)
      dphideta.resize      (n_total_approx_shape_functions);
    
    if (Dim == 3)
      dphidzeta.resize     (n_total_approx_shape_functions);
      
    
    phi_map.resize         (n_total_mapping_shape_functions);
    dphidxi_map.resize     (n_total_mapping_shape_functions);
#ifdef ENABLE_SECOND_DERIVATIVES
    d2phidxi2_map.resize   (n_total_mapping_shape_functions);

    if (Dim > 1)
      {
        d2phidxideta_map.resize   (n_total_mapping_shape_functions);
        d2phideta2_map.resize     (n_total_mapping_shape_functions);
      }

    if (Dim == 3)
      {
        d2phidxidzeta_map.resize  (n_total_mapping_shape_functions);
        d2phidetadzeta_map.resize (n_total_mapping_shape_functions);
        d2phidzeta2_map.resize    (n_total_mapping_shape_functions);
      }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
    
    if (Dim > 1)
      dphideta_map.resize  (n_total_mapping_shape_functions);
    
    if (Dim == 3)
      dphidzeta_map.resize (n_total_mapping_shape_functions);
  }



  // -----------------------------------------------------------------
  // collect all the for loops, where inner vectors are 
  // resized to the appropriate number of quadrature points
  {    
    for (unsigned int i=0; i<n_total_approx_shape_functions; i++)
      {
	phi[i].resize         (n_total_qp);
	dphi[i].resize        (n_total_qp);
	dphidx[i].resize      (n_total_qp);
	dphidy[i].resize      (n_total_qp);
	dphidz[i].resize      (n_total_qp);
	dphidxi[i].resize     (n_total_qp);
#ifdef ENABLE_SECOND_DERIVATIVES
	d2phi[i].resize       (n_total_qp);
	d2phidx2[i].resize    (n_total_qp);
	d2phidxdy[i].resize   (n_total_qp);
	d2phidxdz[i].resize   (n_total_qp);
	d2phidy2[i].resize    (n_total_qp);
	d2phidydz[i].resize   (n_total_qp);
	d2phidy2[i].resize    (n_total_qp);
	d2phidxi2[i].resize   (n_total_qp);

	if (Dim > 1)
	  {
	    d2phidxideta[i].resize   (n_total_qp);
	    d2phideta2[i].resize     (n_total_qp);
	  }
	if (Dim > 2)	     
	  {
	    d2phidxidzeta[i].resize  (n_total_qp);
	    d2phidetadzeta[i].resize (n_total_qp);
	    d2phidzeta2[i].resize    (n_total_qp);
	  }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	   
	if (Dim > 1)
	  dphideta[i].resize  (n_total_qp);
	   
	if (Dim == 3)	     
	  dphidzeta[i].resize (n_total_qp);
	     
      }
       
    for (unsigned int i=0; i<n_total_mapping_shape_functions; i++)
      {
	phi_map[i].resize         (n_total_qp);
	dphidxi_map[i].resize     (n_total_qp);
#ifdef ENABLE_SECOND_DERIVATIVES
	d2phidxi2_map[i].resize   (n_total_qp);
	if (Dim > 1)
	  {
	    d2phidxideta_map[i].resize   (n_total_qp);
	    d2phideta2_map[i].resize     (n_total_qp);
	  }

	if (Dim > 2)
	  {
	    d2phidxidzeta_map[i].resize  (n_total_qp);
	    d2phidetadzeta_map[i].resize (n_total_qp);
	    d2phidzeta2_map[i].resize    (n_total_qp);
	  }
#endif // ifdef ENABLE_SECOND_DERIVATIVES
	   
	if (Dim > 1)
	  dphideta_map[i].resize  (n_total_qp);
	   
	if (Dim == 3)
	  dphidzeta_map[i].resize (n_total_qp);
      }
  }



  {
    // -----------------------------------------------------------------
    // (a) compute scalar values at _all_ quadrature points  -- for uniform 
    //     access from the outside to these fields
    // (b) form a std::vector<Real> which contains the appropriate weights
    //     of the combined quadrature rule!
    const std::vector<Point>&  radial_qp = radial_qrule->get_points();
    libmesh_assert (radial_qp.size() == n_radial_qp);

    const std::vector<Real>&   radial_qw = radial_qrule->get_weights();
    const std::vector<Real>&   base_qw   = base_qrule->get_weights();
    libmesh_assert (radial_qw.size() == n_radial_qp);
    libmesh_assert (base_qw.size()   == n_base_qp);

    for (unsigned int rp=0; rp<n_radial_qp; rp++)
      for (unsigned int bp=0; bp<n_base_qp; bp++)
        {
          weight   [ bp+rp*n_base_qp ] = Radial::D       (radial_qp[rp](0)); 
	  dweightdv[ bp+rp*n_base_qp ] = Radial::D_deriv (radial_qp[rp](0));

	  _total_qrule_weights[  bp+rp*n_base_qp ] = radial_qw[rp] * base_qw[bp];
	}
  }

 
  /**
   * Stop logging the radial shape function initialization
   */
  STOP_LOG("init_shape_functions()", "InfFE");

}





template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::combine_base_radial(const Elem* inf_elem)
{
  libmesh_assert (inf_elem != NULL);
  // at least check whether the base element type is correct.
  // otherwise this version of computing dist would give problems
  libmesh_assert (base_elem->type() == Base::get_elem_type(inf_elem->type()));


  /**
   * Start logging the combination of radial and base parts
   */
  START_LOG("combine_base_radial()", "InfFE");


  // zero  the phase, since it is to be summed up
  std::fill (dphasedxi.begin(),   dphasedxi.end(),   0.);
  std::fill (dphasedeta.begin(),  dphasedeta.end(),  0.);
  std::fill (dphasedzeta.begin(), dphasedzeta.end(), 0.);


  const unsigned int n_base_mapping_sf   = dist.size();
  const Point origin                     = inf_elem->origin();

  // for each new infinite element, compute the radial distances
  for (unsigned int n=0; n<n_base_mapping_sf; n++)
      dist[n] =  Point(base_elem->point(n) - origin).size();


  switch (Dim)
    {

      //------------------------------------------------------------
      // 1D
    case 1:
      {
	std::cout << "ERROR: Not implemented." << std::endl;
	libmesh_error();
	break;
      }


      
      //------------------------------------------------------------
      // 2D
    case 2:
      {
	std::cout << "ERROR: Not implemented." << std::endl;
	libmesh_error();
     	break;
      }


      
      //------------------------------------------------------------
      // 3D
    case 3:
      {
	// fast access to the approximation and mapping shapes of base_fe
	const std::vector<std::vector<Real> >& S  = base_fe->phi;
	const std::vector<std::vector<Real> >& Ss = base_fe->dphidxi;
	const std::vector<std::vector<Real> >& St = base_fe->dphideta;
	const std::vector<std::vector<Real> >& S_map  = base_fe->phi_map;
	const std::vector<std::vector<Real> >& Ss_map = base_fe->dphidxi_map;
	const std::vector<std::vector<Real> >& St_map = base_fe->dphideta_map;

	const unsigned int n_radial_qp         = radial_qrule->n_points();
	const unsigned int n_base_qp           = base_qrule->  n_points();

	const unsigned int n_total_mapping_sf  = radial_map.size() * n_base_mapping_sf;

	const unsigned int n_total_approx_sf   = Radial::n_dofs(fe_type.radial_order) *  base_fe->n_shape_functions();


	// compute the phase term derivatives
	{
  	  unsigned int tp=0;
	  for (unsigned int rp=0; rp<n_radial_qp; rp++)  // over radial qp's
	    for (unsigned int bp=0; bp<n_base_qp; bp++)  // over base qp's
	      {
	        // sum over all base shapes, to get the average distance
    	        for (unsigned int i=0; i<n_base_mapping_sf; i++)
	          {
	            dphasedxi[tp]   += Ss_map[i][bp] * dist[i] * radial_map   [1][rp];
		    dphasedeta[tp]  += St_map[i][bp] * dist[i] * radial_map   [1][rp];
		    dphasedzeta[tp] += S_map [i][bp] * dist[i] * dradialdv_map[1][rp];
		  }

		tp++;

	      } // loop radial and base qp's

	}




	libmesh_assert (phi.size()       == n_total_approx_sf);
	libmesh_assert (dphidxi.size()   == n_total_approx_sf);
	libmesh_assert (dphideta.size()  == n_total_approx_sf);
	libmesh_assert (dphidzeta.size() == n_total_approx_sf);
	
	// compute the overall approximation shape functions,
	// pick the appropriate radial and base shapes through using
	// _base_shape_index and _radial_shape_index
	for (unsigned int rp=0; rp<n_radial_qp; rp++)  // over radial qp's
	  for (unsigned int bp=0; bp<n_base_qp; bp++)  // over base qp's
	    for (unsigned int ti=0; ti<n_total_approx_sf; ti++)  // over _all_ approx_sf
	      {
		// let the index vectors take care of selecting the appropriate base/radial shape
		const unsigned int bi = _base_shape_index  [ti];
		const unsigned int ri = _radial_shape_index[ti];
		phi      [ti][bp+rp*n_base_qp] = S [bi][bp] * mode[ri][rp] * som[rp];
		dphidxi  [ti][bp+rp*n_base_qp] = Ss[bi][bp] * mode[ri][rp] * som[rp];
		dphideta [ti][bp+rp*n_base_qp] = St[bi][bp] * mode[ri][rp] * som[rp];
		dphidzeta[ti][bp+rp*n_base_qp] = S [bi][bp] 
		    * (dmodedv[ri][rp] * som[rp] + mode[ri][rp] * dsomdv[rp]);
	      }

	
	libmesh_assert (phi_map.size()       == n_total_mapping_sf);
	libmesh_assert (dphidxi_map.size()   == n_total_mapping_sf);
	libmesh_assert (dphideta_map.size()  == n_total_mapping_sf);
	libmesh_assert (dphidzeta_map.size() == n_total_mapping_sf);
	
	// compute the overall mapping functions,
	// pick the appropriate radial and base entries through using
	// _base_node_index and _radial_node_index
	for (unsigned int rp=0; rp<n_radial_qp; rp++)  // over radial qp's
	  for (unsigned int bp=0; bp<n_base_qp; bp++)  // over base qp's
	    for (unsigned int ti=0; ti<n_total_mapping_sf; ti++)  // over all mapping shapes
	      {
		// let the index vectors take care of selecting the appropriate base/radial mapping shape
		const unsigned int bi = _base_node_index  [ti];
		const unsigned int ri = _radial_node_index[ti];
		phi_map      [ti][bp+rp*n_base_qp] = S_map [bi][bp] * radial_map   [ri][rp];
		dphidxi_map  [ti][bp+rp*n_base_qp] = Ss_map[bi][bp] * radial_map   [ri][rp];
		dphideta_map [ti][bp+rp*n_base_qp] = St_map[bi][bp] * radial_map   [ri][rp];
		dphidzeta_map[ti][bp+rp*n_base_qp] = S_map [bi][bp] * dradialdv_map[ri][rp];
	      }
			

	break;
      }


    default:
      libmesh_error();
    }


  /**
   * Start logging the combination of radial and base parts
   */
  STOP_LOG("combine_base_radial()", "InfFE");

}






template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
void InfFE<Dim,T_radial,T_map>::compute_shape_functions(const Elem*)
{
  libmesh_assert (radial_qrule != NULL);


  
  // Start logging the overall computation of shape functions
  START_LOG("compute_shape_functions()", "InfFE");


  const unsigned int n_total_qp  = _n_total_qp;


  //-------------------------------------------------------------------------
  // Compute the shape function values (and derivatives)
  // at the Quadrature points.  Note that the actual values
  // have already been computed via init_shape_functions

  // Compute the value of the derivative shape function i at quadrature point p
  switch (dim)
    {
      
    case 1:
      {
	std::cout << "ERROR: Not implemented." << std::endl;
	libmesh_error();
	break;
      }

    case 2:
      {
	std::cout << "ERROR: Not implemented." << std::endl;
	libmesh_error();
	break;
      }
    
    case 3:
      {
	// These are _all_ shape functions of this infinite element
	for (unsigned int i=0; i<phi.size(); i++)
	  for (unsigned int p=0; p<n_total_qp; p++)
	    {
	      // dphi/dx    = (dphi/dxi)*(dxi/dx) + (dphi/deta)*(deta/dx) + (dphi/dzeta)*(dzeta/dx);
	      dphi[i][p](0) =
		dphidx[i][p] = (dphidxi[i][p]*dxidx_map[p] +
				dphideta[i][p]*detadx_map[p] +
				dphidzeta[i][p]*dzetadx_map[p]);
		
	      // dphi/dy    = (dphi/dxi)*(dxi/dy) + (dphi/deta)*(deta/dy) + (dphi/dzeta)*(dzeta/dy);
	      dphi[i][p](1) =
		dphidy[i][p] = (dphidxi[i][p]*dxidy_map[p] +
				dphideta[i][p]*detady_map[p] +
				dphidzeta[i][p]*dzetady_map[p]);
		
	      // dphi/dz    = (dphi/dxi)*(dxi/dz) + (dphi/deta)*(deta/dz) + (dphi/dzeta)*(dzeta/dz);
	      dphi[i][p](2) =
		dphidz[i][p] = (dphidxi[i][p]*dxidz_map[p] +
				dphideta[i][p]*detadz_map[p] +
				dphidzeta[i][p]*dzetadz_map[p]);	      
	    }


	// This is the derivative of the phase term of this infinite element
	for (unsigned int p=0; p<n_total_qp; p++)
	  {
	    // the derivative of the phase term
	    dphase[p](0) = (dphasedxi[p]   * dxidx_map[p] +
			    dphasedeta[p]  * detadx_map[p] +
			    dphasedzeta[p] * dzetadx_map[p]);
		
	    dphase[p](1) = (dphasedxi[p]   * dxidy_map[p] +
			    dphasedeta[p]  * detady_map[p] +
			    dphasedzeta[p] * dzetady_map[p]);
		
	    dphase[p](2) = (dphasedxi[p]   * dxidz_map[p] +
			    dphasedeta[p]  * detadz_map[p] +
			    dphasedzeta[p] * dzetadz_map[p]);

	    // the derivative of the radial weight - varies only in radial direction,
	    // therefore dweightdxi = dweightdeta = 0.
	    dweight[p](0) = dweightdv[p] * dzetadx_map[p];
		
	    dweight[p](1) = dweightdv[p] * dzetady_map[p];
		
	    dweight[p](2) = dweightdv[p] * dzetadz_map[p];

	  }

	break;
      }



    default:
      {
	libmesh_error();
      }
    }


  
  // Stop logging the overall computation of shape functions
  STOP_LOG("compute_shape_functions()", "InfFE");

}



template <unsigned int Dim, FEFamily T_radial, InfMapType T_map>
bool InfFE<Dim,T_radial,T_map>::shapes_need_reinit() const 
{ 
  return false; 
}


//--------------------------------------------------------------
// Explicit instantiations
#include "inf_fe_instantiate_1D.h"
#include "inf_fe_instantiate_2D.h"
#include "inf_fe_instantiate_3D.h"



#endif //ifdef ENABLE_INFINITE_ELEMENTS