fe_lagrange_shape_3d.c

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

		      case 2:
			return dzeta2deta;
		 
		      case 3:
			return dzeta3deta;
		 
		      default:
			libmesh_error();
		      }
		  }
	    
		  // d()/dzeta
		case 2:
		  {
		    switch(i)
		      {
		      case 0:
			return dzeta0dzeta;
		 
		      case 1:
			return dzeta1dzeta;
		 
		      case 2:
			return dzeta2dzeta;
		 
		      case 3:
			return dzeta3dzeta;
		 
		      default:
			libmesh_error();
		      }
		  }
		}
	    }
	    
	    // linear prism shape functions	    
	  case PRISM6:
	  case PRISM15:
	  case PRISM18:
	    {
	      libmesh_assert (i<6);
	
	      // Compute prism shape functions as a tensor-product
	      // of a triangle and an edge
	
	      Point p2d(p(0),p(1));
	      Point p1d(p(2));
	
	      //                                0  1  2  3  4  5  
	      static const unsigned int i0[] = {0, 0, 0, 1, 1, 1};
	      static const unsigned int i1[] = {0, 1, 2, 0, 1, 2};

	      switch (j)
		{
		  // d()/dxi
		case 0:
		  return (FE<2,LAGRANGE>::shape_deriv(TRI3,  FIRST, i1[i], 0, p2d)*
			  FE<1,LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));
		  
		  // d()/deta
		case 1:
		  return (FE<2,LAGRANGE>::shape_deriv(TRI3,  FIRST, i1[i], 1, p2d)*
			  FE<1,LAGRANGE>::shape(EDGE2, FIRST, i0[i], p1d));

		  // d()/dzeta
		case 2:
		  return (FE<2,LAGRANGE>::shape(TRI3,  FIRST, i1[i], p2d)*
			  FE<1,LAGRANGE>::shape_deriv(EDGE2, FIRST, i0[i], 0, p1d));
		}
	    }

	    // linear pyramid shape functions
	  case PYRAMID5:
	    {
	      libmesh_assert (i<5);
	      	
	      const Real xi   = p(0);
	      const Real eta  = p(1);
	      const Real zeta = p(2);
	      const Real eps  = 1.e-35;

	      switch (j)
		{
		  // d/dxi
		case 0:	   
		  switch(i)
		    {
		    case 0:
		      return  .25*(zeta + eta - 1.)/((1. - zeta) + eps);
		
		    case 1:
		      return -.25*(zeta + eta - 1.)/((1. - zeta) + eps);
		
		    case 2:
		      return -.25*(zeta - eta - 1.)/((1. - zeta) + eps);
		
		    case 3:
		      return  .25*(zeta - eta - 1.)/((1. - zeta) + eps);
		
		    case 4:
		      return 0;
		
		    default:
		      libmesh_error();
		    }


		  // d/deta
		case 1:	   
		  switch(i)
		    {
		    case 0:
		      return  .25*(zeta + xi - 1.)/((1. - zeta) + eps);
		
		    case 1:
		      return  .25*(zeta - xi - 1.)/((1. - zeta) + eps);
		
		    case 2:
		      return -.25*(zeta - xi - 1.)/((1. - zeta) + eps);
		
		    case 3:
		      return -.25*(zeta + xi - 1.)/((1. - zeta) + eps);
		
		    case 4:
		      return 0;
	    
		    default:
		      libmesh_error();
		    }


		  // d/dzeta
		case 2:	    
		  switch(i)
		    {
		    case 0:
		      {
			const Real a=1.;
			const Real b=1.;
		  
			return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
				     (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
				    ((1. - zeta)*(1. - zeta) + eps));
		      }
		
		    case 1:
		      {
			const Real a=-1.;
			const Real b=1.;
		  
			return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
				     (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
				    ((1. - zeta)*(1. - zeta) + eps));
		      }
		
		    case 2:
		      {
			const Real a=-1.;
			const Real b=-1.;
		  
			return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
				     (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
				    ((1. - zeta)*(1. - zeta) + eps));
		      }
		
		    case 3:
		      {
			const Real a=1.;
			const Real b=-1.;
		  
			return .25*(((zeta + a*xi - 1.)*(zeta + b*eta - 1.) +
				     (1. - zeta)*((zeta + a*xi -1.) + (zeta + b*eta - 1.)))/
				    ((1. - zeta)*(1. - zeta) + eps));
		      }
		
		    case 4:
		      return 1.;
		
		    default:
		      libmesh_error();
		    }
	    

		default:
		  libmesh_error();
		}
	    }
	    
	    
	  default:
	    {
	      std::cerr << "ERROR: Unsupported 3D element type!: " << type
			<< std::endl;
	      libmesh_error();
	    }
	  }
      }

      
      // quadratic Lagrange shape functions
    case SECOND:
      {
	switch (type)
	  {

	    // serendipity hexahedral quadratic shape functions
	  case HEX20:
	    {
	      libmesh_assert (i<20);
	
	      const Real xi   = p(0);
	      const Real eta  = p(1);
	      const Real zeta = p(2);

	      // these functions are defined for (x,y,z) in [0,1]^3
	      // so transform the locations
	      const Real x = .5*(xi   + 1.);
	      const Real y = .5*(eta  + 1.);
	      const Real z = .5*(zeta + 1.);

	      // and don't forget the chain rule!
	
	      switch (j)
		{

		  // d/dx*dx/dxi
		case 0:
		  switch (i)
		    {
		    case 0:
		      return .5*(1. - y)*(1. - z)*((1. - x)*(-2.) +
						   (-1.)*(1. - 2.*x - 2.*y - 2.*z));
		      libmesh_error(); // done to here
		
		    case 1:
		      return .5*(1. - y)*(1. - z)*(x*(2.) +
						   (1.)*(2.*x - 2.*y - 2.*z - 1.));
		
		    case 2:
		      return .5*y*(1. - z)*(x*(2.) +
					    (1.)*(2.*x + 2.*y - 2.*z - 3.));
		
		    case 3:
		      return .5*y*(1. - z)*((1. - x)*(-2.) +
					    (-1.)*(2.*y - 2.*x - 2.*z - 1.));
		
		    case 4:
		      return .5*(1. - y)*z*((1. - x)*(-2.) +
					    (-1.)*(2.*z - 2.*x - 2.*y - 1.));
		
		    case 5:
		      return .5*(1. - y)*z*(x*(2.) +
					    (1.)*(2.*x - 2.*y + 2.*z - 3.));
		
		    case 6:
		      return .5*y*z*(x*(2.) +
				     (1.)*(2.*x + 2.*y + 2.*z - 5.));
		
		    case 7:
		      return .5*y*z*((1. - x)*(-2.) +
				     (-1.)*(2.*y - 2.*x + 2.*z - 3.));
		
		    case 8:
		      return 2.*(1. - y)*(1. - z)*(1. - 2.*x);
		
		    case 9:
		      return 2.*y*(1. - y)*(1. - z);
		
		    case 10:
		      return 2.*y*(1. - z)*(1. - 2.*x);
		
		    case 11:
		      return 2.*y*(1. - y)*(1. - z)*(-1.);
		
		    case 12:
		      return 2.*(1. - y)*z*(1. - z)*(-1.);
		
		    case 13:
		      return 2.*(1. - y)*z*(1. - z);
		
		    case 14:
		      return 2.*y*z*(1. - z);
		
		    case 15:
		      return 2.*y*z*(1. - z)*(-1.);
		
		    case 16:
		      return 2.*(1. - y)*z*(1. - 2.*x);
		
		    case 17:
		      return 2.*y*(1. - y)*z;
		
		    case 18:
		      return 2.*y*z*(1. - 2.*x);
		
		    case 19:
		      return 2.*y*(1. - y)*z*(-1.);
		
		    default:
		      libmesh_error();
		    }


		  // d/dy*dy/deta
		case 1:
		  switch (i)
		    {
		    case 0:
		      return .5*(1. - x)*(1. - z)*((1. - y)*(-2.) +
						   (-1.)*(1. - 2.*x - 2.*y - 2.*z));
		
		    case 1:
		      return .5*x*(1. - z)*((1. - y)*(-2.) +
					    (-1.)*(2.*x - 2.*y - 2.*z - 1.));
		
		    case 2:
		      return .5*x*(1. - z)*(y*(2.) +
					    (1.)*(2.*x + 2.*y - 2.*z - 3.));
		
		    case 3:
		      return .5*(1. - x)*(1. - z)*(y*(2.) +
						   (1.)*(2.*y - 2.*x - 2.*z - 1.));
		
		    case 4:
		      return .5*(1. - x)*z*((1. - y)*(-2.) +
					    (-1.)*(2.*z - 2.*x - 2.*y - 1.));
		
		    case 5:
		      return .5*x*z*((1. - y)*(-2.) +
				     (-1.)*(2.*x - 2.*y + 2.*z - 3.));
		
		    case 6:
		      return .5*x*z*(y*(2.) +
				     (1.)*(2.*x + 2.*y + 2.*z - 5.));
		
		    case 7:
		      return .5*(1. - x)*z*(y*(2.) +
					    (1.)*(2.*y - 2.*x + 2.*z - 3.));
		
		    case 8:
		      return 2.*x*(1. - x)*(1. - z)*(-1.);
		
		    case 9:
		      return 2.*x*(1. - z)*(1. - 2.*y);
		
		    case 10:
		      return 2.*x*(1. - x)*(1. - z);
		
		    case 11:
		      return 2.*(1. - x)*(1. - z)*(1. - 2.*y);
		
		    case 12:
		      return 2.*(1. - x)*z*(1. - z)*(-1.);
		
		    case 13:
		      return 2.*x*z*(1. - z)*(-1.);
		
		    case 14:
		      return 2.*x*z*(1. - z);
		
		    case 15:
		      return 2.*(1. - x)*z*(1. - z);
		
		    case 16:
		      return 2.*x*(1. - x)*z*(-1.);
		
		    case 17:
		      return 2.*x*z*(1. - 2.*y);
		
		    case 18:
		      return 2.*x*(1. - x)*z;
		
		    case 19:
		      return 2.*(1. - x)*z*(1. - 2.*y);
		
		    default:
		      libmesh_error();
		    }


		  // d/dz*dz/dzeta
		case 2:
		  switch (i)
		    {
		    case 0:
		      return .5*(1. - x)*(1. - y)*((1. - z)*(-2.) +
						   (-1.)*(1. - 2.*x - 2.*y - 2.*z));
		
		    case 1:
		      return .5*x*(1. - y)*((1. - z)*(-2.) +
					    (-1.)*(2.*x - 2.*y - 2.*z - 1.));
		
		    case 2:
		      return .5*x*y*((1. - z)*(-2.) +
				     (-1.)*(2.*x + 2.*y - 2.*z - 3.));
		
		    case 3:
		      return .5*(1. - x)*y*((1. - z)*(-2.) +
					    (-1.)*(2.*y - 2.*x - 2.*z - 1.));
		
		    case 4:
		      return .5*(1. - x)*(1. - y)*(z*(2.) +
						   (1.)*(2.*z - 2.*x - 2.*y - 1.));
		
		    case 5:
		      return .5*x*(1. - y)*(z*(2.) +
					    (1.)*(2.*x - 2.*y + 2.*z - 3.));
		
		    case 6:
		      return .5*x*y*(z*(2.) +
				     (1.)*(2.*x + 2.*y + 2.*z - 5.));
		
		    case 7:
		      return .5*(1. - x)*y*(z*(2.) +
					    (1.)*(2.*y - 2.*x + 2.*z - 3.));
		
		    case 8:
		      return 2.*x*(1. - x)*(1. - y)*(-1.);
		
		    case 9:
		      return 2.*x*y*(1. - y)*(-1.);
		
		    case 10:
		      return 2.*x*(1. - x)*y*(-1.);
		
		    case 11:
		      return 2.*(1. - x)*y*(1. - y)*(-1.);
		
		    case 12:
		      return 2.*(1. - x)*(1. - y)*(1. - 2.*z);
		
		    case 13:
		      return 2.*x*(1. - y)*(1. - 2.*z);
		
		    case 14:
		      return 2.*x*y*(1. - 2.*z);
		
		    case 15:
		      return 2.*(1. - x)*y*(1. - 2.*z);
		
		    case 16:
		      return 2.*x*(1. - x)*(1. - y);
		
		    case 17:
		      return 2.*x*y*(1. - y);
		
		    case 18:
		      return 2.*x*(1. - x)*y;
		
		    case 19:
		      return 2.*(1. - x)*y*(1. - y);
		
		    default:
		      libmesh_error();
		    }
		}
	    }

	    // triquadraic hexahedral shape funcions	    
	  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)
		{
		case 0:
		  return (FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, xi)*
			  FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i1[i], eta)*
			  FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i2[i], zeta));

		case 1:     
		  return (FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i0[i], xi)*
			  FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i1[i], 0, eta)*
			  FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i2[i], zeta));

		case 2:     
		  return (FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i0[i], xi)*
			  FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i1[i], eta)*
			  FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i2[i], 0, zeta));
		  
		default:
		  {
		    libmesh_error();
		  }
		}