fe_hierarchic_shape_3d.c

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

	    {
	      // Case 1
	      xi   = xi_saved;
	      zeta = zeta_saved;
	    }
	  else
	    {
	      // Case 2
	      xi   = zeta_saved;
	      zeta = xi_saved;
	    }

	else if (elem->point(7) == min_point)
	  if (elem->point(3) == std::min(elem->point(3), elem->point(6)))
	    {
	      // Case 3
	      xi   = -zeta_saved;
	      zeta = xi_saved;
	    }
	  else
	    {
	      // Case 4
	      xi   = xi_saved;
	      zeta = -zeta_saved;
	    }

	else if (elem->point(6) == min_point)
	  if (elem->point(7) == std::min(elem->point(7), elem->point(2)))
	    {
	      // Case 5
	      xi   = -xi_saved;
	      zeta = -zeta_saved;
	    }
	  else
	    {
	      // Case 6
	      xi   = -zeta_saved;
	      zeta = -xi_saved;
	    }

	else if (elem->point(2) == min_point)
	  {
	    if (elem->point(6) == std::min(elem->point(3), elem->point(6)))
	      {
		// Case 7
		xi   = zeta_saved;
		zeta = -xi_saved;
	      }
	    else
	      {
		// Case 8
		xi   = -xi_saved;
	      zeta = zeta_saved;
	      }
	  }
      }
    // Face 4
    else if (i < 8 + 12*e + 5*e*e)
      {
	unsigned int basisnum = i - 8 - 12*e - 4*e*e;
	i0 = 0;
	i1 = square_number_row[basisnum] + 2;
	i2 = square_number_column[basisnum] + 2;
	const Point min_point = get_min_point(elem, 3, 0, 4, 7);

	if (elem->point(0) == min_point)
	  if (elem->point(3) == std::min(elem->point(3), elem->point(4)))
	    {
	      // Case 1
	      eta  = eta_saved;
	      zeta = zeta_saved;
	    }
	  else
	    {
	      // Case 2
	      eta  = zeta_saved;
	      zeta = eta_saved;
	    }

	else if (elem->point(4) == min_point)
	  if (elem->point(0) == std::min(elem->point(0), elem->point(7)))
	    {
	      // Case 3
	      eta  = -zeta_saved;
	      zeta = eta_saved;
	    }
	  else
	    {
	      // Case 4
	      eta  = eta_saved;
	      zeta = -zeta_saved;
	    }

	else if (elem->point(7) == min_point)
	  if (elem->point(4) == std::min(elem->point(4), elem->point(3)))
	    {
	      // Case 5
	      eta  = -eta_saved;
	      zeta = -zeta_saved;
	    }
	  else
	    {
	      // Case 6
	      eta  = -zeta_saved;
	      zeta = -eta_saved;
	    }

	else if (elem->point(3) == min_point)
	  {
	    if (elem->point(7) == std::min(elem->point(7), elem->point(0)))
	      {
		// Case 7
		eta   = zeta_saved;
		zeta = -eta_saved;
	      }
	    else
	      {
		// Case 8
		eta  = -eta_saved;
		zeta = zeta_saved;
	      }
	  }
      }
    // Face 5
    else if (i < 8 + 12*e + 6*e*e)
      {
	unsigned int basisnum = i - 8 - 12*e - 5*e*e;
	i0 = square_number_row[basisnum] + 2;
	i1 = square_number_column[basisnum] + 2;
	i2 = 1;
	const Point min_point = get_min_point(elem, 4, 5, 6, 7);

	if (elem->point(4) == min_point)
	  if (elem->point(5) == std::min(elem->point(5), elem->point(7)))
	    {
	      // Case 1
	      xi  = xi_saved;
	      eta = eta_saved;
	    }
	  else
	    {
	      // Case 2
	      xi  = eta_saved;
	      eta = xi_saved;
	    }

	else if (elem->point(5) == min_point)
	  if (elem->point(6) == std::min(elem->point(6), elem->point(4)))
	    {
	      // Case 3
	      xi  = eta_saved;
	      eta = -xi_saved;
	    }
	  else
	    {
	      // Case 4
	      xi  = -xi_saved;
	      eta = eta_saved;
	    }

	else if (elem->point(6) == min_point)
	  if (elem->point(7) == std::min(elem->point(7), elem->point(5)))
	    {
	      // Case 5
	      xi  = -xi_saved;
	      eta = -eta_saved;
	    }
	  else
	    {
	      // Case 6
	      xi  = -eta_saved;
	      eta = -xi_saved;
	    }

	else if (elem->point(7) == min_point)
	  {
	    if (elem->point(4) == std::min(elem->point(4), elem->point(6)))
	      {
		// Case 7
		xi  = -eta_saved;
		eta = xi_saved;
	      }
	    else
	      {
		// Case 8
		xi  = xi_saved;
		eta = eta_saved;
	      }
	  }
      }

    // Internal DoFs
    else 
      {
	unsigned int basisnum = i - 8 - 12*e - 6*e*e;
	i0 = cube_number_column[basisnum] + 2;
	i1 = cube_number_row[basisnum] + 2;
	i2 = cube_number_page[basisnum] + 2;
      }
  }
} // end anonymous namespace




template <>
Real FE<3,HIERARCHIC>::shape(const ElemType,
                             const Order,
                             const unsigned int,
                             const Point&)
{
  std::cerr << "Hierarchic polynomials require the element type\n"
            << "because edge and face orientation is needed."
            << std::endl;
  
  libmesh_error();
  return 0.;
}



template <>
Real FE<3,HIERARCHIC>::shape(const Elem* elem,
                             const Order order,
                             const unsigned int i,
                             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());
  
  switch (type)
    {
    case HEX8:
    case HEX20:
      libmesh_assert(totalorder < 2);
    case HEX27:
      {
        libmesh_assert (i<(totalorder+1u)*(totalorder+1u)*(totalorder+1u));
        
        // Compute hex shape functions as a tensor-product
        Real xi   = p(0);
        Real eta  = p(1);
        Real zeta = p(2);
        
        unsigned int i0, i1, i2;
        
	cube_indices(elem, totalorder, i, xi, eta, zeta, i0, i1, i2);

        return (FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i0, xi)*
                FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i1, eta)*
                FE<1,HIERARCHIC>::shape(EDGE3, totalorder, i2, zeta));
      }

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




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



template <>
Real FE<3,HIERARCHIC>::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);

  libmesh_assert (j < 3);
  
  // cheat by using finite difference approximations:
  const Real eps = 1.e-6;
  Point pp, pm;
        
  switch (j)
  {
    // d()/dxi
  case 0:
    {
      pp = Point(p(0)+eps, p(1), p(2));
      pm = Point(p(0)-eps, p(1), p(2));
      break;
    }

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

    // d()/dzeta
  case 2:
    {
      pp = Point(p(0), p(1), p(2)+eps);
      pm = Point(p(0), p(1), p(2)-eps);
      break;
    }

  default:
    libmesh_error();
  }

  return (FE<3,HIERARCHIC>::shape(elem, order, i, pp) -
          FE<3,HIERARCHIC>::shape(elem, order, i, pm))/2./eps;
#endif
  
  libmesh_error();
  return 0.;
}



template <>
Real FE<3,HIERARCHIC>::shape_second_deriv(const ElemType,
                                          const Order,
                                          const unsigned int,
                                          const unsigned int,
                                          const Point& )
{
  std::cerr << "Hierarchic polynomials require the element type\n"
            << "because edge and face orientation is needed."
            << std::endl;
  libmesh_error();
  
  return 0.;
}



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

  const Real eps = 1.e-6;
  Point pp, pm;
  unsigned int prevj = libMesh::invalid_uint;

  switch (j)
  {
    //  d^2()/dxi^2
    case 0:
      {
        pp = Point(p(0)+eps, p(1), p(2));
        pm = Point(p(0)-eps, p(1), p(2));
        prevj = 0;
        break;
      }

    //  d^2()/dxideta
    case 1:
      {
        pp = Point(p(0), p(1)+eps, p(2));
        pm = Point(p(0), p(1)-eps, p(2));
        prevj = 0;
        break;
      }

    //  d^2()/deta^2
    case 2:
      {
        pp = Point(p(0), p(1)+eps, p(2));
        pm = Point(p(0), p(1)-eps, p(2));
        prevj = 1;
        break;
      }

    //  d^2()/dxidzeta
    case 3:
      {
        pp = Point(p(0), p(1), p(2)+eps);
        pm = Point(p(0), p(1), p(2)-eps);
        prevj = 0;
        break;
      }

    //  d^2()/detadzeta
    case 4:
      {
        pp = Point(p(0), p(1), p(2)+eps);
        pm = Point(p(0), p(1), p(2)-eps);
        prevj = 1;
        break;
      }

    //  d^2()/dzeta^2
    case 5:
      {
        pp = Point(p(0), p(1), p(2)+eps);
        pm = Point(p(0), p(1), p(2)-eps);
        prevj = 2;
        break;
      }
    default:
      libmesh_error();
  }
  return (FE<3,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pp) -
          FE<3,HIERARCHIC>::shape_deriv(elem, order, i, prevj, pm))
          / 2. / eps;
}