fe_interface.c

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

Real FEInterface::shape(const unsigned int dim,
			const FEType& fe_t,
			const ElemType t,
			const unsigned int i,
			const Point& p)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(t) )
    return ifem_shape(dim, fe_t, t, i, p);

#endif

  const Order o = fe_t.order;
  
  switch (dim)
    {
      // 1D
    case 1:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    return FE<1,CLOUGH>::shape(t,o,i,p);
	    
	  case HERMITE:
	    return FE<1,HERMITE>::shape(t,o,i,p);
	    
	  case HIERARCHIC:
	    return FE<1,HIERARCHIC>::shape(t,o,i,p);
	    
	  case LAGRANGE:
	    return FE<1,LAGRANGE>::shape(t,o,i,p);
	    
	  case MONOMIAL:
	    return FE<1,MONOMIAL>::shape(t,o,i,p);

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<1,BERNSTEIN>::shape(t,o,i,p);

	  case SZABAB:
	    return FE<1,SZABAB>::shape(t,o,i,p);

#endif
	    
	  case XYZ:
	    return FEXYZ<1>::shape(t,o,i,p);

	  default:
	    libmesh_error();
	  }
      }

      
      // 2D
    case 2:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    return FE<2,CLOUGH>::shape(t,o,i,p);
	    
	  case HERMITE:
	    return FE<2,HERMITE>::shape(t,o,i,p);
	    
	  case HIERARCHIC:
	    return FE<2,HIERARCHIC>::shape(t,o,i,p);
	    
	  case LAGRANGE:
	    return FE<2,LAGRANGE>::shape(t,o,i,p);
	    
	  case MONOMIAL:
	    return FE<2,MONOMIAL>::shape(t,o,i,p);

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<2,BERNSTEIN>::shape(t,o,i,p);

	  case SZABAB:
	    return FE<2,SZABAB>::shape(t,o,i,p);

#endif
	    
	  case XYZ:
	    return FEXYZ<2>::shape(t,o,i,p);

	  default:
	    libmesh_error();
	  }
      }

      
      // 3D
    case 3:
      {
	switch (fe_t.family)
	  {
	  case HERMITE:
	    return FE<3,HERMITE>::shape(t,o,i,p);
	    
	  case HIERARCHIC:
	    return FE<3,HIERARCHIC>::shape(t,o,i,p);
	    
	  case LAGRANGE:
	    return FE<3,LAGRANGE>::shape(t,o,i,p);
	    
	  case MONOMIAL:
	    return FE<3,MONOMIAL>::shape(t,o,i,p);

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<3,BERNSTEIN>::shape(t,o,i,p);

	  case SZABAB:
	    return FE<3,SZABAB>::shape(t,o,i,p);

#endif
	    
	  case XYZ:
	    return FEXYZ<3>::shape(t,o,i,p);

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  
  libmesh_error();
  return 0.;
}




Real FEInterface::shape(const unsigned int dim,
			const FEType& fe_t,
			const Elem* elem,
			const unsigned int i,
			const Point& p)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(elem->type()) )
    return ifem_shape(dim, fe_t, elem, i, p);

#endif

  const Order o = fe_t.order;

  switch (dim)
    {
      // 1D
    case 1:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    return FE<1,CLOUGH>::shape(elem,o,i,p);
	    
	  case HERMITE:
	    return FE<1,HERMITE>::shape(elem,o,i,p);
	    
	  case HIERARCHIC:
	    return FE<1,HIERARCHIC>::shape(elem,o,i,p);
	    
	  case LAGRANGE:
	    return FE<1,LAGRANGE>::shape(elem,o,i,p);
	    
	  case MONOMIAL:
	    return FE<1,MONOMIAL>::shape(elem,o,i,p);

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<1,BERNSTEIN>::shape(elem,o,i,p);

	  case SZABAB:
	    return FE<1,SZABAB>::shape(elem,o,i,p);

#endif
	    
	  case XYZ:
	    return FEXYZ<1>::shape(elem,o,i,p);

	  default:
	    libmesh_error();
	  }
      }

      
      // 2D
    case 2:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    return FE<2,CLOUGH>::shape(elem,o,i,p);
	    
	  case HERMITE:
	    return FE<2,HERMITE>::shape(elem,o,i,p);
	    
	  case HIERARCHIC:
	    return FE<2,HIERARCHIC>::shape(elem,o,i,p);
	    
	  case LAGRANGE:
	    return FE<2,LAGRANGE>::shape(elem,o,i,p);
	    
	  case MONOMIAL:
	    return FE<2,MONOMIAL>::shape(elem,o,i,p);

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<2,BERNSTEIN>::shape(elem,o,i,p);

	  case SZABAB:
	    return FE<2,SZABAB>::shape(elem,o,i,p);

#endif
	    
	  case XYZ:
	    return FEXYZ<2>::shape(elem,o,i,p);

	  default:
	    libmesh_error();
	  }
      }

      
      // 3D
    case 3:
      {
	switch (fe_t.family)
	  {
	  case HERMITE:
	    return FE<3,HERMITE>::shape(elem,o,i,p);
	    
	  case HIERARCHIC:
	    return FE<3,HIERARCHIC>::shape(elem,o,i,p);
	    
	  case LAGRANGE:
	    return FE<3,LAGRANGE>::shape(elem,o,i,p);
	    
	  case MONOMIAL:
	    return FE<3,MONOMIAL>::shape(elem,o,i,p);

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<3,BERNSTEIN>::shape(elem,o,i,p);

	  case SZABAB:
	    return FE<3,SZABAB>::shape(elem,o,i,p);

#endif
	    
	  case XYZ:
	    return FEXYZ<3>::shape(elem,o,i,p);

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  
  libmesh_error();
  return 0.;
}




void FEInterface::compute_data(const unsigned int dim,
			       const FEType& fe_t,
			       const Elem* elem,
			       FEComputeData& data)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(elem->type()) )
    {
      data.init();
      ifem_compute_data(dim, fe_t, elem, data);
      return;
    }

#endif

  FEType p_refined = fe_t;
  p_refined.order = static_cast<Order>(p_refined.order + elem->p_level());

  const unsigned int n_dof = n_dofs (dim, p_refined, elem->type());
  const Point&       p     = data.p;
  data.shape.resize(n_dof);

  // set default values for all the output fields
  data.init();

  for (unsigned int n=0; n<n_dof; n++)
      data.shape[n] = shape(dim, p_refined, elem, n, p);

   return;
}



#ifdef ENABLE_AMR

void FEInterface::compute_constraints (DofConstraints &constraints,
				       DofMap &dof_map,
				       const unsigned int variable_number,
				       const Elem* elem)
{
  libmesh_assert (elem != NULL);

  const FEType& fe_t = dof_map.variable_type(variable_number);
  
  switch (elem->dim())
    {
    case 1:
      {
	// No constraints in 1D.
	return;
      }

      
    case 2:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    FE<2,CLOUGH>::compute_constraints (constraints,
					       dof_map,
					       variable_number,
					       elem); return;

	  case HERMITE:
	    FE<2,HERMITE>::compute_constraints (constraints,
					        dof_map,
					        variable_number,
					        elem); return;

	  case LAGRANGE:
	    FE<2,LAGRANGE>::compute_constraints (constraints,
						 dof_map,
						 variable_number,
						 elem); return;

	  case HIERARCHIC:
	    FE<2,HIERARCHIC>::compute_constraints (constraints,
						   dof_map,
						   variable_number,
						   elem); return;

	  default:
	    return;
	  }
      }


    case 3:
      {
	switch (fe_t.family)
	  {
	  case HERMITE:
	    FE<3,HERMITE>::compute_constraints (constraints,
					        dof_map,
					        variable_number,
					        elem); return;

	  case LAGRANGE:
	    FE<3,LAGRANGE>::compute_constraints (constraints,
					         dof_map,
						 variable_number,
						 elem); return;      

	  case HIERARCHIC:
	    FE<3,HIERARCHIC>::compute_constraints (constraints,
						   dof_map,
						   variable_number,
						   elem); return;
	  default:
	    return;
	  }
      }

      
    default:
      libmesh_error();
    }
}

#endif // #ifdef ENABLE_AMR



#ifdef ENABLE_PERIODIC

void FEInterface::compute_periodic_constraints (DofConstraints &constraints,
				                DofMap &dof_map,
                                                PeriodicBoundaries &boundaries,
						const MeshBase &mesh,
				                const unsigned int variable_number,
				                const Elem* elem)
{
  // No element-specific optimizations currently exist
  FEBase::compute_periodic_constraints (constraints,
                                        dof_map,
                                        boundaries,
                                        mesh,
				        variable_number,
				        elem);
}

#endif // #ifdef ENABLE_PERIODIC
  


unsigned int FEInterface::max_order(const FEType& fe_t,
			            const ElemType& el_t)
{
  // Yeah, I know, infinity is much larger than 11, but our
  // solvers don't seem to like high degree polynomials, and our
  // quadrature rules and number_lookups tables
  // need to go up higher.
  const unsigned int unlimited = 11;

  // If we used 0 as a default, then elements missing from this 
  // table (e.g. infinite elements) would be considered broken.
  const unsigned int unknown = unlimited;

  switch (fe_t.family)
    {
      case LAGRANGE:
	switch (el_t)
	  {
	    case EDGE2:
	    case EDGE3:
	    case EDGE4:
	      return 3;
	    case TRI3:
	      return 1;
	    case TRI6:
	      return 2;
	    case QUAD4:
	      return 1;
	    case QUAD8:
	    case QUAD9:
	      return 2;
	    case TET4:
	      return 1;
	    case TET10:
	      return 2;
	    case HEX8:
	      return 1;
	    case HEX20:
	    case HEX27:
	      return 2;
	    case PRISM6:
	    case PRISM15:
	      return 1;
	    case PRISM18:
	      return 2;
	    case PYRAMID5:
	      return 1;
	    default:
	      return unknown;
	  }
	break;
      case MONOMIAL:
	switch (el_t)
	  {
	    case EDGE2:
	    case EDGE3:
	    case EDGE4:
	    case TRI3:
	    case TRI6:
	    case QUAD4:
	    case QUAD8:
	    case QUAD9:
	    case TET4:
	    case TET10:
	    case HEX8:
	    case HEX20:
	    case HEX27:
	    case PRISM6:
	    case PRISM15:
	    case PRISM18:
	    case PYRAMID5:
	      return unlimited;
	    default:
	      return unknown;
	  }
	break;
#ifdef ENABLE_HIGHER_ORDER_SHAPES
      case BERNSTEIN:
	switch (el_t)
	  {
	    case EDGE2:
	    case EDGE3:
	    case EDGE4:
	      return unlimited;
	    case TRI3:
	      return 0;
	    case TRI6:
	      return 6;
	    case QUAD4:
	      return 0;
	    case QUAD8:
	    case QUAD9:
	      return unlimited;
	    case TET4:
	      return 1;
	    case TET10:
	      return 2;
	    case HEX8:
	      return 0;
	    case HEX20:
	      return 2;
	    case HEX27:
	      return 4;
	    case PRISM6:
	    case PRISM15:
	    case PRISM18:
	    case PYRAMID5:
	      return 0;
	    default:
	      return unknown;
	  }
	break;
      case SZABAB:
	switch (el_t)
	  {
	    case EDGE2:
	    case EDGE3:
	    case EDGE4:
	      return 7;
	    case TRI3:
	      return 0;
	    case TRI6:
	      return 7;
	    case QUAD4:
	      return 0;
	    case QUAD8:
	    case QUAD9:
	      return 7;
	    case TET4:
	    case TET10:
	    case HEX8:
	    case HEX20:
	    case HEX27:
	    case PRISM6:
	    case PRISM15:
	    case PRISM18:
	    case PYRAMID5:
	      return 0;
	    default:
	      return unknown;
	  }
	break;
#endif
      case XYZ:
	switch (el_t)
	  {
	    case EDGE2:
	    case EDGE3:
	    case EDGE4:
	    case TRI3:
	    case TRI6:
	    case QUAD4:
	    case QUAD8:
	    case QUAD9:
	    case TET4:
	    case TET10:
	    case HEX8:
	    case HEX20:
	    case HEX27:
	    case PRISM6:
	    case PRISM15:
	    case PRISM18:
	    case PYRAMID5:
	      return unlimited;
	    default:
	      return unknown;
	  }
	break;
      case CLOUGH:
	switch (el_t)
	  {
	    case EDGE2:
	    case EDGE3:
	      return 3;
	    case EDGE4:
	    case TRI3:
	      return 0;
	    case TRI6:
	      return 3;
	    case QUAD4:
	    case QUAD8:
	    case QUAD9:
	    case TET4:
	    case TET10:
	    case HEX8:
	    case HEX20:
	    case HEX27:
	    case PRISM6:
	    case PRISM15:
	    case PRISM18:
	    case PYRAMID5:
	      return 0;
	    default:
	      return unknown;
	  }
	break;
      case HERMITE:
	switch (el_t)
	  {
	    case EDGE2:
	    case EDGE3:
	      return unlimited;
	    case EDGE4:
	    case TRI3:
	    case TRI6:
	      return 0;
	    case QUAD4:
	      return 3;
	    case QUAD8:
	    case QUAD9:
	      return unlimited;
	    case TET4:
	    case TET10:
	      return 0;
	    case HEX8:
	      return 3;
	    case HEX20:
	    case HEX27:
	      return unlimited;
	    case PRISM6:
	    case PRISM15:
	    case PRISM18:
	    case PYRAMID5:
	      return 0;
	    default:
	      return unknown;
	  }
	break;
      case HIERARCHIC:
	switch (el_t)
	  {
	    case EDGE2:
	    case EDGE3:
	    case EDGE4:
	      return unlimited;
	    case TRI3:
	      return 1;
	    case TRI6:
	      return unlimited;
	    case QUAD4:
	      return 1;
	    case QUAD8:
	    case QUAD9:
	      return unlimited;
	    case TET4:
	    case TET10:
	      return 0;
	    case HEX8:
	    case HEX20:
	      return 1;
	    case HEX27:
	      return unlimited;
	    case PRISM6:
	    case PRISM15:
	    case PRISM18:
	    case PYRAMID5:
	      return 0;
	    default:
	      return unknown;
	  }
	break;
      default:
	return 0;
	break;
    }
}

  

bool FEInterface::extra_hanging_dofs(const FEType& fe_t)
{
  switch (fe_t.family)
    {
      case LAGRANGE:
      case MONOMIAL:
#ifdef ENABLE_HIGHER_ORDER_SHAPES
      case BERNSTEIN:
      case SZABAB:
#endif
      case XYZ:
	return false;
      case CLOUGH:
      case HERMITE:
      case HIERARCHIC:
      default:
	return true;
    }
}