fe_szabab_shape_2d.c

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

// $Id: fe_szabab_shape_2D.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



// C++ includes
#include <cmath> // for std::sqrt


// Local includes
#include "libmesh_config.h"

#ifdef ENABLE_HIGHER_ORDER_SHAPES

#include "fe.h"
#include "elem.h"
#include "utility.h"


// Anonymous namespace to hold static std::sqrt values
namespace
{
  static const Real sqrt2  = std::sqrt(2.);
  static const Real sqrt6  = std::sqrt(6.);
  static const Real sqrt10 = std::sqrt(10.);
  static const Real sqrt14 = std::sqrt(14.);
  static const Real sqrt22 = std::sqrt(22.);
  static const Real sqrt26 = std::sqrt(26.);
}


template <>
Real FE<2,SZABAB>::shape(const ElemType,
			 const Order,
			 const unsigned int,
			 const Point&)
{
  std::cerr << "Szabo-Babuska polynomials require the element type\n"
	    << "because edge orientation is needed."
	    << std::endl;
  
  libmesh_error();
  return 0.;
}



template <>
Real FE<2,SZABAB>::shape(const Elem* elem,
			 const Order order,
			 const unsigned int i,
			 const Point& p)
{
  libmesh_assert (elem != NULL);
  
  const ElemType type = elem->type();

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

  // Declare that we are using our own special power function
  // from the Utility namespace.  This saves typing later.
  using Utility::pow;
  
  switch (totalorder)
    {      
      // 1st & 2nd-order Szabo-Babuska.
    case FIRST:
    case SECOND:
      {
	switch (type)
	  {

	    // Szabo-Babuska shape functions on the triangle.
	  case TRI6:
	    {
	      const Real l1 = 1-p(0)-p(1);
	      const Real l2 = p(0);
	      const Real l3 = p(1);	      
	      
	      libmesh_assert (i<6);
	      
	      switch (i)
		{
		case 0: return l1;		
		case 1: return l2;		  
		case 2: return l3;
		  
		case 3: return l1*l2*(-4.*sqrt6);			  
		case 4: return l2*l3*(-4.*sqrt6);		  
		case 5: return l3*l1*(-4.*sqrt6);	
		  
		default:
		  libmesh_error();
		}
	    }
    

	    // Szabo-Babuska shape functions on the quadrilateral.
	  case QUAD8:
	  case QUAD9:
	    {
	      // Compute quad shape functions as a tensor-product
	      const Real xi  = p(0);
	      const Real eta = p(1);
	      
	      libmesh_assert (i < 9);
	      
	      //                                0  1  2  3  4  5  6  7  8
	      static const unsigned int i0[] = {0, 1, 1, 0, 2, 1, 2, 0, 2};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 2, 1, 2, 2};
	      
	      return (FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
		      FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
	      
	    }

	    
	  default:
	    libmesh_error();
	  }
      }
	   

      // 3rd-order Szabo-Babuska.
    case THIRD:
      {
	switch (type)
	  {
	 
	    // Szabo-Babuska shape functions on the triangle.
	  case TRI6:
	    {
	      Real l1 = 1-p(0)-p(1);
	      Real l2 = p(0);
	      Real l3 = p(1);
	      
	      Real f=1;
	      
	      libmesh_assert (i<10);


	      if (i==4 && (elem->point(0) > elem->point(1)))f=-1;
	      if (i==6 && (elem->point(1) > elem->point(2)))f=-1;     
	      if (i==8 && (elem->point(2) > elem->point(0)))f=-1;


	      switch (i)
		{
		  //nodal modes
		case 0: return l1;		
		case 1: return l2;		  
		case 2: return l3;
		  
		  //side modes
		case 3: return   l1*l2*(-4.*sqrt6);			  
		case 4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);

		case 5: return   l2*l3*(-4.*sqrt6);		  
		case 6: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
	  
		case 7: return   l3*l1*(-4.*sqrt6);
		case 8: return f*l3*l1*(-4.*sqrt10)*(l1-l3);	
		
		  //internal modes
		case 9: return l1*l2*l3;
	
		default:
		  libmesh_error();
		}
	    }
	    

	    // Szabo-Babuska shape functions on the quadrilateral.
	  case QUAD8:
	  case QUAD9:
	    {
	      // Compute quad shape functions as a tensor-product
	      const Real xi  = p(0);
	      const Real eta = p(1);
	      
	      libmesh_assert (i < 16);

	      //                                0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15
	      static const unsigned int i0[] = {0,  1,  1,  0,  2,  3,  1,  1,  2,  3,  0,  0,  2,  3,  2,  3};
	      static const unsigned int i1[] = {0,  0,  1,  1,  0,  0,  2,  3,  1,  1,  2,  3,  2,  2,  3,  3};
	      	
	      Real f=1.;
	      	      
	      // take care of edge orientation, this is needed at
	      // edge shapes with (y=0)-asymmetric 1D shapes, these have
	      // one 1D shape index being 0 or 1, the other one being odd and >=3
	      	     
	      switch(i)
		{
		case  5: // edge 0 points    			
		  if (elem->point(0) > elem->point(1))f = -1.;
		  break;
	      	case  7: // edge 1 points
		  if (elem->point(1) > elem->point(2))f = -1.;
		  break;
	        case  9: // edge 2 points
		  if (elem->point(3) > elem->point(2))f = -1.;
		  break;
	        case 11: // edge 3 points
		  if (elem->point(0) > elem->point(3))f = -1.;
		  break;
		}	      
	      
	      return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
			FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
	    }

	  default:
	    libmesh_error();
	  }
      }
	   

      

      // 4th-order Szabo-Babuska.
    case FOURTH:
      {
	switch (type)
	  {
	    // Szabo-Babuska shape functions on the triangle.
	  case TRI6:
	    {
	      Real l1 = 1-p(0)-p(1);
	      Real l2 = p(0);
	      Real l3 = p(1);
	      
	      Real f=1;
	      
	      libmesh_assert (i<15);
	      
	      
	      if (i== 4 && (elem->point(0) > elem->point(1)))f=-1;
	      if (i== 7 && (elem->point(1) > elem->point(2)))f=-1;     
	      if (i==10 && (elem->point(2) > elem->point(0)))f=-1;
	      

	      switch (i)
		{
		  //nodal modes
		case  0: return l1;		
		case  1: return l2;		  
		case  2: return l3;
		  
		  //side modes
		case  3: return   l1*l2*(-4.*sqrt6);			  
		case  4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
		case  5: return   l1*l2*(-sqrt14)*(5.*pow<2>(l2-l1)-1);		  
		  
		case  6: return   l2*l3*(-4.*sqrt6);	
		case  7: return f*l2*l3*(-4.*sqrt10)*(l3-l2);	  
		case  8: return   l2*l3*(-sqrt14)*(5.*pow<2>(l3-l2)-1);	  
		  
		case  9: return   l3*l1*(-4.*sqrt6);		  
		case 10: return f*l3*l1*(-4.*sqrt10)*(l1-l3);		
		case 11: return   l3*l1*(-sqrt14)*(5.*pow<2>(l1-l3)-1);
		  
		  //internal modes
		case 12: return l1*l2*l3;
		  
		case 13: return l1*l2*l3*(l2-l1);
		case 14: return l1*l2*l3*(2*l3-1);
	
		default:
		  libmesh_error();
		}
	    }
	  

	    // Szabo-Babuska shape functions on the quadrilateral.
	  case QUAD8:
	  case QUAD9:
	    {
	      // Compute quad shape functions as a tensor-product
	      const Real xi  = p(0);
	      const Real eta = p(1);
	      
	      libmesh_assert (i < 25);
	      
	      //                                0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
	      static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
	      
	      Real f=1.;	      	      
	      
	      switch(i)
		{
		case  5: // edge 0 points    			
		  if (elem->point(0) > elem->point(1))f = -1.;
		  break;
	      	case  8: // edge 1 points
		  if (elem->point(1) > elem->point(2))f = -1.;
		  break;
	        case 11: // edge 2 points
		  if (elem->point(3) > elem->point(2))f = -1.;
		  break;
	        case 14: // edge 3 points
		  if (elem->point(0) > elem->point(3))f = -1.;
		  break;
		}	      
   	      
	      return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
			FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));
	    }
	    
	  default:
	    libmesh_error();
	  }
      }
      

      

      // 5th-order Szabo-Babuska.
    case FIFTH:
      {
	switch (type)
	  {
	    // Szabo-Babuska shape functions on the triangle.
	  case TRI6:
	    {
	      Real l1 = 1-p(0)-p(1);
	      Real l2 = p(0);
	      Real l3 = p(1);

	      const Real x=l2-l1;
	      const Real y=2.*l3-1;	
	      
	      Real f=1;	      
	      
	      libmesh_assert (i<21);      	      
	      

	      if ((i== 4||i== 6) && (elem->point(0) > elem->point(1)))f=-1;
	      if ((i== 8||i==10) && (elem->point(1) > elem->point(2)))f=-1;     
	      if ((i==12||i==14) && (elem->point(2) > elem->point(0)))f=-1;
	      

	      switch (i)
		{
		  //nodal modes
		case  0: return l1;		
		case  1: return l2;		  
		case  2: return l3;
		  
		  //side modes
		case  3: return   l1*l2*(-4.*sqrt6);			  
		case  4: return f*l1*l2*(-4.*sqrt10)*(l2-l1);
		case  5: return   -sqrt14*l1*l2*(5.0*l1*l1-1.0+(-10.0*l1+5.0*l2)*l2);
		case  6: return f*(-sqrt2)*l1*l2*((9.-21.*l1*l1)*l1+(-9.+63.*l1*l1+(-63.*l1+21.*l2)*l2)*l2);
	  
		case  7: return   l2*l3*(-4.*sqrt6);
		case  8: return f*l2*l3*(-4.*sqrt10)*(l3-l2);
		case  9: return   -sqrt14*l2*l3*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l2)*l2);
		case 10: return -f*sqrt2*l2*l3*((-9.0+21.0*l3*l3)*l3+(-63.0*l3*l3+9.0+(63.0*l3-21.0*l2)*l2)*l2);
		  
		case 11: return   l3*l1*(-4.*sqrt6);		  
		case 12: return f*l3*l1*(-4.*sqrt10)*(l1-l3); 
		case 13: return -sqrt14*l3*l1*(5.0*l3*l3-1.0+(-10.0*l3+5.0*l1)*l1);
		case 14: return f*(-sqrt2)*l3*l1*((9.0-21.0*l3*l3)*l3+(-9.0+63.0*l3*l3+(-63.0*l3+21.0*l1)*l1)*l1);
		  
		  //internal modes
		case 15: return l1*l2*l3;
		  
		case 16: return l1*l2*l3*x;	
		case 17: return l1*l2*l3*y;
		  
		case 18: return l1*l2*l3*(1.5*l1*l1-.5+(-3.0*l1+1.5*l2)*l2);
		case 19: return l1*l2*l3*x*y;	
		case 20: return l1*l2*l3*(1.0+(-6.0+6.0*l3)*l3);
		  
		default:
		  libmesh_error();
		}
	    } // case TRI6
	    
	    // Szabo-Babuska shape functions on the quadrilateral.
	  case QUAD8:
	  case QUAD9:
	    {
	      // Compute quad shape functions as a tensor-product
	      const Real xi  = p(0);
	      const Real eta = p(1);
	      
	      libmesh_assert (i < 36);
	      
	      //                                0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
	      static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 0, 0, 0, 0, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5, 2, 3, 4, 5};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 2, 3, 4, 5, 1, 1, 1, 1, 2, 3, 4, 5, 2, 2, 2, 2, 3, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5};
	      
	      Real f=1.;
	      
	      switch(i)
		{
		case  5: // edge 0 points
		case  7:    			
		  if (elem->point(0) > elem->point(1))f = -1.;
		  break;
		case  9: // edge 1 points
	      	case 11:	      	
		  if (elem->point(1) > elem->point(2))f = -1.;
		  break;
	        case 13: // edge 2 points
	        case 15:
		  if (elem->point(3) > elem->point(2))f = -1.;
		  break;
	        case 14: // edge 3 points
	        case 19:
		  if (elem->point(0) > elem->point(3))f = -1.;
		  break;
		}	     
	      
	      return f*(FE<1,SZABAB>::shape(EDGE3, totalorder, i0[i], xi)*
			FE<1,SZABAB>::shape(EDGE3, totalorder, i1[i], eta));	      

	    } // case QUAD8/QUAD9

	  default:
	    libmesh_error();

	  } // switch type

      } // case FIFTH
      
      // 6th-order Szabo-Babuska.
    case SIXTH:
      {
	switch (type)
	  {
	    // Szabo-Babuska shape functions on the triangle.
	  case TRI6:
	    {
	      Real l1 = 1-p(0)-p(1);
	      Real l2 = p(0);
	      Real l3 = p(1);

	      const Real x=l2-l1;
	      const Real y=2.*l3-1;	      
	      
	      Real f=1;

	      libmesh_assert (i<28);
	      
	      
	      if ((i== 4||i== 6) && (elem->point(0) > elem->point(1)))f=-1;
	      if ((i== 9||i==11) && (elem->point(1) > elem->point(2)))f=-1;     
	      if ((i==14||i==16) && (elem->point(2) > elem->point(0)))f=-1;

	      
	      switch (i)
		{