fe_clough_shape_2d.c

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

            case 1:
              return - 8 + 8*xi + 12*eta;
            case 2:
              return 4 - 8*xi - 8*eta;
          }
      case 11:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 4 - 8*xi - 8*eta;
            case 1:
              return -8 + 12*xi + 8*eta;
            case 2:
              return -4*xi;
          }
    }

  // second derivative in eta-eta direction
  case 2:
  switch (basis_num)
    {
      case 0:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -6 - 6*xi + 18*eta;
            case 1:
              return -30 + 42*xi + 42*eta;
            case 2:
              return -6 + 12*eta;
          }
      case 1:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 3*xi - 3*eta;
            case 1:
              return 12 - 21*xi - 15*eta;
            case 2:
              return 0;
          }
      case 2:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 6 + 3*xi - 15*eta;
            case 1:
              return 18 - 21.*xi - 27*eta;
            case 2:
              return 6 - 12*eta;
          }
      case 3:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -3 - xi + 14*eta;
            case 1:
              return 1 + 3*xi - 2*eta;
            case 2:
              return 4*xi;
          }
      case 4:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -1 + 4*xi - 7*eta;
            case 1:
              return -9 + 8*xi + 13*eta;
            case 2:
              return -4 + 6*eta;
          }
      case 5:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return - 1./2.*xi + 1./2.*eta;
            case 1:
              return -2 + 7./2.*xi + 5./2.*eta;
            case 2:
              return 0;
          }
      case 6:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 1 + 7./2.*xi - 13./2.*eta;
            case 1:
              return 7 - 25./2.*xi - 17./2.*eta;
            case 2:
              return 0;
          }
      case 7:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -1 + 9./2.*xi + 5./2.*eta;
            case 1:
              return 9 - 23./2.*xi - 23./2.*eta;
            case 2:
              return 4*xi;
          }
      case 8:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -2 - 1./2.*xi + 13./2.*eta;
            case 1:
              return -4 + 7./2.*xi + 17./2.*eta;
            case 2:
              return -2 + 6*eta;
          }
      case 9:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return std::sqrt(2.) * (2*xi - 2*eta);
            case 1:
              return std::sqrt(2.) * (8 - 14*xi - 10*eta);
            case 2:
              return 0;
          }
      case 10:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 4 - 4*xi - 16*eta;
            case 1:
              return -12 + 12*xi + 16*eta;
            case 2:
              return -8*xi;
          }
      case 11:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -4 - 8*xi + 20*eta;
            case 1:
              return -4 + 8*xi + 4*eta;
            case 2:
              return 0;
          }
    }
  }

  libmesh_error();
  return 0.;
}



Real clough_raw_shape_deriv(const unsigned int basis_num,
                            const unsigned int deriv_type,
                            const Point& p)
{
  Real xi = p(0), eta = p(1);

  switch (deriv_type)
  {
  // derivative in xi direction
  case 0:
  switch (basis_num)
    {
      case 0:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -6*xi + 6*xi*xi
                - 3*eta*eta;
            case 1:
              return 9 - 30*xi + 21*xi*xi
                - 30*eta + 42*xi*eta
                + 21*eta*eta;
            case 2:
              return -6*xi + 9*xi*xi
                - 6*xi*eta;
          }
      case 1:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 6*xi - 6*xi*xi
                + 3./2.*eta*eta;
            case 1:
              return - 9./2. + 18*xi - 27./2.*xi*xi
                + 15*eta - 21*xi*eta
                - 21./2.*eta*eta;
            case 2:
              return 6*xi - 15./2.*xi*xi
                + 3*xi*eta;
          }
      case 2:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 3./2.*eta*eta;
            case 1:
              return - 9./2. + 12*xi - 15./2.*xi*xi
                + 15*eta - 21*xi*eta
                - 21./2.*eta*eta;
            case 2:
              return -3./2.*xi*xi
                + 3*xi*eta;
          }
      case 3:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 1 - 4*xi + 3*xi*xi
                - 1./2.*eta*eta;
            case 1:
              return 5./2. - 9*xi + 13./2.*xi*xi
                - 4*eta + 8*xi*eta
                + 3./2.*eta*eta;
            case 2:
              return 1 - xi - 7./2.*xi*xi
                - 3*eta + 4*xi*eta
                + 2*eta*eta;
          }
      case 4:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return - 3*eta + 4*xi*eta
                + 2*eta*eta;
            case 1:
              return xi - xi*xi
                - 4*eta + 3*xi*eta
                + 4*eta*eta;
            case 2:
              return -3*xi + 7*xi*xi
                - xi*eta;
          }
      case 5:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -2*xi + 3*xi*xi
                - 1./4.*eta*eta;
            case 1:
              return 3./4. - 4*xi + 17./4.*xi*xi
                - 5./2.*eta + 7./2.*xi*eta
                + 7./4.*eta*eta;
            case 2:
              return -2*xi + 13./4.*xi*xi
                - 1./2.*xi*eta;
          }
      case 6:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -eta + 4*xi*eta
                + 7./4.*eta*eta;
            case 1:
              return -13./4. + 9*xi - 23./4.*xi*xi
                + 19./2.*eta - 23./2.*xi*eta
                - 25./4.*eta*eta;
            case 2:
              return -xi + 5./4.*xi*xi
                + 9./2.*xi*eta;
          }
      case 7:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 9./4.*eta*eta;
            case 1:
              return - 11./4. + 7*xi - 17./4.*xi*xi
                + 19./2.*eta - 25./2.*xi*eta
                - 23./4.*eta*eta;
            case 2:
              return xi - 13./4.*xi*xi
                - eta + 7./2.*xi*eta + 2*eta*eta;
          }
      case 8:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return - 1./4.*eta*eta;
            case 1:
              return 3./4. - 2*xi + 5./4.*xi*xi
                - 5./2.*eta + 7./2.*xi*eta
                + 7./4.*eta*eta;
            case 2:
              return 1./4.*xi*xi
                - 1./2.*xi*eta;
          }
      case 9:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return std::sqrt(2.) * eta*eta;
            case 1:
              return std::sqrt(2.) * (-3 + 8*xi - 5*xi*xi
                + 10*eta - 14*xi*eta
                - 7*eta*eta);
            case 2:
              return std::sqrt(2.) * (-xi*xi + 2*xi*eta);
          }
      case 10:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -2*eta*eta;
            case 1:
              return 2 - 4*xi + 2*xi*xi
                - 8*eta + 8*xi*eta
                + 6*eta*eta;
            case 2:
              return -4*xi + 10*xi*xi
                + 4*eta - 8*xi*eta
                - 4*eta*eta;
          }
      case 11:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 4*eta - 8*xi*eta
                - 4*eta*eta;
            case 1:
              return 4 - 12*xi + 8*xi*xi
                - 8*eta + 12*xi*eta
                + 4*eta*eta;
            case 2:
              return 4*xi - 8*xi*xi
                - 4*xi*eta;
          }
    }

  // derivative in eta direction
  case 1:
  switch (basis_num)
    {
      case 0:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return - 6*eta - 6*xi*eta + 9*eta*eta;
            case 1:
              return 9 - 30*xi + 21*xi*xi
                - 30*eta + 42*xi*eta + 21*eta*eta;
            case 2:
              return - 3*xi*xi
                - 6*eta + 6*eta*eta;
          }
      case 1:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return + 3*xi*eta
                - 3./2.*eta*eta;
            case 1:
              return - 9./2. + 15*xi - 21./2.*xi*xi
                + 12*eta - 21*xi*eta - 15./2.*eta*eta;
            case 2:
              return + 3./2.*xi*xi;
          }
      case 2:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 6*eta + 3*xi*eta - 15./2.*eta*eta;
            case 1:
              return - 9./2. + 15*xi - 21./2.*xi*xi
                + 18*eta - 21.*xi*eta - 27./2.*eta*eta;
            case 2:
              return 3./2.*xi*xi
                + 6*eta - 6*eta*eta;
          }
      case 3:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return - 3*eta - xi*eta + 7*eta*eta;
            case 1:
              return - 4*xi + 4*xi*xi
                + eta + 3*xi*eta - eta*eta;
            case 2:
              return - 3*xi + 2*xi*xi
                + 4*xi*eta;
          }
      case 4:
        switch (subtriangle_lookup(p))
          {
            case 0: