fe_bernstein_shape_2d.c

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

// 		  libmesh_error();
// 		}
// 	    }
	    
	
	    
// 	    // Bernstein 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_reg[] = {0,  1,  1,  0,  2,  3,  1,  1,  2,  3,  0,  0,  2,  3,  2,  3};
// 	      static const unsigned int i1_reg[] = {0,  0,  1,  1,  0,  0,  2,  3,  1,  1,  2,  3,  2,  2,  3,  3};
	      
// 	      unsigned int i0=i0_reg[i];
// 	      unsigned int i1=i1_reg[i];

// 	      if((i== 4||i== 5) && elem->node(0) > elem->node(1)) i0=5-i0;
// 	      if((i== 6||i== 7) && elem->node(1) > elem->node(2)) i1=5-i1;
// 	      if((i== 8||i== 9) && elem->node(3) > elem->node(2)) i0=5-i0;
// 	      if((i==10||i==11) && elem->node(0) > elem->node(3)) i1=5-i1;      
	      
// 	      switch (j)
// 		{
// 			  // d()/dxi
// 		case 0:		  		  
// 		  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
// 			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1,    eta));
		  
// 		  // d()/deta
// 		case 1:		  		  
// 		  return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0,    xi)*
// 			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
		  
// 		default:
// 		  libmesh_error();
// 		}
// 	    }
	    
// 	  default:
// 	    libmesh_error();
// 	  }
//       }
      
      
      

//       // 4th-order Bernsteins.
//     case FOURTH:
//       {
// 	switch (type)
// 	  {
// 	    // Bernstein shape functions on the triangle.
// 	  case TRI6:
// 	    {
// 	      // I have been lazy here and am using finite differences
// 	      // to compute the derivatives!
// 	      const Real eps = 1.e-6;
	      
// 	      libmesh_assert (i < 15);
// 	      libmesh_assert (j < 2);
	      
// 	      unsigned int shape=i;
	      
// 	      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,BERNSTEIN>::shape(elem, order, shape, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, order, shape, 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,BERNSTEIN>::shape(elem, order, shape, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
// 		  }
		  
		  
// 		default:
// 		  libmesh_error();
// 		}
// 	    }

	    
	    
// 	    // Bernstein 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_reg[] = {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_reg[] = {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};
	      	
// 	      	unsigned int i0=i0_reg[i];
// 	      	unsigned int i1=i1_reg[i];
		
// 		if((i== 4||i== 6) && elem->node(0) > elem->node(1)) i0=6-i0;	
// 		if((i== 7||i== 9) && elem->node(1) > elem->node(2)) i1=6-i1;
// 		if((i==10||i==12) && elem->node(3) > elem->node(2)) i0=6-i0;
// 	        if((i==13||i==15) && elem->node(0) > elem->node(3)) i1=6-i1;      
		
	      
// 		switch (j)
// 		  {
// 		    // d()/dxi
// 		  case 0:		  		  
// 		    return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
// 			    FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1,    eta));
		    
// 		  // d()/deta
// 		  case 1:		  		  
// 		    return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0,    xi)*
// 			    FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
		    
// 		default:
// 		  libmesh_error();
// 		}
// 	    }
	    
// 	  default:
// 	    libmesh_error();
// 	  }
//       }
      

      

//       // 5th-order Bernsteins.
//     case FIFTH:
//       {
// 	// Bernstein shape functions on the quadrilateral.
// 	switch (type)
// 	  {
// 	    // Bernstein shape functions on the triangle.
// 	  case TRI6:
// 	    {
// 	      // I have been lazy here and am using finite differences
// 	      // to compute the derivatives!
// 	      const Real eps = 1.e-6;
	      
// 	      libmesh_assert (i < 21);
// 	      libmesh_assert (j < 2);
	      
// 	      unsigned int shape=i;
	      
// 	      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,BERNSTEIN>::shape(elem, order, shape, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, order, shape, 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,BERNSTEIN>::shape(elem, order, shape, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
// 		  }
		  
		  
// 		default:
// 		  libmesh_error();
// 		}
// 	    } // case TRI6
	    
	    
	    
// 	  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_reg[] = {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_reg[] = {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};
		    
// 		    unsigned int i0=i0_reg[i];
// 		    unsigned int i1=i1_reg[i];
		    
// 		    if((i>= 4&&i<= 7) && elem->node(0) > elem->node(1)) i0=7-i0;
// 		    if((i>= 8&&i<=11) && elem->node(1) > elem->node(2)) i1=7-i1;
// 		    if((i>=12&&i<=15) && elem->node(3) > elem->node(2)) i0=7-i0;
// 		    if((i>=16&&i<=19) && elem->node(0) > elem->node(3)) i1=7-i1;      
		    
		    
// 		    switch (j)
// 		      {
// 			// d()/dxi
// 		      case 0:		  		  
// 			return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
// 				FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1,    eta));
			
// 			// d()/deta
// 		      case 1:		  		  
// 			return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0,    xi)*
// 				FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));
			
// 		      default:
// 			libmesh_error();
// 		      }
// 	    } // case QUAD8/9
	    
// 	  default:
// 	    libmesh_error();
// 	  } 
	
//       }
      
      
//       // 6th-order Bernsteins.
//     case SIXTH:
//       {
// 	// Bernstein shape functions on the quadrilateral.
// 	switch (type)
// 	  {
// 	    // Bernstein shape functions on the triangle.
// 	  case TRI6:
// 	    {
// 	      // I have been lazy here and am using finite differences
// 	      // to compute the derivatives!
// 	      const Real eps = 1.e-6;
	      
// 	      libmesh_assert (i < 28);
// 	      libmesh_assert (j < 2);
	      
// 	      unsigned int shape=i;
		      
// 	      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,BERNSTEIN>::shape(elem, order, shape, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, order, shape, 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,BERNSTEIN>::shape(elem, order, shape, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, order, shape, pm))/2./eps;
// 		  }
		  
		  
// 		default:
// 		  libmesh_error();
// 		}
// 	    } // case TRI6
	    
	    
	    
// 	  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_reg[] = {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_reg[] = {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};
	      
// 	      unsigned int i0=i0_reg[i];
// 	      unsigned int i1=i1_reg[i];
	      
// 	      if((i>= 4&&i<= 8) && elem->node(0) > elem->node(1)) i0=8-i0;
// 	      if((i>= 9&&i<=13) && elem->node(1) > elem->node(2)) i1=8-i1;
// 	      if((i>=14&&i<=18) && elem->node(3) > elem->node(2)) i0=8-i0;
// 	      if((i>=19&&i<=23) && elem->node(0) > elem->node(3)) i1=8-i1;      
	      
	      
// 	      switch (j)
// 		{
// 		  // d()/dxi
// 		case 0:		  		  
// 		  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0, 0, xi)*
// 			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1,    eta));
		  
// 		  // d()/deta
// 		case 1:		  		  
// 		  return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0,    xi)*
// 			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1, 0, eta));

// 		default:
// 		  libmesh_error();
// 		}
// 	    } // case QUAD8/9
	    
// 	  default:
// 	    libmesh_error();
// 	  } 

// 	// 7th-order Bernstein.
//       case SEVENTH:
// 	{
// 	  // Szabo-Babuska shape functions on the quadrilateral.
// 	  switch (type)
// 	    {
	      
	      
// 	    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 nodes
// 		  case  7:
// 		  case  9:    			
// 		    if (elem->node(0) > elem->node(1))f = -1.;
// 		    break;
// 		  case 11: // edge 1 nodes
// 		  case 13:	      	
// 		  case 15:
// 		    if (elem->node(1) > elem->node(2))f = -1.;
// 		    break;
// 		  case 17: // edge 2 nodes
// 		  case 19:
// 		  case 21:
// 		    if (elem->node(3) > elem->node(2))f = -1.;
// 		    break;
// 		  case 23: // edge 3 nodes
// 		  case 25:
// 		  case 27:
// 		    if (elem->node(0) > elem->node(3))f = -1.;
// 		  break;
// 		  }	     	       
		
		
// 		switch (j)
// 		  {
// 		    // d()/dxi
// 		  case 0:		  		  
// 		    return f*(FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
// 			      FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[i],    eta));
		    
// 		    // d()/deta
// 		  case 1:		  		  
// 		    return f*(FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[i],    xi)*
// 			      FE<1,BERNSTEIN>::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,BERNSTEIN>::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 Bernstein elements "
            << "are not yet implemented!"
            << std::endl;

  warning_given = true;
  return 0.;
}



template <>
Real FE<2,BERNSTEIN>::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 Bernstein elements "
            << "are not yet implemented!"
            << std::endl;

  warning_given = true;
  return 0.;
}


#endif	// ENABLE_HIGHER_ORDER_SHAPES