fe_interface.c

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

// $Id: fe_interface.C 2789 2008-04-13 02:24:40Z roystgnr $

// The libMesh Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA



// Local includes
#include "fe_interface.h"
#include "elem.h"
#include "fe.h"
#include "fe_compute_data.h"
#include "dof_map.h"

//------------------------------------------------------------
//FEInterface class members
FEInterface::FEInterface()
{
  std::cerr << "ERROR: Do not define an object of this type." 
	    << std::endl;
  libmesh_error();
}




unsigned int FEInterface::n_shape_functions(const unsigned int dim,
					    const FEType& fe_t,
					    const ElemType t)
{

#ifdef ENABLE_INFINITE_ELEMENTS
  /*
   * Since the FEType, stored in DofMap/(some System child), has to
   * be the _same_ for InfFE and FE, we have to catch calls
   * to infinite elements through the element type.
   */

  if ( is_InfFE_elem(t) )
    return ifem_n_shape_functions(dim, fe_t, t);

#endif

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<1,BERNSTEIN>::n_shape_functions(t, o);

	  case SZABAB:
	    return FE<1,SZABAB>::n_shape_functions(t, o);

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


	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<2,BERNSTEIN>::n_shape_functions(t, o);

	  case SZABAB:
	    return FE<2,SZABAB>::n_shape_functions(t, o);

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

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<3,BERNSTEIN>::n_shape_functions(t, o);

	  case SZABAB:
	    return FE<3,SZABAB>::n_shape_functions(t, o);

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

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  
  libmesh_error();
  return 0;
}





unsigned int FEInterface::n_dofs(const unsigned int dim,
				 const FEType& fe_t,
				 const ElemType t)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(t) )
    return ifem_n_dofs(dim, fe_t, t);

#endif

  const Order o = fe_t.order;

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<1,BERNSTEIN>::n_dofs(t, o);

	  case SZABAB:
	    return FE<1,SZABAB>::n_dofs(t, o);

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

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<2,BERNSTEIN>::n_dofs(t, o);

	  case SZABAB:
	    return FE<2,SZABAB>::n_dofs(t, o);

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

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<3,BERNSTEIN>::n_dofs(t, o);

	  case SZABAB:
	    return FE<3,SZABAB>::n_dofs(t, o);

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

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  
  libmesh_error();
  return 0;
}

		


unsigned int FEInterface::n_dofs_at_node(const unsigned int dim,
					 const FEType& fe_t,
					 const ElemType t,
					 const unsigned int n)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(t) )
    return ifem_n_dofs_at_node(dim, fe_t, t, n);

#endif

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<1,BERNSTEIN>::n_dofs_at_node(t, o, n);

	  case SZABAB:
	    return FE<1,SZABAB>::n_dofs_at_node(t, o, n);

#endif
	    
	  case XYZ:
	    return FEXYZ<1>::n_dofs_at_node(t, o, n);

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<2,BERNSTEIN>::n_dofs_at_node(t, o, n);

	  case SZABAB:
	    return FE<2,SZABAB>::n_dofs_at_node(t, o, n);

#endif
	    
	  case XYZ:
	    return FEXYZ<2>::n_dofs_at_node(t, o, n);

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<3,BERNSTEIN>::n_dofs_at_node(t, o, n);

	  case SZABAB:
	    return FE<3,SZABAB>::n_dofs_at_node(t, o, n);

#endif
	    
	  case XYZ:
	    return FEXYZ<3>::n_dofs_at_node(t, o, n);

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  
  libmesh_error();
  return 0;
}





unsigned int FEInterface::n_dofs_per_elem(const unsigned int dim,
					  const FEType& fe_t,
					  const ElemType t)
{
#ifdef ENABLE_INFINITE_ELEMENTS

  if ( is_InfFE_elem(t) )
    return ifem_n_dofs_per_elem(dim, fe_t, t);

#endif

  const Order o = fe_t.order;

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<1,BERNSTEIN>::n_dofs_per_elem(t, o);

	  case SZABAB:
	    return FE<1,SZABAB>::n_dofs_per_elem(t, o);

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

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<2,BERNSTEIN>::n_dofs_per_elem(t, o);

	  case SZABAB:
	    return FE<2,SZABAB>::n_dofs_per_elem(t, o);

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

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

	  case BERNSTEIN:
	    return FE<3,BERNSTEIN>::n_dofs_per_elem(t, o);

	  case SZABAB:
	    return FE<3,SZABAB>::n_dofs_per_elem(t, o);

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

	  default:
	    libmesh_error();
	  }
      }


    default:
      libmesh_error();
    }

  
  libmesh_error();
  return 0;
}




void FEInterface::dofs_on_side(const Elem* const elem,
			       const unsigned int dim,
			       const FEType& fe_t,
			       unsigned int s,
			       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_side(elem, o, s, di);
            return;
	    
	  case HERMITE:
	    FE<1,HERMITE>::dofs_on_side(elem, o, s, di);
            return;
	    
	  case HIERARCHIC:
	    FE<1,HIERARCHIC>::dofs_on_side(elem, o, s, di);
            return;
	    
	  case LAGRANGE:
	    FE<1,LAGRANGE>::dofs_on_side(elem, o, s, di);
            return;
	    
	  case MONOMIAL:
	    FE<1,MONOMIAL>::dofs_on_side(elem, o, s, di);
            return;

#ifdef ENABLE_HIGHER_ORDER_SHAPES

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

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

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

	  default:
	    libmesh_error();
	  }
      }

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

#ifdef ENABLE_HIGHER_ORDER_SHAPES

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

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

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

	  default:
	    libmesh_error();
	  }
      }

      
      // 3D
    case 3:
      {
	switch (fe_t.family)
	  {
	  case HERMITE:
	    FE<3,HERMITE>::dofs_on_side(elem, o, s, di);
            return;