fe_bernstein_shape_2d.c

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

// 		  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;
// 		}	     
	      
// 	      return f*(FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
// 			FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta));	      

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

// 	  } // switch type

//       } // case SEVENTH

      
//       // by default throw an error
//     default:
//       std::cerr << "ERROR: Unsupported polynomial order!" << std::endl;
//       libmesh_error();
//     }
  
  libmesh_error();
  return 0.;
}



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



template <>
Real FE<2,BERNSTEIN>::shape_deriv(const Elem* elem,
				  const Order order,
				  const unsigned int i,
				  const unsigned int j,
				  const Point& p)
{
  libmesh_assert (elem != NULL);
  
  const ElemType type = elem->type();
  
  const Order totalorder = static_cast<Order>(order + elem->p_level());

  switch (type)
    {      
      // Hierarchic shape functions on the quadrilateral.
    case QUAD4:
    case QUAD9:
      {
        // Compute quad shape functions as a tensor-product
        const Real xi  = p(0);
        const Real eta = p(1);
      
        libmesh_assert (i < (totalorder+1u)*(totalorder+1u));
	      
        unsigned int i0, i1;

        // Vertex DoFs
        if (i == 0)
          { i0 = 0; i1 = 0; }
        else if (i == 1)
          { i0 = 1; i1 = 0; }
        else if (i == 2)
          { i0 = 1; i1 = 1; }
        else if (i == 3)
          { i0 = 0; i1 = 1; }


        // Edge DoFs
        else if (i < totalorder + 3u)
          { i0 = i - 2; i1 = 0; }
        else if (i < 2u*totalorder + 2)
          { i0 = 1; i1 = i - totalorder - 1; }
        else if (i < 3u*totalorder + 1)
          { i0 = i - 2u*totalorder; i1 = 1; }
        else if (i < 4u*totalorder)
          { i0 = 0; i1 = i - 3u*totalorder + 1; }
        // Interior DoFs
        else
          {
	    unsigned int basisnum = i - 4*totalorder;
	    i0 = square_number_column[basisnum] + 2;
	    i1 = square_number_row[basisnum] + 2;
          }
	
	
	// Flip odd degree of freedom values if necessary
	// to keep continuity on sides
	if     ((i>= 4                 && i<= 4+  totalorder-2u) && elem->point(0) > elem->point(1)) i0=totalorder+2-i0;	//
	else if((i>= 4+  totalorder-1u && i<= 4+2*totalorder-3u) && elem->point(1) > elem->point(2)) i1=totalorder+2-i1;    
	else if((i>= 4+2*totalorder-2u && i<= 4+3*totalorder-4u) && elem->point(3) > elem->point(2)) i0=totalorder+2-i0;
	else if((i>= 4+3*totalorder-3u && i<= 4+4*totalorder-5u) && elem->point(0) > elem->point(3)) i1=totalorder+2-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));
	    
	  }
      }
      
      // Bernstein shape functions on the 8-noded quadrilateral
      // is handled separately.
    case QUAD8:
      {
	libmesh_assert (totalorder < 3);

	const Real xi  = p(0);
	const Real eta = p(1);
	      
	libmesh_assert (i < 8);
	      
	//                                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};
	static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
	switch (j)
	  {
	    // d()/dxi
	  case 0:		      
	    return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
		    FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[i],    eta)
		    +scal[i]*
		    FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
		    FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[8],    eta));
	    
	    // d()/deta
	  case 1:		      
	    return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[i],    xi)*
		    FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
		    +scal[i]*
		    FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[8],    xi)*
		    FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));
	    
	  default:
	    libmesh_error();
	  }	      
      }
      
    case TRI3:
      libmesh_assert (totalorder<2);
    case TRI6:
      {
	// I have been lazy here and am using finite differences
	// to compute the derivatives!
	const Real eps = 1.e-6;
		
	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, totalorder, i, pp) -
		      FE<2,BERNSTEIN>::shape(elem, totalorder, 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,BERNSTEIN>::shape(elem, totalorder, i, pp) -
		      FE<2,BERNSTEIN>::shape(elem, totalorder, i, pm))/2./eps;
	    }
	    
	    
	  default:
	    libmesh_error();
	  }
      }
      
    default:
      {
	std::cerr << "ERROR: Unsupported element type!" << std::endl;
	libmesh_error();
      }
      
    }
  
  // old code
//   switch (totalorder)
//     {            
      
//       // 1st order Bernsteins are aquivalent to Lagrange elements.
//     case FIRST:
//       {
// 	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 < 3);
// 	      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,BERNSTEIN>::shape(elem, order, i, pp) -
// 			    FE<2,BERNSTEIN>::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,BERNSTEIN>::shape(elem, order, i, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, order, i, 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 < 4);
	      
// 	      //                                0  1  2  3
// 	      static const unsigned int i0[] = {0, 1, 1, 0};
// 	      static const unsigned int i1[] = {0, 0, 1, 1};
	      
// 	      switch (j)
// 		{
// 		  // d()/dxi
// 		case 0:		      
// 		  return (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 (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();
// 	  }
//       }

//       // second order Bernsteins.
//     case SECOND:
//       {
// 	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 < 6);
// 	      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,BERNSTEIN>::shape(elem, order, i, pp) -
// 			    FE<2,BERNSTEIN>::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,BERNSTEIN>::shape(elem, order, i, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
// 		  }
		  

// 		default:
// 		  libmesh_error();
// 		}
// 	    }
	    
// 	    // Bernstein shape functions on the 8-noded quadrilateral.
// 	  case QUAD8:
// 	    {
// 	      const Real xi  = p(0);
// 	      const Real eta = p(1);
	      
// 	      libmesh_assert (i < 8);
	      
// 	      //                                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};
// 	      static const Real scal[] = {-0.25, -0.25, -0.25, -0.25, 0.5, 0.5, 0.5, 0.5};
// 	      switch (j)
// 		{
// 		  // d()/dxi
// 		case 0:		      
// 		  return (FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[i], 0, xi)*
// 			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[i],    eta)
// 			  +scal[i]*
// 			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[8], 0, xi)*
// 			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[8],    eta));
		  
// 		  // d()/deta
// 		case 1:		      
// 		  return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[i],    xi)*
// 			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[i], 0, eta)
// 			  +scal[i]*
// 			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[8],    xi)*
// 			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[8], 0, eta));

// 		default:
// 		  libmesh_error();
// 		}	      
// 	    }


// 	    // Bernstein shape functions on the 9-noded quadrilateral.
// 	  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};
	      
// 	      switch (j)
// 		{
// 		  // d()/dxi
// 		case 0:		      
// 		  return (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 (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();
// 	  }
//       }
      

//       // 3rd-order Bernsteins.
//     case THIRD:
//       {
// 	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 < 10);
// 	      libmesh_assert (j < 2);
	      
	      
// 	      unsigned int shape=i;

// 	      /**
// 	       * Note, since the derivatives are computed using FE::shape
// 	       * the continuity along neighboring elements (shape_flip)
// 	       * is already considered. 
// 	       */

	      
// 	      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, totalorder, shape, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, totalorder, 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, totalorder, shape, pp) -
// 			    FE<2,BERNSTEIN>::shape(elem, totalorder, shape, pm))/2./eps;
// 		  }
		  
		  
// 		default: