fe_lagrange.c

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

// $Id: fe_lagrange.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 "dof_map.h"
#include "fe.h"
#include "fe_macro.h"
#include "fe_interface.h"
#include "elem.h"
#include "threads.h"




// ------------------------------------------------------------
// Lagrange-specific implementations
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();

  const Order totalorder = static_cast<Order>(order+elem->p_level());

  nodal_soln.resize(n_nodes);


  
  switch (totalorder)
    {
      // linear Lagrange shape functions
    case FIRST:
      {
	switch (type)
	  {
	  case EDGE3:
	    {
	      libmesh_assert (elem_soln.size()  == 2);
	      libmesh_assert (nodal_soln.size() == 3);

	      nodal_soln[0] = elem_soln[0];
	      nodal_soln[1] = elem_soln[1];
	      nodal_soln[2] = .5*(elem_soln[0] + elem_soln[1]);

	      return;
	    }

          case EDGE4:
            {
              libmesh_assert(elem_soln.size() == 2);
              libmesh_assert(nodal_soln.size() == 4);

              nodal_soln[0] = elem_soln[0];
              nodal_soln[1] = elem_soln[1];
              nodal_soln[2] = (2.*elem_soln[0] + elem_soln[1])/3.;
              nodal_soln[3] = (elem_soln[0] + 2.*elem_soln[1])/3.;

              return;
            }

	    
	  case TRI6:
	    {
	      libmesh_assert (elem_soln.size()  == 3);
	      libmesh_assert (nodal_soln.size() == 6);
	      
	      nodal_soln[0] = elem_soln[0];
	      nodal_soln[1] = elem_soln[1];
	      nodal_soln[2] = elem_soln[2];
	      nodal_soln[3] = .5*(elem_soln[0] + elem_soln[1]);	      
	      nodal_soln[4] = .5*(elem_soln[1] + elem_soln[2]);
	      nodal_soln[5] = .5*(elem_soln[2] + elem_soln[0]);

	      return;
	    }
	    

	  case QUAD8:
	  case QUAD9:
	    {
	      libmesh_assert (elem_soln.size()  == 4);
	      
	      if (type == QUAD8)
		libmesh_assert (nodal_soln.size() == 8);
	      else
		libmesh_assert (nodal_soln.size() == 9);

	      
	      nodal_soln[0] = elem_soln[0];
	      nodal_soln[1] = elem_soln[1];
	      nodal_soln[2] = elem_soln[2];
	      nodal_soln[3] = elem_soln[3];
	      nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
	      nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
	      nodal_soln[6] = .5*(elem_soln[2] + elem_soln[3]);
	      nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
					       
	      if (type == QUAD9)
		nodal_soln[8] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);

	      return;
	    }

	    
	  case TET10:
	    {
	      libmesh_assert (elem_soln.size()  == 4);
	      libmesh_assert (nodal_soln.size() == 10);
	      
	      nodal_soln[0] = elem_soln[0];
	      nodal_soln[1] = elem_soln[1];
	      nodal_soln[2] = elem_soln[2];
	      nodal_soln[3] = elem_soln[3];
	      nodal_soln[4] = .5*(elem_soln[0] + elem_soln[1]);
	      nodal_soln[5] = .5*(elem_soln[1] + elem_soln[2]);
	      nodal_soln[6] = .5*(elem_soln[2] + elem_soln[0]);
	      nodal_soln[7] = .5*(elem_soln[3] + elem_soln[0]);
	      nodal_soln[8] = .5*(elem_soln[3] + elem_soln[1]);
	      nodal_soln[9] = .5*(elem_soln[3] + elem_soln[2]);
	      
	      return;
	    }

	    
	  case HEX20:
	  case HEX27:
	    {
	      libmesh_assert (elem_soln.size()  == 8);
	      
	      if (type == HEX20)
		libmesh_assert (nodal_soln.size() == 20);
	      else
		libmesh_assert (nodal_soln.size() == 27);
	      
	      nodal_soln[0]  = elem_soln[0];
	      nodal_soln[1]  = elem_soln[1];
	      nodal_soln[2]  = elem_soln[2];
	      nodal_soln[3]  = elem_soln[3];
	      nodal_soln[4]  = elem_soln[4];
	      nodal_soln[5]  = elem_soln[5];
	      nodal_soln[6]  = elem_soln[6];
	      nodal_soln[7]  = elem_soln[7];
	      nodal_soln[8]  = .5*(elem_soln[0] + elem_soln[1]);
	      nodal_soln[9]  = .5*(elem_soln[1] + elem_soln[2]);
	      nodal_soln[10] = .5*(elem_soln[2] + elem_soln[3]);
	      nodal_soln[11] = .5*(elem_soln[3] + elem_soln[0]);
	      nodal_soln[12] = .5*(elem_soln[0] + elem_soln[4]);
	      nodal_soln[13] = .5*(elem_soln[1] + elem_soln[5]);
	      nodal_soln[14] = .5*(elem_soln[2] + elem_soln[6]);
	      nodal_soln[15] = .5*(elem_soln[3] + elem_soln[7]);
	      nodal_soln[16] = .5*(elem_soln[4] + elem_soln[5]);
	      nodal_soln[17] = .5*(elem_soln[5] + elem_soln[6]);
	      nodal_soln[18] = .5*(elem_soln[6] + elem_soln[7]);
	      nodal_soln[19] = .5*(elem_soln[4] + elem_soln[7]);

	      if (type == HEX27)
		{
		  nodal_soln[20] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
		  nodal_soln[21] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[5]);
		  nodal_soln[22] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[6]);
		  nodal_soln[23] = .25*(elem_soln[2] + elem_soln[3] + elem_soln[6] + elem_soln[7]);
		  nodal_soln[24] = .25*(elem_soln[3] + elem_soln[0] + elem_soln[7] + elem_soln[4]);
		  nodal_soln[25] = .25*(elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);

		  nodal_soln[26] = .125*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3] +
					 elem_soln[4] + elem_soln[5] + elem_soln[6] + elem_soln[7]);
		}
	      
	      return;
	    }
	    
	    
	  case PRISM15:
	  case PRISM18:
	    {
	      libmesh_assert (elem_soln.size()  == 6);

	      if (type == PRISM15)
		libmesh_assert (nodal_soln.size() == 15);
	      else
		libmesh_assert (nodal_soln.size() == 18);
	      
	      nodal_soln[0]  = elem_soln[0];
	      nodal_soln[1]  = elem_soln[1];
	      nodal_soln[2]  = elem_soln[2];
	      nodal_soln[3]  = elem_soln[3];
	      nodal_soln[4]  = elem_soln[4];
	      nodal_soln[5]  = elem_soln[5];
	      nodal_soln[6]  = .5*(elem_soln[0] + elem_soln[1]);
	      nodal_soln[7]  = .5*(elem_soln[1] + elem_soln[2]);
	      nodal_soln[8]  = .5*(elem_soln[0] + elem_soln[2]);
	      nodal_soln[9]  = .5*(elem_soln[0] + elem_soln[3]);
	      nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
	      nodal_soln[11] = .5*(elem_soln[2] + elem_soln[5]);
	      nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
	      nodal_soln[13] = .5*(elem_soln[4] + elem_soln[5]);
	      nodal_soln[14] = .5*(elem_soln[3] + elem_soln[5]);

	      if (type == PRISM18)
		{
		  nodal_soln[15] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[4] + elem_soln[3]);
		  nodal_soln[16] = .25*(elem_soln[1] + elem_soln[2] + elem_soln[5] + elem_soln[4]);
		  nodal_soln[17] = .25*(elem_soln[2] + elem_soln[0] + elem_soln[3] + elem_soln[5]);
		}
	      	      
	      return;
	    }


	    
	  default:
	    {
	      // By default the element solution _is_ nodal,
	      // so just copy it.
	      nodal_soln = elem_soln;
	      
	      return;
	    }
	  }
      }

    case SECOND:
      {
	switch (type)
	  {
	  case EDGE4:
            {
              libmesh_assert(elem_soln.size()  == 3);
              libmesh_assert(nodal_soln.size() == 4);

              // Project quadratic solution onto cubic element nodes
              nodal_soln[0] = elem_soln[0];
              nodal_soln[1] = elem_soln[1];
              nodal_soln[2] = (2.*elem_soln[0] - elem_soln[1] + 
                               8.*elem_soln[2])/9.;
              nodal_soln[3] = (-elem_soln[0] + 2.*elem_soln[1] +
                               8.*elem_soln[2])/9.;
              return;
            }

          default:
            {
	      // By default the element solution _is_ nodal,
	      // so just copy it.
	      nodal_soln = elem_soln;
	      
	      return;
            }
          }
        }



      
    default:
      {
	// By default the element solution _is_ nodal,
	// so just copy it.
	nodal_soln = elem_soln;
	
	return;
      }      
    }
}




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

      // linear Lagrange shape functions
    case FIRST:
      {
	switch (t)
	  {
	  case NODEELEM:
	    return 1;

	  case EDGE2:
	  case EDGE3:
          case EDGE4:
	    return 2;

	  case TRI3:
	  case TRI6:
	    return 3;

	  case QUAD4:
	  case QUAD8:
	  case QUAD9:
	    return 4;

	  case TET4:
	  case TET10:
	    return 4;

	  case HEX8:
	  case HEX20:
	  case HEX27:
	    return 8;

	  case PRISM6:
	  case PRISM15:
	  case PRISM18:
	    return 6;

	  case PYRAMID5:
	    return 5;

	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for FIRST order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
	  }
      }

      
      // quadratic Lagrange shape functions
    case SECOND:
      {
	switch (t)
	  {
	  case NODEELEM:
	    return 1;

	  case EDGE3:
	    return 3;

	  case TRI6:
	    return 6;

	  case QUAD8:
	    return 8;
	    
	  case QUAD9:
	    return 9;

	  case TET10:
	    return 10;

	  case HEX20:
	    return 20;
	    
	  case HEX27:
	    return 27;

	  case PRISM15:
	    return 15;

	  case PRISM18:
	    return 18;

	  default:
	    {
#ifdef DEBUG
	      std::cerr << "ERROR: Bad ElemType = " << t
			<< " for SECOND order approximation!" 
			<< std::endl;
#endif
	      libmesh_error();	    
	    }
	  }
      }

    case THIRD:
      {
        switch (t)
        {
	  case NODEELEM:
	    return 1;

          case EDGE4: