fe_xyz.c

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

	    return 15;

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

      
    default:
      {
        const unsigned int order = static_cast<unsigned int>(o);
	switch (t)
	  {
	  case EDGE2:
	  case EDGE3:
	    return (order+1);

	  case TRI3:
	  case TRI6:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    return (order+1)*(order+2)/2;

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



template <unsigned int Dim, FEFamily T>
unsigned int FE<Dim,T>::n_dofs_at_node(const ElemType,
				       const Order,
				       const unsigned int)
{
  // Monomials elements have no dofs at nodes
  // (just on the element)
  return 0;
}



template <unsigned int Dim, FEFamily T>
unsigned int FE<Dim,T>::n_dofs_per_elem(const ElemType t,
					const Order o)
{
  switch (o)
    {
      // constant shape functions always have 1 DOF per element
    case CONSTANT:
      return 1;

      
      // Discontinuous linear shape functions
      // expressed in the XYZ monomials.
    case FIRST:
      {
	switch (t)
	  {
	    // 1D linears have 2 DOFs per element
	  case EDGE2:
	  case EDGE3:
	  case EDGE4:
	    return 2;

	    // 2D linears have 3 DOFs per element
	  case TRI3:
	  case TRI6:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    return 3;

	    // 3D linears have 4 DOFs per element
 	  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)
	  {
	    // 1D quadratics have 3 DOFs per element
	  case EDGE2:
	  case EDGE3:
	  case EDGE4:
	    return 3;

	    // 2D quadratics have 6 DOFs per element
	  case TRI3:
	  case TRI6:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    return 6;

	    // 3D quadratics have 10 DOFs per element
	  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:
	  case EDGE4:
	    return 5;

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

	  case TET4:
	  case TET10:
	  case HEX8:
	  case HEX20:
	  case HEX27:
	  case PRISM6:
	  case PRISM15:
	  case PRISM18:
	  case PYRAMID5:
	    return 35;
	    
	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for " << o << "th order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
	  }
      }
      
    default:
      {
        const unsigned int order = static_cast<unsigned int>(o);
	switch (t)
	  {
	  case EDGE2:
	  case EDGE3:
	    return (order+1);

	  case TRI3:
	  case TRI6:
	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    return (order+1)*(order+2)/2;

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



template <unsigned int Dim, FEFamily T>
FEContinuity FE<Dim,T>::get_continuity() const
{
  return DISCONTINUOUS;
}



template <unsigned int Dim, FEFamily T>
bool FE<Dim,T>::is_hierarchic() const
{
  return true;
}



#ifdef ENABLE_AMR
template <unsigned int Dim, FEFamily T>
void FE<Dim,T>::compute_constraints (DofConstraints &,
				     DofMap &,
				     const unsigned int,
				     const Elem*)
{
  // Monomials are discontinuous...  No constraints.
  return;
}
#endif // #ifdef ENABLE_AMR



template <unsigned int Dim, FEFamily T>
bool FE<Dim,T>::shapes_need_reinit() const
{
  return false;
}



//--------------------------------------------------------------
// Explicit instantiation of member functions
INSTANTIATE_MBRF(1,XYZ);
INSTANTIATE_MBRF(2,XYZ);
INSTANTIATE_MBRF(3,XYZ);
template void  FEXYZ<1>::init_shape_functions(const std::vector<Point>&,
					    const Elem*);
template void  FEXYZ<2>::init_shape_functions(const std::vector<Point>&,
					    const Elem*);
template void  FEXYZ<3>::init_shape_functions(const std::vector<Point>&,
					    const Elem*);
template void  FEXYZ<1>::compute_shape_functions(const Elem*);
template void  FEXYZ<2>::compute_shape_functions(const Elem*);
template void  FEXYZ<3>::compute_shape_functions(const Elem*);

#ifdef ENABLE_AMR
template void FE<2,XYZ>::compute_constraints(DofConstraints&, DofMap&, 
					     const unsigned int, const Elem*);
template void FE<3,XYZ>::compute_constraints(DofConstraints&, DofMap&, 
					     const unsigned int, const Elem*);
#endif // #ifdef ENABLE_AMR