fe_lagrange_shape_3d.c

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

	    }
	    
	    // quadratic tetrahedral shape functions	    
	  case TET10:
	    {
	      libmesh_assert (i<10);
	
	      // Area coordinates, pg. 205, Vol. I, Carey, Oden, Becker FEM
	      const Real zeta1 = p(0);
	      const Real zeta2 = p(1);
	      const Real zeta3 = p(2);
	      const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
	
	      const Real dzeta0dxi = -1.;
	      const Real dzeta1dxi =  1.;
	      const Real dzeta2dxi =  0.;
	      const Real dzeta3dxi =  0.;
	
	      const Real dzeta0deta = -1.;
	      const Real dzeta1deta =  0.;
	      const Real dzeta2deta =  1.;
	      const Real dzeta3deta =  0.;
	
	      const Real dzeta0dzeta = -1.;
	      const Real dzeta1dzeta =  0.;
	      const Real dzeta2dzeta =  0.;
	      const Real dzeta3dzeta =  1.;

	      switch (j)
		{
		  // d()/dxi
		case 0:
		  {
		    switch(i)
		      {
		      case 0:
			return (4.*zeta0 - 1.)*dzeta0dxi;
		 
		      case 1:
			return (4.*zeta1 - 1.)*dzeta1dxi;
		 
		      case 2:
			return (4.*zeta2 - 1.)*dzeta2dxi;
		 
		      case 3:
			return (4.*zeta3 - 1.)*dzeta3dxi;
		 
		      case 4:
			return 4.*(zeta0*dzeta1dxi + dzeta0dxi*zeta1);
		 
		      case 5:
			return 4.*(zeta1*dzeta2dxi + dzeta1dxi*zeta2);
		 
		      case 6:
			return 4.*(zeta0*dzeta2dxi + dzeta0dxi*zeta2);
		 
		      case 7:
			return 4.*(zeta0*dzeta3dxi + dzeta0dxi*zeta3);
		 
		      case 8:
			return 4.*(zeta1*dzeta3dxi + dzeta1dxi*zeta3);
		 
		      case 9:
			return 4.*(zeta2*dzeta3dxi + dzeta2dxi*zeta3);
		 
		      default:
			libmesh_error();
		      }
		  }
		      
		  // d()/deta
		case 1:
		  {
		    switch(i)
		      {
		      case 0:
			return (4.*zeta0 - 1.)*dzeta0deta;
		 
		      case 1:
			return (4.*zeta1 - 1.)*dzeta1deta;
		 
		      case 2:
			return (4.*zeta2 - 1.)*dzeta2deta;
		 
		      case 3:
			return (4.*zeta3 - 1.)*dzeta3deta;
		 
		      case 4:
			return 4.*(zeta0*dzeta1deta + dzeta0deta*zeta1);
		 
		      case 5:
			return 4.*(zeta1*dzeta2deta + dzeta1deta*zeta2);
		 
		      case 6:
			return 4.*(zeta0*dzeta2deta + dzeta0deta*zeta2);
		 
		      case 7:
			return 4.*(zeta0*dzeta3deta + dzeta0deta*zeta3);
		 
		      case 8:
			return 4.*(zeta1*dzeta3deta + dzeta1deta*zeta3);
		 
		      case 9:
			return 4.*(zeta2*dzeta3deta + dzeta2deta*zeta3);
		 
		      default:
			libmesh_error();
		      }
		  }
		      
		  // d()/dzeta
		case 2:
		  {
		    switch(i)
		      {
		      case 0:
			return (4.*zeta0 - 1.)*dzeta0dzeta;
		 
		      case 1:
			return (4.*zeta1 - 1.)*dzeta1dzeta;
		 
		      case 2:
			return (4.*zeta2 - 1.)*dzeta2dzeta;
		 
		      case 3:
			return (4.*zeta3 - 1.)*dzeta3dzeta;
		 
		      case 4:
			return 4.*(zeta0*dzeta1dzeta + dzeta0dzeta*zeta1);
		 
		      case 5:
			return 4.*(zeta1*dzeta2dzeta + dzeta1dzeta*zeta2);
		 
		      case 6:
			return 4.*(zeta0*dzeta2dzeta + dzeta0dzeta*zeta2);
		 
		      case 7:
			return 4.*(zeta0*dzeta3dzeta + dzeta0dzeta*zeta3);
		 
		      case 8:
			return 4.*(zeta1*dzeta3dzeta + dzeta1dzeta*zeta3);
		 
		      case 9:
			return 4.*(zeta2*dzeta3dzeta + dzeta2dzeta*zeta3);
		 
		      default:
			libmesh_error();
		      }
		  }

		default:
		  libmesh_error();
		}
	    }


	    
	    // quadradic prism shape functions	    
	  case PRISM18:
	    {
	      libmesh_assert (i<18);
	
	      // 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  6  7  8  9 10 11 12 13 14 15 16 17 
	      static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
	      static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};

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

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

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

#endif
  
  libmesh_error();
  return 0.;  
}



template <>
Real FE<3,LAGRANGE>::shape_deriv(const Elem* elem,
				 const Order order,
				 const unsigned int i,
				 const unsigned int j,
				 const Point& p)
{
  libmesh_assert (elem != NULL);
      
  // call the orientation-independent shape function derivatives
  return FE<3,LAGRANGE>::shape_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
}



template <>
Real FE<3,LAGRANGE>::shape_second_deriv(const ElemType type,
				        const Order order,
				        const unsigned int i,
				        const unsigned int j,
				        const Point& p)
{
#if DIM == 3
  
  libmesh_assert (j<6);
  
  switch (order)
    {
      // linear Lagrange shape functions
    case FIRST:
      {
	return 0.;
      }
      
      // quadratic Lagrange shape functions
    case SECOND:
      {
	switch (type)
	  {

	    // serendipity hexahedral quadratic shape functions
	  case HEX20:
	    {
              static bool warning_given_HEX20 = false;

              if (!warning_given_HEX20)
              std::cerr << "Second derivatives for 2D Lagrangian HEX20"
                        << " elements are not yet implemented!"
                        << std::endl;
              warning_given_HEX20 = true;
	    }

	    // 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)
		{
                // d^2()/dxi^2
		case 0:
		  return (FE<1,LAGRANGE>::shape_second_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));

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

                // d^2()/deta^2
		case 2:     
		  return (FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i0[i], xi)*
			  FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i1[i], 0, eta)*
			  FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i2[i], zeta));

                // d^2()/dxidzeta
		case 3:     
		  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_deriv(EDGE3, SECOND, i2[i], 0, zeta));

                // d^2()/detadzeta
		case 4:     
		  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_deriv(EDGE3, SECOND, i2[i], 0, zeta));

                // d^2()/dzeta^2
		case 5:     
		  return (FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i0[i], xi)*
			  FE<1,LAGRANGE>::shape      (EDGE3, SECOND, i1[i], eta)*
			  FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i2[i], 0, zeta));
		  
		default:
		  {
		    libmesh_error();
		  }
		}
	    }
	    
	    // quadratic tetrahedral shape functions	    
	  case TET10:
	    {
	      // The area coordinates are the same as used for the
	      // shape() and shape_deriv() functions.
	      // const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
	      // const Real zeta1 = p(0);
	      // const Real zeta2 = p(1);
	      // const Real zeta3 = p(2);
	      static const Real dzetadxi[4][3] =
		{
		  {-1., -1., -1.},
		  {1.,   0.,  0.},
		  {0.,   1.,  0.},
		  {0.,   0.,  1.}
		};

	      // Convert from j -> (j,k) indices for independent variable
	      // (0=xi, 1=eta, 2=zeta)
	      static const unsigned short int independent_var_indices[6][2] =
		{
		  {0, 0}, // d^2 phi / dxi^2
		  {0, 1}, // d^2 phi / dxi deta
		  {1, 1}, // d^2 phi / deta^2
		  {0, 2}, // d^2 phi / dxi dzeta
		  {1, 2}, // d^2 phi / deta dzeta
		  {2, 2}  // d^2 phi / dzeta^2
		};

	      // Convert from i -> zeta indices.  Each quadratic shape
	      // function for the Tet10 depends on up to two of the zeta
	      // area coordinate functions (see the shape() function above).
	      // This table just tells which two area coords it uses.
	      static const unsigned short int zeta_indices[10][2] =
		{
		  {0, 0}, 
		  {1, 1}, 
		  {2, 2}, 
		  {3, 3},
		  {0, 1}, 		  
		  {1, 2}, 
		  {2, 0},
		  {0, 3},
		  {1, 3},
		  {2, 3},  		  		  
		};

	      // Look up the independent variable indices for this value of j.
	      const unsigned int my_j = independent_var_indices[j][0];
	      const unsigned int my_k = independent_var_indices[j][1];
	      
	      if (i<4)
		{
		  return 4.*dzetadxi[i][my_j]*dzetadxi[i][my_k];
		}
	      
	      else if (i<10)
		{
		  const unsigned short int my_m = zeta_indices[i][0];
		  const unsigned short int my_n = zeta_indices[i][1];

		  return 4.*(dzetadxi[my_n][my_j]*dzetadxi[my_m][my_k] +
			     dzetadxi[my_m][my_j]*dzetadxi[my_n][my_k] );
		}
	      else
		{
		  std::cerr<< "Invalid shape function index " << i << std::endl;
		  libmesh_error();
		}
	    }

	    
	    // quadradic prism shape functions	    
	  case PRISM18:
	    {
	      libmesh_assert (i<18);
	
	      // 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  6  7  8  9 10 11 12 13 14 15 16 17 
	      static const unsigned int i0[] = {0, 0, 0, 1, 1, 1, 0, 0, 0, 2, 2, 2, 1, 1, 1, 2, 2, 2};
	      static const unsigned int i1[] = {0, 1, 2, 0, 1, 2, 3, 4, 5, 0, 1, 2, 3, 4, 5, 3, 4, 5};

	      switch (j)
		{
		  // d^2()/dxi^2
		case 0:
		  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 0, p2d)*
			  FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));

		  // d^2()/dxideta
		case 1:
		  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 1, p2d)*
			  FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));
		  
		  // d^2()/deta^2
		case 2:
		  return (FE<2,LAGRANGE>::shape_second_deriv(TRI6, SECOND, i1[i], 2, p2d)*
			  FE<1,LAGRANGE>::shape(EDGE3, SECOND, i0[i], p1d));

		  // d^2()/dxidzeta
		case 3:
		  return (FE<2,LAGRANGE>::shape_deriv(TRI6,  SECOND, i1[i], 0, p2d)*
			  FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));

		  // d^2()/detadzeta
		case 4:
		  return (FE<2,LAGRANGE>::shape_deriv(TRI6,  SECOND, i1[i], 1, p2d)*
			  FE<1,LAGRANGE>::shape_deriv(EDGE3, SECOND, i0[i], 0, p1d));

		  // d^2()/dzeta^2
		case 5:
		  return (FE<2,LAGRANGE>::shape(TRI6,  SECOND, i1[i], p2d)*
			  FE<1,LAGRANGE>::shape_second_deriv(EDGE3, SECOND, i0[i], 0, p1d));
		}
	    }

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

#endif
  
  libmesh_error();
  return 0.;  
}



template <>
Real FE<3,LAGRANGE>::shape_second_deriv(const Elem* elem,
				 const Order order,
				 const unsigned int i,
				 const unsigned int j,
				 const Point& p)
{
  libmesh_assert (elem != NULL);
      
  // call the orientation-independent shape function derivatives
  return FE<3,LAGRANGE>::shape_second_deriv(elem->type(), static_cast<Order>(order + elem->p_level()), i, j, p);
}