fe_szabab_shape_2d.c

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

	  case TRI6:
	    {
	      // Here we use finite differences to compute the derivatives!
	      const Real eps = 1.e-6;
	      
	      libmesh_assert (i < 15);
	      libmesh_assert (j < 2);
	      
	      switch (j)
		{
		  //  d()/dxi
		case 0:
		  {
		    const Point pp(p(0)+eps, p(1));
		    const Point pm(p(0)-eps, p(1));

		    return (FE<2,SZABAB>::shape(elem, order, i, pp) -
			    FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
		  }

		  // d()/deta
		case 1:
		  {
		    const Point pp(p(0), p(1)+eps);
		    const Point pm(p(0), p(1)-eps);

		    return (FE<2,SZABAB>::shape(elem, order, i, pp) -
			    FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
		  }
		  

		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;
		}	      

	      
	      switch (j)
		{
		  // d()/dxi
		case 0:		  		  
		  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
			    FE<1,SZABAB>::shape      (EDGE3, totalorder, i1[i],    eta));
	      
		  // d()/deta
		case 1:		  		  
		  return f*(FE<1,SZABAB>::shape      (EDGE3, totalorder, i0[i],    xi)*
			    FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));

		default:
		  libmesh_error();
		}
	    }

	  default:
	    libmesh_error();
	  }
      }
	   

      

      // 5th-order Szabo-Babuska.
    case FIFTH:
      {
	// Szabo-Babuska shape functions on the quadrilateral.
	switch (type)
	  {
	 
	    // Szabo-Babuska shape functions on the triangle.
	  case TRI6:
	    {
	      // Here we use finite differences to compute the derivatives!
	      const Real eps = 1.e-6;
	      
	      libmesh_assert (i < 21);
	      libmesh_assert (j < 2);
	      
	      switch (j)
		{
		  //  d()/dxi
		case 0:
		  {
		    const Point pp(p(0)+eps, p(1));
		    const Point pm(p(0)-eps, p(1));

		    return (FE<2,SZABAB>::shape(elem, order, i, pp) -
			    FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
		  }

		  // d()/deta
		case 1:
		  {
		    const Point pp(p(0), p(1)+eps);
		    const Point pm(p(0), p(1)-eps);

		    return (FE<2,SZABAB>::shape(elem, order, i, pp) -
			    FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
		  }
		  

		default:
		  libmesh_error();
		}
	    }
	
	    

	  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;
		}	     
	      
	      
	      switch (j)
		{
		  // d()/dxi
		case 0:		  		  
		  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
			    FE<1,SZABAB>::shape      (EDGE3, totalorder, i1[i],    eta));
	      
		  // d()/deta
		case 1:		  		  
		  return f*(FE<1,SZABAB>::shape      (EDGE3, totalorder, i0[i],    xi)*
			    FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));

		default:
		  libmesh_error();
		}
	    }

	  default:
	    libmesh_error();
	  }
      }


    // 6th-order Szabo-Babuska.
    case SIXTH:
      {
	// Szabo-Babuska shape functions on the quadrilateral.
	switch (type)
	  {
	 
	    // Szabo-Babuska shape functions on the triangle.
	  case TRI6:
	    {
	      // Here we use finite differences to compute the derivatives!
	      const Real eps = 1.e-6;
	      
	      libmesh_assert (i < 28);
	      libmesh_assert (j < 2);
	      
	      switch (j)
		{
		  //  d()/dxi
		case 0:
		  {
		    const Point pp(p(0)+eps, p(1));
		    const Point pm(p(0)-eps, p(1));

		    return (FE<2,SZABAB>::shape(elem, order, i, pp) -
			    FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
		  }

		  // d()/deta
		case 1:
		  {
		    const Point pp(p(0), p(1)+eps);
		    const Point pm(p(0), p(1)-eps);

		    return (FE<2,SZABAB>::shape(elem, order, i, pp) -
			    FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
		  }
		  

		default:
		  libmesh_error();
		}
	    }
	
	    

	  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 < 49);

	      //                                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 36 37 38 39 40 41 42 43 44 45 46 47 48
	      static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6, 2, 3, 4, 5, 6};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6};
	      
	      Real f=1.;
	     
	      switch(i)
		{
	      	case  5: // edge 0 points
	      	case  7:    			
		  if (elem->point(0) > elem->point(1))f = -1.;
		  break;
	      	case 10: // edge 1 points
	      	case 12:	      	
		  if (elem->point(1) > elem->point(2))f = -1.;
		  break;
	        case 15: // edge 2 points
	        case 17:
		  if (elem->point(3) > elem->point(2))f = -1.;
		  break;
	        case 20: // edge 3 points
	        case 22:
		  if (elem->point(0) > elem->point(3))f = -1.;
		  break;
		}	       

	      
	      switch (j)
		{
		  // d()/dxi
		case 0:		  		  
		  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
			    FE<1,SZABAB>::shape      (EDGE3, totalorder, i1[i],    eta));
	      
		  // d()/deta
		case 1:		  		  
		  return f*(FE<1,SZABAB>::shape      (EDGE3, totalorder, i0[i],    xi)*
			    FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));

		default:
		  libmesh_error();
		}
	    }

	  default:
	    libmesh_error();
	  }
      }


    // 7th-order Szabo-Babuska.
    case SEVENTH:
      {
	// Szabo-Babuska shape functions on the quadrilateral.
	switch (type)
	  {
	 
	    // Szabo-Babuska shape functions on the triangle.
	  case TRI6:
	    {
	      // Here we use finite differences to compute the derivatives!
	      const Real eps = 1.e-6;
	      
	      libmesh_assert (i < 36);
	      libmesh_assert (j < 2);
	      
	      switch (j)
		{
		  //  d()/dxi
		case 0:
		  {
		    const Point pp(p(0)+eps, p(1));
		    const Point pm(p(0)-eps, p(1));

		    return (FE<2,SZABAB>::shape(elem, order, i, pp) -
			    FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
		  }

		  // d()/deta
		case 1:
		  {
		    const Point pp(p(0), p(1)+eps);
		    const Point pm(p(0), p(1)-eps);

		    return (FE<2,SZABAB>::shape(elem, order, i, pp) -
			    FE<2,SZABAB>::shape(elem, order, i, pm))/2./eps;
		  }
		  

		default:
		  libmesh_error();
		}
	    }
	
	    

	  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 < 64);

	      //                                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 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
	      static const unsigned int i0[] = {0, 1, 1, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7, 2, 3, 4, 5, 6, 7};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 6, 7, 1, 1, 1, 1, 1, 1, 2, 3, 4, 5, 6, 7, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7};
	      
	      Real f=1.;
	     
	      switch(i)
		{
	      	case  5: // edge 0 points
	      	case  7:
		case  9:    			
		  if (elem->point(0) > elem->point(1))f = -1.;
		  break;
	      	case 11: // edge 1 points
	      	case 13:	      	
		case 15:
		  if (elem->point(1) > elem->point(2))f = -1.;
		  break;
	        case 17: // edge 2 points
	        case 19:
		case 21:
		  if (elem->point(3) > elem->point(2))f = -1.;
		  break;
	        case 23: // edge 3 points
	        case 25:
		case 27:
		  if (elem->point(0) > elem->point(3))f = -1.;
		  break;
		}	     	       
	      
	      
	      switch (j)
		{
		  // d()/dxi
		case 0:		  		  
		  return f*(FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
			    FE<1,SZABAB>::shape      (EDGE3, totalorder, i1[i],    eta));
		  
		  // d()/deta
		case 1:		  		  
		  return f*(FE<1,SZABAB>::shape      (EDGE3, totalorder, i0[i],    xi)*
			    FE<1,SZABAB>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta));
		  
		default:
		  libmesh_error();
		}
	    }
	    
	  default:
	    libmesh_error();
	  }
      }
      

      
      // by default throw an error;call the orientation-independent shape functions
    default:
      std::cerr << "ERROR: Unsupported polynomial order!" << std::endl;
      libmesh_error();
    }

  
  libmesh_error();
  return 0.;
}



template <>
Real FE<2,SZABAB>::shape_second_deriv(const ElemType,
			              const Order,
			              const unsigned int,
			              const unsigned int,
			              const Point&)
{
  static bool warning_given = false;

  if (!warning_given)
  std::cerr << "Second derivatives for Szabab elements "
            << " are not yet implemented!"
            << std::endl;

  warning_given = true;
  return 0.;
}



template <>
Real FE<2,SZABAB>::shape_second_deriv(const Elem*,
			              const Order,
			              const unsigned int,
			              const unsigned int,
			              const Point&)
{
  static bool warning_given = false;

  if (!warning_given)
  std::cerr << "Second derivatives for Szabab elements "
            << " are not yet implemented!"
            << std::endl;

  warning_given = true;
  return 0.;
}


#endif	// ENABLE_HIGHER_ORDER_SHAPES