fe_xyz.c

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

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

  // Start logging the shape function computation
  START_LOG("compute_shape_functions()", "FE");

  const std::vector<Point>& xyz_qp = this->get_xyz();
  
  // Compute the value of the derivative shape function i at quadrature point p
  switch (this->dim)
    {
      
    case 1:
      {
        if (this->calculate_phi)
	  for (unsigned int i=0; i<this->phi.size(); i++)
	    for (unsigned int p=0; p<this->phi[i].size(); p++)
	      this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz_qp[p]);
        if (this->calculate_dphi)
	  for (unsigned int i=0; i<this->dphi.size(); i++)
	    for (unsigned int p=0; p<this->dphi[i].size(); p++)
	      {
	        this->dphi[i][p](0) =
		  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
	      
	        this->dphi[i][p](1) = this->dphidy[i][p] = 0.;
	        this->dphi[i][p](2) = this->dphidz[i][p] = 0.;
	      }
#ifdef ENABLE_SECOND_DERIVATIVES
        if (this->calculate_d2phi)
	  for (unsigned int i=0; i<this->d2phi.size(); i++)
	    for (unsigned int p=0; p<this->d2phi[i].size(); p++)
	      {
	        this->d2phi[i][p](0,0) =
		  this->d2phidx2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
	      
#if DIM>1
	        this->d2phi[i][p](0,1) = this->d2phidxdy[i][p] =
	        this->d2phi[i][p](1,0) = 0.;
	        this->d2phi[i][p](1,1) = this->d2phidy2[i][p] = 0.;
#if DIM>2
	        this->d2phi[i][p](0,2) = this->d2phidxdz[i][p] =
	        this->d2phi[i][p](2,0) = 0.;
	        this->d2phi[i][p](1,2) = this->d2phidydz[i][p] =
	        this->d2phi[i][p](2,1) = 0.;
	        this->d2phi[i][p](2,2) = this->d2phidz2[i][p] = 0.;
#endif
#endif
	      }
#endif // ifdef ENABLE_SECOND_DERIVATIVES

	// All done
	break;
      }

    case 2:
      {
        if (this->calculate_phi)
	  for (unsigned int i=0; i<this->phi.size(); i++)
	    for (unsigned int p=0; p<this->phi[i].size(); p++)
	      this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz_qp[p]);
        if (this->calculate_dphi)
	  for (unsigned int i=0; i<this->dphi.size(); i++)
	    for (unsigned int p=0; p<this->dphi[i].size(); p++)
	      {
	        this->dphi[i][p](0) =
		  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
	      
	        this->dphi[i][p](1) =
		  this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz_qp[p]);
	      
#if DIM == 3  
	        this->dphi[i][p](2) = // can only assign to the Z component if DIM==3
#endif
		this->dphidz[i][p] = 0.;
	      }
#ifdef ENABLE_SECOND_DERIVATIVES
        if (this->calculate_d2phi)
	  for (unsigned int i=0; i<this->d2phi.size(); i++)
	    for (unsigned int p=0; p<this->d2phi[i].size(); p++)
	      {
	        this->d2phi[i][p](0,0) =
		  this->d2phidx2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
	      
	        this->d2phi[i][p](0,1) = this->d2phidxdy[i][p] =
	        this->d2phi[i][p](1,0) = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 1, xyz_qp[p]);
	        this->d2phi[i][p](1,1) = 
                  this->d2phidy2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 2, xyz_qp[p]);
#if DIM>2
	        this->d2phi[i][p](0,2) = this->d2phidxdz[i][p] =
	        this->d2phi[i][p](2,0) = 0.;
	        this->d2phi[i][p](1,2) = this->d2phidydz[i][p] =
	        this->d2phi[i][p](2,1) = 0.;
	        this->d2phi[i][p](2,2) = this->d2phidz2[i][p] = 0.;
#endif
	      }
#endif // ifdef ENABLE_SECOND_DERIVATIVES

	// All done
	break;
      }
    
    case 3:
      {
        if (this->calculate_dphi)
	  for (unsigned int i=0; i<this->phi.size(); i++)
	    for (unsigned int p=0; p<this->phi[i].size(); p++)
	      this->phi[i][p] = FE<Dim,XYZ>::shape (elem, this->fe_type.order, i, xyz_qp[p]);
	       
        if (this->calculate_dphi)
	  for (unsigned int i=0; i<this->dphi.size(); i++)
	    for (unsigned int p=0; p<this->dphi[i].size(); p++)
	      {
	        this->dphi[i][p](0) =
		  this->dphidx[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
		
	        this->dphi[i][p](1) =
		  this->dphidy[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 1, xyz_qp[p]);
		
	        this->dphi[i][p](2) =
		  this->dphidz[i][p] = FE<Dim,XYZ>::shape_deriv (elem, this->fe_type.order, i, 2, xyz_qp[p]);	      
	      }
#ifdef ENABLE_SECOND_DERIVATIVES
        if (this->calculate_d2phi)
	  for (unsigned int i=0; i<this->d2phi.size(); i++)
	    for (unsigned int p=0; p<this->d2phi[i].size(); p++)
	      {
	        this->d2phi[i][p](0,0) =
		  this->d2phidx2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 0, xyz_qp[p]);
	      
	        this->d2phi[i][p](0,1) = this->d2phidxdy[i][p] =
	        this->d2phi[i][p](1,0) = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 1, xyz_qp[p]);
	        this->d2phi[i][p](1,1) = 
                  this->d2phidy2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 2, xyz_qp[p]);
	        this->d2phi[i][p](0,2) = this->d2phidxdz[i][p] =
	        this->d2phi[i][p](2,0) = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 3, xyz_qp[p]);
	        this->d2phi[i][p](1,2) = this->d2phidydz[i][p] =
	        this->d2phi[i][p](2,1) = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 4, xyz_qp[p]);
	        this->d2phi[i][p](2,2) = this->d2phidz2[i][p] = FE<Dim,XYZ>::shape_second_deriv (elem, this->fe_type.order, i, 5, xyz_qp[p]);
	      }
#endif // ifdef ENABLE_SECOND_DERIVATIVES

	// All done
	break;
      }

    default:
      {
	libmesh_error();
      }
    }
  
  // Stop logging the shape function computation
  STOP_LOG("compute_shape_functions()", "FE");
}



template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::nodal_soln(const Elem* elem,
			   const Order order,
			   const std::vector<Number>& elem_soln,
			   std::vector<Number>&       nodal_soln)
{
  const unsigned int n_nodes = elem->n_nodes();
  
  const ElemType type = elem->type();

  nodal_soln.resize(n_nodes);

  const Order totalorder = static_cast<Order>(order + elem->p_level());
  
  switch (totalorder)
    {
      // Constant shape functions
    case CONSTANT:
      {
	libmesh_assert (elem_soln.size() == 1);
	
	const Number val = elem_soln[0];
	
	for (unsigned int n=0; n<n_nodes; n++)
	  nodal_soln[n] = val;
	
	return;
      }


      // For other bases do interpolation at the nodes
      // explicitly.
    default:
      {

	const unsigned int n_sf =
	  FE<Dim,T>::n_shape_functions(type, totalorder);
	
	for (unsigned int n=0; n<n_nodes; n++)
	  {
	    libmesh_assert (elem_soln.size() == n_sf);

	    // Zero before summation
	    nodal_soln[n] = 0;

	    // u_i = Sum (alpha_i phi_i)
	    for (unsigned int i=0; i<n_sf; i++)
	      nodal_soln[n] += elem_soln[i]*FE<Dim,T>::shape(elem,
							     order,
							     i,
							     elem->point(n));
	  }

	return;
      }
    }
}



template <unsigned int Dim, FEFamily T>
unsigned int FE<Dim,T>::n_dofs(const ElemType t, const Order o)
{
  switch (o)
    {

      // constant shape functions
      // no matter what shape there is only one DOF.
    case CONSTANT:
      return 1;


      // Discontinuous linear shape functions
      // expressed in the XYZ monomials.
    case FIRST:
      {
	switch (t)
	  {
	  case EDGE2:
	  case EDGE3:
	  case EDGE4:
	    return 2;

	  case TRI3:
	  case TRI6:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    return 3;

	  case TET4:
	  case TET10:
	  case HEX8:
	  case HEX20:
	  case HEX27:
	  case PRISM6:
	  case PRISM15:
	  case PRISM18:
	  case PYRAMID5:
	    return 4;
	    
	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for " << o << "th order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
	  }
      }


      // Discontinuous quadratic shape functions
      // expressed in the XYZ monomials.
    case SECOND:
      {
	switch (t)
	  {
	  case EDGE2:
	  case EDGE3:
	  case EDGE4:
	    return 3;

	  case TRI3:
	  case TRI6:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    return 6;

	  case TET4:
	  case TET10:
	  case HEX8:
	  case HEX20:
	  case HEX27:
	  case PRISM6:
	  case PRISM15:
	  case PRISM18:
	  case PYRAMID5:
	    return 10;
	    
	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for " << o << "th order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
	  }
      }


      // Discontinuous cubic shape functions
      // expressed in the XYZ monomials.
    case THIRD:
      {
	switch (t)
	  {
	  case EDGE2:
	  case EDGE3:
	  case EDGE4:
	    return 4;

	  case TRI3:
	  case TRI6:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    return 10;

	  case TET4:
	  case TET10:
	  case HEX8:
	  case HEX20:
	  case HEX27:
	  case PRISM6:
	  case PRISM15:
	  case PRISM18:
	  case PYRAMID5:
	    return 20;
	    
	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for " << o << "th order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
	  }
      }


      // Discontinuous quartic shape functions
      // expressed in the XYZ monomials.
    case FOURTH:
      {
	switch (t)
	  {
	  case EDGE2:
	  case EDGE3:
	    return 5;

	  case TRI3:
	  case TRI6:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9: