fe_bernstein_shape_3d.c

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

	      }
	      
	      
	      return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi_mapped)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta_mapped)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta_mapped));
	    }

	    
	  default:
	    libmesh_error();
	  }	
      }
      
      
    default:
      libmesh_error();
    }
  
#endif
  
  libmesh_error();
  return 0.;
}




template <>
Real FE<3,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 and face orientation is needed."
	    << std::endl;
  libmesh_error();
  
  return 0.;
}



template <>
Real FE<3,BERNSTEIN>::shape_deriv(const Elem* elem,
				  const Order order,
				  const unsigned int i,
				  const unsigned int j,
				  const Point& p)
{

#if DIM == 3
  libmesh_assert (elem != NULL);
  const ElemType type = elem->type();

  const Order totalorder = static_cast<Order>(order + elem->p_level());
  
  libmesh_assert (j < 3);
  
  switch (totalorder)
    {
      
      
      // 1st order Bernstein.
    case FIRST:
      {
	switch (type)
	  {
	    // Bernstein shape functions on the tetrahedron.
	  case TET4:
	  case TET10:
	    {
	      // I have been lazy here and am using finite differences
	      // to compute the derivatives!
	      const Real eps = 1.e-6;
	      
	      libmesh_assert (i < 4);
	      libmesh_assert (j < 3);
	      
	      
	      switch (j)
		{
		  //  d()/dxi
		case 0:
		  {
		    const Point pp(p(0)+eps, p(1), p(2));
		    const Point pm(p(0)-eps, p(1), p(2));
		    
		    return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
			    FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
		  }

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

		    return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
			    FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
		  }
		  // d()/dzeta
		case 2:
		  {
		    const Point pp(p(0), p(1), p(2)+eps);
		    const Point pm(p(0), p(1), p(2)-eps);

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




	    // Bernstein shape functions on the hexahedral.
	  case HEX20:
	  case HEX27:
	    {
	      libmesh_assert (i<8);
	
	      // Compute hex shape functions as a tensor-product
	      const Real xi   = p(0);
	      const Real eta  = p(1);
	      const Real zeta = p(2);
	
	      // The only way to make any sense of this
	      // is to look at the mgflo/mg2/mgf documentation
	      // and make the cut-out cube!
	      //                                0  1  2  3  4  5  6  7
	      static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
	      static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 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)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[i],    zeta));

		  // 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)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[i],    zeta));

		  // d()/dzeta
		case 2:
		  return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[i],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[i],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));

		default:
		  libmesh_error();
		}
	    }


	  default:
	    libmesh_error();
	  }
      }




    case SECOND:
      {
	switch (type)
	  {
        	// Bernstein shape functions on the tetrahedron.
	  case TET10:
		{
	      // 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 < 3);


		switch (j)
		{
		  //  d()/dxi
		case 0:
		  {
		    const Point pp(p(0)+eps, p(1), p(2));
		    const Point pm(p(0)-eps, p(1), p(2));

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

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

		    return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
			    FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
		  }
		  // d()/dzeta
		case 2:
		  {
		    const Point pp(p(0), p(1), p(2)+eps);
		    const Point pm(p(0), p(1), p(2)-eps);

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

		// Bernstein shape functions on the hexahedral.
	  case HEX20:
	    {
	      libmesh_assert (i<20);
	
	      // Compute hex shape functions as a tensor-product
	      const Real xi   = p(0);
	      const Real eta  = p(1);
	      const Real zeta = p(2);
	
	      // The only way to make any sense of this
	      // is to look at the mgflo/mg2/mgf documentation
	      // and make the cut-out cube!
	      //                                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
	      static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
	      static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
	      static const Real scal20[] =     {-0.25, -0.25, -0.25, -0.25, 0,     0,     0,     0,     0.5,   0.5,   0.5,   0.5,   0,     0,     0,     0,     0,     0,     0,     0};
	      static const Real scal21[] =     {-0.25, -0.25, 0,     0,     -0.25, -0.25, 0,     0,     0.5,   0,     0,     0,     0.5,   0.5,   0,     0,     0.5,   0,     0,     0};
	      static const Real scal22[] =     {0,     -0.25, -0.25, 0,     0,     -0.25, -0.25, 0,     0,     0.5,   0,     0,     0,     0.5,   0.5,   0,     0,     0.5,   0,     0};
	      static const Real scal23[] =     {0,     0,     -0.25, -0.25, 0,     0,     -0.25, -0.25, 0,     0,     0.5,   0,     0,     0,     0.5,   0.5,   0,     0,     0.5,   0};
	      static const Real scal24[] =     {-0.25, 0,     0,     -0.25, -0.25, 0,     0,     -0.25, 0,     0,     0,     0.5,   0.5,   0,     0,     0.5,   0,     0,     0,     0.5};
	      static const Real scal25[] =     {0,     0,     0,     0,     -0.25, -0.25, -0.25, -0.25, 0,     0,     0,     0,     0,     0,     0,     0,     0.5,   0.5,   0.5,   0.5};
	      static const Real scal26[] =     {-0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, 0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25};

	      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)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[i],    zeta)
			  +scal20[i]*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[20], 0, xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[20],    eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[20],    zeta)
			  +scal21[i]*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[21], 0, xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[21],    eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[21],    zeta)
			  +scal22[i]*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[22], 0, xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[22],    eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[22],    zeta)
			  +scal23[i]*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[23], 0, xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[23],    eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[23],    zeta)
			  +scal24[i]*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[24], 0, xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[24],    eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[24],    zeta)
			  +scal25[i]*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[25], 0, xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[25],    eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[25],    zeta)
			  +scal26[i]*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i0[26], 0, xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[26],    eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[26],    zeta));

		  // 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)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[i],    zeta)
			  +scal20[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[20],     xi)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[20], 0, eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[20],    zeta)
			  +scal21[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[21],     xi)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[21], 0, eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[21],    zeta)
			  +scal22[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[22],     xi)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[22], 0, eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[22],    zeta)
			  +scal23[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[23],     xi)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[23], 0, eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[23],    zeta)
			  +scal24[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[24],     xi)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[24], 0, eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[24],    zeta)
			  +scal25[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[25],     xi)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[25], 0, eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[25],    zeta)
			  +scal26[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[26],     xi)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i1[26], 0, eta)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[26],    zeta));

		  // d()/dzeta
		case 2:
		  return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[i],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[i],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta)
			  +scal20[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[20],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[20],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[20], 0, zeta)
			  +scal21[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[21],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[21],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[21], 0, zeta)
			  +scal22[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[22],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[22],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[22], 0, zeta)
			  +scal23[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[23],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[23],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[23], 0, zeta)
			  +scal24[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[24],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[24],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[24], 0, zeta)
			  +scal25[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[25],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[25],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[25], 0, zeta)
			  +scal26[i]*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[26],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[26],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[26], 0, zeta));

		default:
		  libmesh_error();
		}
	    }

	    // Bernstein shape functions on the hexahedral.
	  case HEX27:
	    {
	      libmesh_assert (i<27);
	
	      // Compute hex shape functions as a tensor-product
	      const Real xi   = p(0);
	      const Real eta  = p(1);
	      const Real zeta = p(2);
	
	      // The only way to make any sense of this
	      // is to look at the mgflo/mg2/mgf documentation
	      // and make the cut-out cube!
	      //                                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
	      static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
	      static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 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)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[i],    zeta));

		  // 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)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i2[i],    zeta));

		  // d()/dzeta
		case 2:
		  return (FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i0[i],    xi)*
			  FE<1,BERNSTEIN>::shape      (EDGE3, totalorder, i1[i],    eta)*
			  FE<1,BERNSTEIN>::shape_deriv(EDGE3, totalorder, i2[i], 0, zeta));

		default:
		  libmesh_error();
		}
	    }


	  default:
	    libmesh_error();
	  }
      }

      

      // 3rd-order Bernstein.
    case THIRD:
      {
	switch (type)
	  {

// 	    // Bernstein shape functions derivatives.
// 	  case TET10:
// 	    {
// 	      // I have been lazy here and am using finite differences
// 	      // to compute the derivatives!
// 	      const Real eps = 1.e-6;
	      
// 	      libmesh_assert (i < 20);
// 	      libmesh_assert (j < 3);    
 
// 		switch (j)
// 		{
// 		  //  d()/dxi
// 		case 0:
// 		  {
// 		    const Point pp(p(0)+eps, p(1), p(2));
// 		    const Point pm(p(0)-eps, p(1), p(2));

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

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

// 		    return (FE<3,BERNSTEIN>::shape(elem, order, i, pp) -
// 			    FE<3,BERNSTEIN>::shape(elem, order, i, pm))/2./eps;
// 		  }
// 		  // d()/dzeta
// 		case 2:
// 		  {
// 		    const Point pp(p(0), p(1), p(2)+eps);
// 		    const Point pm(p(0), p(1), p(2)-eps);