fe_interface.c

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

	    
	  case HIERARCHIC:
	    FE<3,HIERARCHIC>::dofs_on_side(elem, o, s, di);
            return;
	    
	  case LAGRANGE:
	    FE<3,LAGRANGE>::dofs_on_side(elem, o, s, di);
            return;
	    
	  case MONOMIAL:
	    FE<3,MONOMIAL>::dofs_on_side(elem, o, s, di);
            return;

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<3,BERNSTEIN>::dofs_on_side(elem, o, s, di);
            return;

	  case SZABAB:
	    FE<3,SZABAB>::dofs_on_side(elem, o, s, di);
            return;

#endif
	    
	  case XYZ:
	    FEXYZ<3>::dofs_on_side(elem, o, s, di);
            return;

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  libmesh_error();
}



void FEInterface::dofs_on_edge(const Elem* const elem,
			       const unsigned int dim,
			       const FEType& fe_t,
			       unsigned int e,
			       std::vector<unsigned int>& di)
{
  const Order o = fe_t.order;

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<1,BERNSTEIN>::dofs_on_edge(elem, o, e, di);
            return;

	  case SZABAB:
	    FE<1,SZABAB>::dofs_on_edge(elem, o, e, di);
            return;

#endif
	    
	  case XYZ:
	    FEXYZ<1>::dofs_on_edge(elem, o, e, di);
            return;

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<2,BERNSTEIN>::dofs_on_edge(elem, o, e, di);
            return;

	  case SZABAB:
	    FE<2,SZABAB>::dofs_on_edge(elem, o, e, di);
            return;

#endif
	    
	  case XYZ:
	    FEXYZ<2>::dofs_on_edge(elem, o, e, di);
            return;

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<3,BERNSTEIN>::dofs_on_edge(elem, o, e, di);
            return;

	  case SZABAB:
	    FE<3,SZABAB>::dofs_on_edge(elem, o, e, di);
            return;

#endif
	    
	  case XYZ:
	    FEXYZ<3>::dofs_on_edge(elem, o, e, di);
            return;

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  libmesh_error();
}




void FEInterface::nodal_soln(const unsigned int dim,
			     const FEType& fe_t,
			     const Elem* elem,
			     const std::vector<Number>& elem_soln,
			     std::vector<Number>&       nodal_soln)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(elem->type()) )
  {
    ifem_nodal_soln(dim, fe_t, elem, elem_soln, nodal_soln);
    return;
  }

#endif

  const Order order = fe_t.order;
  
  switch (dim)
    {
      // 1D
    case 1:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    FE<1,CLOUGH>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;
	    
	  case HERMITE:
	    FE<1,HERMITE>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;
	    
	  case HIERARCHIC:
	    FE<1,HIERARCHIC>::nodal_soln(elem, order,
					 elem_soln, nodal_soln);
	    return;

	  case LAGRANGE:
	    FE<1,LAGRANGE>::nodal_soln(elem, order,
				       elem_soln, nodal_soln);
	    return;
   
	  case MONOMIAL:
	    FE<1,MONOMIAL>::nodal_soln(elem, order,
				       elem_soln, nodal_soln);
	    return;

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<1,BERNSTEIN>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;

	  case SZABAB:
	    FE<1,SZABAB>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;

#endif

	  case XYZ:
	    FEXYZ<1>::nodal_soln(elem, order,
				 elem_soln, nodal_soln);
	    return;

	  default:
	    libmesh_error();
	  }
      }

      
      // 2D
    case 2:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    FE<2,CLOUGH>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;
	    
	  case HERMITE:
	    FE<2,HERMITE>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;
	    
	  case HIERARCHIC:
	    FE<2,HIERARCHIC>::nodal_soln(elem, order,
					 elem_soln, nodal_soln);
	    return;
	    
	  case LAGRANGE:
	    FE<2,LAGRANGE>::nodal_soln(elem, order,
				       elem_soln, nodal_soln);
	    return;
   
	  case MONOMIAL:
	    FE<2,MONOMIAL>::nodal_soln(elem, order,
				       elem_soln, nodal_soln);
	    return;

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<2,BERNSTEIN>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;

	  case SZABAB:
	    FE<2,SZABAB>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;

#endif

	  case XYZ:
	    FEXYZ<2>::nodal_soln(elem, order,
				 elem_soln, nodal_soln);
	    return;

	  default:
	    libmesh_error();
	  }
      }

      
      // 3D
    case 3:
      {
	switch (fe_t.family)
	  {
	  case HERMITE:
	    FE<3,HERMITE>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;
	    
	  case HIERARCHIC:
	    FE<3,HIERARCHIC>::nodal_soln(elem, order,
					 elem_soln, nodal_soln);
	    return;
	    
	  case LAGRANGE:
	    FE<3,LAGRANGE>::nodal_soln(elem, order,
				       elem_soln, nodal_soln);
	    return;
	    
	  case MONOMIAL:
	    FE<3,MONOMIAL>::nodal_soln(elem, order,
				       elem_soln, nodal_soln);
	    return;

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<3,BERNSTEIN>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;

	  case SZABAB:
	    FE<3,SZABAB>::nodal_soln(elem, order,
				     elem_soln, nodal_soln);
	    return;

#endif

	  case XYZ:
	    FEXYZ<3>::nodal_soln(elem, order,
				 elem_soln, nodal_soln);
	    return;

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }
}




Point FEInterface::inverse_map (const unsigned int dim,
				const FEType& fe_t,
				const Elem* elem,
				const Point& p,
				const Real tolerance,
				const bool secure)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(elem->type()) )
    return ifem_inverse_map(dim, fe_t, elem, p,tolerance, secure);

#endif

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<1,BERNSTEIN>::inverse_map(elem, p, tolerance, secure);

	  case SZABAB:
	    return FE<1,SZABAB>::inverse_map(elem, p, tolerance, secure);

#endif
	  case XYZ:
	    return FEXYZ<1>::inverse_map(elem, p, tolerance, secure);


	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<2,BERNSTEIN>::inverse_map(elem, p, tolerance, secure);

	  case SZABAB:
	    return FE<2,SZABAB>::inverse_map(elem, p, tolerance, secure);

#endif
	  case XYZ:
	    return FEXYZ<2>::inverse_map(elem, p, tolerance, secure);

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<3,BERNSTEIN>::inverse_map(elem, p, tolerance, secure);

	  case SZABAB:
	    return FE<3,SZABAB>::inverse_map(elem, p, tolerance, secure);

#endif
	  case XYZ:
	    return FEXYZ<3>::inverse_map(elem, p, tolerance, secure);

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  
  libmesh_error();
  Point pt;
  return pt;
}




void FEInterface::inverse_map (const unsigned int dim,
			       const FEType& fe_t,
			       const Elem* elem,
			       const std::vector<Point>& physical_points,
			       std::vector<Point>&       reference_points,
			       const Real tolerance,
			       const bool secure)
{
  const unsigned int n_pts = physical_points.size();

  // Resize the vector
  reference_points.resize(n_pts);
  
  if (n_pts == 0)
    {
      std::cerr << "WARNING: empty vector physical_points!"
		<< std::endl;
      here();
      return;
    }
  



  
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(elem->type()) )
    {
      ifem_inverse_map(dim, fe_t, elem, physical_points, reference_points, tolerance, secure);
      return;

//       std::cerr << "ERROR: Not implemented!"
// 		<< std::endl;
//       libmesh_error();
    }
  
#endif

  switch (dim)
    {
      // 1D
    case 1:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    FE<1,CLOUGH>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case HERMITE:
	    FE<1,HERMITE>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case HIERARCHIC:
	    FE<1,HIERARCHIC>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case LAGRANGE:
	    FE<1,LAGRANGE>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case MONOMIAL:
	    FE<1,MONOMIAL>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<1,BERNSTEIN>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;

	  case SZABAB:
	    FE<1,SZABAB>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;

#endif
	    
	  case XYZ:
	    FEXYZ<1>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;

	  default:
	    libmesh_error();
	  }
      }

      
      // 2D
    case 2:
      {
	switch (fe_t.family)
	  {
	  case CLOUGH:
	    FE<2,CLOUGH>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case HERMITE:
	    FE<2,HERMITE>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case HIERARCHIC:
	    FE<2,HIERARCHIC>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case LAGRANGE:
	    FE<2,LAGRANGE>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case MONOMIAL:
	    FE<2,MONOMIAL>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<2,BERNSTEIN>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;

	  case SZABAB:
	    FE<2,SZABAB>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;

#endif
	    
	  case XYZ:
	    FEXYZ<2>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;

	  default:
	    libmesh_error();
	  }
      }

      
      // 3D
    case 3:
      {
	switch (fe_t.family)
	  {
	  case HERMITE:
	    FE<3,HERMITE>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case HIERARCHIC:
	    FE<3,HIERARCHIC>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case LAGRANGE:
	    FE<3,LAGRANGE>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  case MONOMIAL:
	    FE<3,MONOMIAL>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    FE<3,BERNSTEIN>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;

	  case SZABAB:
	    FE<3,SZABAB>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;

#endif
	    
	  case XYZ:
	    FEXYZ<3>::inverse_map(elem, physical_points, reference_points, tolerance, secure);
	    return;
	    
	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  
  libmesh_error();
  return;
}



bool FEInterface::on_reference_element(const Point& p,
				       const ElemType t,
				       const Real eps)
{
  return FEBase::on_reference_element(p,t,eps);
}