fe_clough_shape_2d.c

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

  Real d1xd2ndn = d1xd2ndx * N2x + d1xd2ndy * N2y;
  Real d1xd3ndn = d1xd3ndx * N3x + d1xd3ndy * N3y;
  Real d1yd2ndn = d1yd2ndx * N2x + d1yd2ndy * N2y;
  Real d1yd3ndn = d1yd3ndx * N3x + d1yd3ndy * N3y;

  Real d2nd2ndn = d2nd2ndx * N2x + d2nd2ndy * N2y;
  Real d2xd3ndn = d2xd3ndx * N3x + d2xd3ndy * N3y;
  Real d2xd1ndn = d2xd1ndx * N1x + d2xd1ndy * N1y;
  Real d2yd3ndn = d2yd3ndx * N3x + d2yd3ndy * N3y;
  Real d2yd1ndn = d2yd1ndx * N1x + d2yd1ndy * N1y;

  Real d3nd3ndn = d3nd3ndx * N3x + d3nd3ndy * N3y;
  Real d3xd1ndn = d3xd1ndx * N1x + d3xd1ndy * N1y;
  Real d3xd2ndn = d3xd2ndx * N2x + d3xd2ndy * N2y;
  Real d3yd1ndn = d3yd1ndx * N1x + d3yd1ndy * N1y;
  Real d3yd2ndn = d3yd2ndx * N2x + d3yd2ndy * N2y;

  // Calculate midpoint scaling factors
  
  d1nd1n = 1. / d1nd1ndn;
  d2nd2n = 1. / d2nd2ndn;
  d3nd3n = 1. / d3nd3ndn;

  // Calculate midpoint derivative adjustments to nodal value
  // interpolant functions

  d1d2n = -(d1d2ndx * N2x + d1d2ndy * N2y) / d2nd2ndn;
  d1d3n = -(d1d3ndx * N3x + d1d3ndy * N3y) / d3nd3ndn;
  d2d3n = -(d2d3ndx * N3x + d2d3ndy * N3y) / d3nd3ndn;
  d2d1n = -(d2d1ndx * N1x + d2d1ndy * N1y) / d1nd1ndn;
  d3d1n = -(d3d1ndx * N1x + d3d1ndy * N1y) / d1nd1ndn;
  d3d2n = -(d3d2ndx * N2x + d3d2ndy * N2y) / d2nd2ndn;

  // Calculate nodal derivative scaling factors

  d1xd1x = 1. / (d1xd1dx - d1xd1dy * d1yd1dx / d1yd1dy);
  d1xd1y = 1. / (d1yd1dx - d1xd1dx * d1yd1dy / d1xd1dy);
//  d1xd1y = - d1xd1x * (d1xd1dy / d1yd1dy);
  d1yd1y = 1. / (d1yd1dy - d1yd1dx * d1xd1dy / d1xd1dx);
  d1yd1x = 1. / (d1xd1dy - d1yd1dy * d1xd1dx / d1yd1dx);
//  d1yd1x = - d1yd1y * (d1yd1dx / d1xd1dx);
  d2xd2x = 1. / (d2xd2dx - d2xd2dy * d2yd2dx / d2yd2dy);
  d2xd2y = 1. / (d2yd2dx - d2xd2dx * d2yd2dy / d2xd2dy);
//  d2xd2y = - d2xd2x * (d2xd2dy / d2yd2dy);
  d2yd2y = 1. / (d2yd2dy - d2yd2dx * d2xd2dy / d2xd2dx);
  d2yd2x = 1. / (d2xd2dy - d2yd2dy * d2xd2dx / d2yd2dx);
//  d2yd2x = - d2yd2y * (d2yd2dx / d2xd2dx);
  d3xd3x = 1. / (d3xd3dx - d3xd3dy * d3yd3dx / d3yd3dy);
  d3xd3y = 1. / (d3yd3dx - d3xd3dx * d3yd3dy / d3xd3dy);
//  d3xd3y = - d3xd3x * (d3xd3dy / d3yd3dy);
  d3yd3y = 1. / (d3yd3dy - d3yd3dx * d3xd3dy / d3xd3dx);
  d3yd3x = 1. / (d3xd3dy - d3yd3dy * d3xd3dx / d3yd3dx);
//  d3yd3x = - d3yd3y * (d3yd3dx / d3xd3dx);

//  std::cerr << d1xd1dx << ' ';
//  std::cerr << d1xd1dy << ' ';
//  std::cerr << d1yd1dx << ' ';
//  std::cerr << d1yd1dy << ' ';
//  std::cerr << d2xd2dx << ' ';
//  std::cerr << d2xd2dy << ' ';
//  std::cerr << d2yd2dx << ' ';
//  std::cerr << d2yd2dy << ' ';
//  std::cerr << d3xd3dx << ' ';
//  std::cerr << d3xd3dy << ' ';
//  std::cerr << d3yd3dx << ' ';
//  std::cerr << d3yd3dy << std::endl;

  // Calculate midpoint derivative adjustments to nodal derivative
  // interpolant functions
  
  d1xd2n = -(d1xd1x * d1xd2ndn + d1xd1y * d1yd2ndn) / d2nd2ndn;
  d1yd2n = -(d1yd1y * d1yd2ndn + d1yd1x * d1xd2ndn) / d2nd2ndn;
  d1xd3n = -(d1xd1x * d1xd3ndn + d1xd1y * d1yd3ndn) / d3nd3ndn;
  d1yd3n = -(d1yd1y * d1yd3ndn + d1yd1x * d1xd3ndn) / d3nd3ndn;
  d2xd3n = -(d2xd2x * d2xd3ndn + d2xd2y * d2yd3ndn) / d3nd3ndn;
  d2yd3n = -(d2yd2y * d2yd3ndn + d2yd2x * d2xd3ndn) / d3nd3ndn;
  d2xd1n = -(d2xd2x * d2xd1ndn + d2xd2y * d2yd1ndn) / d1nd1ndn;
  d2yd1n = -(d2yd2y * d2yd1ndn + d2yd2x * d2xd1ndn) / d1nd1ndn;
  d3xd1n = -(d3xd3x * d3xd1ndn + d3xd3y * d3yd1ndn) / d1nd1ndn;
  d3yd1n = -(d3yd3y * d3yd1ndn + d3yd3x * d3xd1ndn) / d1nd1ndn;
  d3xd2n = -(d3xd3x * d3xd2ndn + d3xd3y * d3yd2ndn) / d2nd2ndn;
  d3yd2n = -(d3yd3y * d3yd2ndn + d3yd3x * d3xd2ndn) / d2nd2ndn;

  // Cross your fingers
//  std::cerr << d1nd1ndn << ' ';
//  std::cerr << d2xd1ndn << ' ';
//  std::cerr << d2yd1ndn << ' ';
//  std::cerr << std::endl;

//  std::cerr << "Transform variables: ";
//  std::cerr << d1nd1n << ' ';
//  std::cerr << d2nd2n << ' ';
//  std::cerr << d3nd3n << ' ';
//  std::cerr << d1d2n << ' ';
//  std::cerr << d1d3n << ' ';
//  std::cerr << d2d3n << ' ';
//  std::cerr << d2d1n << ' ';
//  std::cerr << d3d1n << ' ';
//  std::cerr << d3d2n << std::endl;
//  std::cerr << d1xd1x << ' ';
//  std::cerr << d1xd1y << ' ';
//  std::cerr << d1yd1x << ' ';
//  std::cerr << d1yd1y << ' ';
//  std::cerr << d2xd2x << ' ';
//  std::cerr << d2xd2y << ' ';
//  std::cerr << d2yd2x << ' ';
//  std::cerr << d2yd2y << ' ';
//  std::cerr << d3xd3x << ' ';
//  std::cerr << d3xd3y << ' ';
//  std::cerr << d3yd3x << ' ';
//  std::cerr << d3yd3y << std::endl;
//  std::cerr << d1xd2n << ' ';
//  std::cerr << d1yd2n << ' ';
//  std::cerr << d1xd3n << ' ';
//  std::cerr << d1yd3n << ' ';
//  std::cerr << d2xd3n << ' ';
//  std::cerr << d2yd3n << ' ';
//  std::cerr << d2xd1n << ' ';
//  std::cerr << d2yd1n << ' ';
//  std::cerr << d3xd1n << ' ';
//  std::cerr << d3yd1n << ' ';
//  std::cerr << d3xd2n << ' ';
//  std::cerr << d3yd2n << ' ';
//  std::cerr << std::endl;
//  std::cerr << std::endl;
}


unsigned char subtriangle_lookup(const Point& p)
{
  if ((p(0) >= p(1)) && (p(0) + 2 * p(1) <= 1))
    return 0;
  if ((p(0) <= p(1)) && (p(1) + 2 * p(0) <= 1))
    return 2;
  return 1;
}

  // Return shape function second derivatives on the unit right
  // triangle
Real clough_raw_shape_second_deriv(const unsigned int basis_num,
                            const unsigned int deriv_type,
                            const Point& p)
{
  Real xi = p(0), eta = p(1);

  switch (deriv_type)
  {

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

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