fe_hermite_shape_3d.c

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

      switch ((i - 64 - 12*4*e) / (2*e*e))
        {
        case 0:
          bases1D[0] = square_number_column[basisnum / 2];
          bases1D[1] = square_number_row[basisnum / 2];
          bases1D[2] = basisnum % 2 * 2;
          if (basisnum % 2)
            coef *= dxdxi[2][0];
          break;
        case 1:
          bases1D[0] = square_number_column[basisnum / 2];
          bases1D[1] = basisnum % 2 * 2;
          bases1D[2] = square_number_row[basisnum / 2];
          if (basisnum % 2)
            coef *= dxdxi[1][0];
          break;
        case 2:
          bases1D[0] = basisnum % 2 * 2 + 1;
          bases1D[1] = square_number_column[basisnum / 2];
          bases1D[2] = square_number_row[basisnum / 2];
          if (basisnum % 2)
            coef *= dxdxi[0][1];
          break;
        case 3:
          bases1D[0] = square_number_column[basisnum / 2];
          bases1D[1] = basisnum % 2 * 2 + 1;
          bases1D[2] = square_number_row[basisnum / 2];
          if (basisnum % 2)
            coef *= dxdxi[1][1];
          break;
        case 4:
          bases1D[0] = basisnum % 2 * 2;
          bases1D[1] = square_number_column[basisnum / 2];
          bases1D[2] = square_number_row[basisnum / 2];
          if (basisnum % 2)
            coef *= dxdxi[0][0];
          break;
        case 5:
          bases1D[0] = square_number_column[basisnum / 2];
          bases1D[1] = square_number_row[basisnum / 2];
          bases1D[2] = basisnum % 2 * 2 + 1;
          if (basisnum % 2)
            coef *= dxdxi[2][1];
          break;
        }
    }
  // Interior
  else
    {
      unsigned int basisnum = i - 64 - 12*4*e;
      bases1D[0] = cube_number_column[basisnum] + 4;
      bases1D[1] = cube_number_row[basisnum] + 4;
      bases1D[2] = cube_number_page[basisnum] + 4;
    }

  // No singular elements
  libmesh_assert(coef);
  return coef;
}


} // end anonymous namespace



template <>
Real FE<3,HERMITE>::shape(const ElemType,
			  const Order,
			  const unsigned int,
			  const Point&)
{
  std::cerr << "Hermite elements require the real element\n"
	    << "to construct gradient-based degrees of freedom."
	    << std::endl;
  
  libmesh_error();
  return 0.;
}



template <>
Real FE<3,HERMITE>::shape(const Elem* elem,
			  const Order order,
			  const unsigned int i,
			  const Point& p)
{
  libmesh_assert (elem != NULL);

  hermite_compute_coefs(elem);

  const ElemType type = elem->type();
  
  const Order totalorder = static_cast<Order>(order + elem->p_level());
  
  switch (totalorder)
    {      
      // 3rd-order tricubic Hermite functions
    case THIRD:
      {
	switch (type)
	  {
	  case HEX8:
	  case HEX20:
	  case HEX27:
	    {
	      libmesh_assert (i<64);

              std::vector<unsigned int> bases1D;

              Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);

	      return coef *
                     FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) *
                     FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) *
                     FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
	    }
	  default:
            std::cerr << "ERROR: Unsupported element type!" << std::endl;
	    libmesh_error();
	  }
      }
      // by default throw an error
    default:
      std::cerr << "ERROR: Unsupported polynomial order!" << std::endl;
      libmesh_error();
    }
  
  libmesh_error();
  return 0.;
}



template <>
Real FE<3,HERMITE>::shape_deriv(const ElemType,
				const Order,			    
				const unsigned int,
				const unsigned int,
				const Point&)
{
  std::cerr << "Hermite elements require the real element\n"
	    << "to construct gradient-based degrees of freedom."
	    << std::endl;

  libmesh_error();
  return 0.;
}



template <>
Real FE<3,HERMITE>::shape_deriv(const Elem* elem,
				const Order order,
				const unsigned int i,
				const unsigned int j,
				const Point& p)
{
  libmesh_assert (elem != NULL);
  libmesh_assert (j == 0 || j == 1 || j == 2);

  hermite_compute_coefs(elem);

  const ElemType type = elem->type();
  
  const Order totalorder = static_cast<Order>(order + elem->p_level());
  
  switch (totalorder)
    {      
      // 3rd-order tricubic Hermite functions
    case THIRD:
      {
	switch (type)
	  {
	  case HEX8:
	  case HEX20:
	  case HEX27:
	    {
	      libmesh_assert (i<64);

              std::vector<unsigned int> bases1D;

              Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);

              switch (j) // Derivative type
		{
		case 0:
                  return coef *
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
                  break;
		case 1:
                  return coef *
                    FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
                  break;
		case 2:
                  return coef *
                    FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
                  break;
                }
                  
	    }
	  default:
            std::cerr << "ERROR: Unsupported element type!" << std::endl;
	    libmesh_error();
	  }
      }
      // by default throw an error
    default:
      std::cerr << "ERROR: Unsupported polynomial order!" << std::endl;
      libmesh_error();
    }
  
  libmesh_error();
  return 0.;
}



template <>
Real FE<3,HERMITE>::shape_second_deriv(const Elem* elem,
                                      const Order order,
                                      const unsigned int i,
                                      const unsigned int j,
                                      const Point& p)
{
  libmesh_assert (elem != NULL);

  hermite_compute_coefs(elem);

  const ElemType type = elem->type();
  
  const Order totalorder = static_cast<Order>(order + elem->p_level());
  
  switch (totalorder)
    {      
      // 3rd-order tricubic Hermite functions
    case THIRD:
      {
	switch (type)
	  {
	  case HEX8:
	  case HEX20:
	  case HEX27:
	    {
	      libmesh_assert (i<64);

              std::vector<unsigned int> bases1D;

              Real coef = hermite_bases_3D(bases1D, dxdxi, totalorder, i);

              switch (j) // Derivative type
		{
		case 0:
                  return coef *
                    FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
                  break;
		case 1:
                  return coef *
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
                  break;
		case 2:
                  return coef *
                    FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[2],p(2));
                  break;
		case 3:
                  return coef *
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
                  break;
		case 4:
                  return coef *
                    FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape_deriv(bases1D[2],p(2));
                  break;
		case 5:
                  return coef *
                    FEHermite<1>::hermite_raw_shape(bases1D[0],p(0)) * 
                    FEHermite<1>::hermite_raw_shape(bases1D[1],p(1)) * 
                    FEHermite<1>::hermite_raw_shape_second_deriv(bases1D[2],p(2));
                  break;
                }
                  
	    }
	  default:
            std::cerr << "ERROR: Unsupported element type!" << std::endl;
	    libmesh_error();
	  }
      }
      // by default throw an error
    default:
      std::cerr << "ERROR: Unsupported polynomial order!" << std::endl;
      libmesh_error();
    }
  
  libmesh_error();
  return 0.;
}