fe_clough_shape_2d.c

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

              return 1 - 3*xi + 2*xi*xi
                - eta + 4*xi*eta - 7./2.*eta*eta;
            case 1:
              return 5./2. - 4*xi + 3./2.*xi*xi
                - 9.*eta + 8*xi*eta + 13./2.*eta*eta;
            case 2:
              return 1 - 1./2.*xi*xi - 4*eta + 3*eta*eta;
          }
      case 5:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return - 1./2.*xi*eta + 1./4.*eta*eta;
            case 1:
              return 3./4. - 5./2.*xi + 7./4.*xi*xi
                - 2*eta + 7./2.*xi*eta + 5./4.*eta*eta;
            case 2:
              return - 1./4.*xi*xi;
          }
      case 6:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -xi + 2*xi*xi
                + eta + 7./2.*xi*eta - 13./4.*eta*eta;
            case 1:
              return - 11./4. + 19./2.*xi - 23./4.*xi*xi
                + 7*eta - 25./2.*xi*eta - 17./4.*eta*eta;
            case 2:
              return 9./4.*xi*xi;
          }
      case 7:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -eta + 9./2.*xi*eta + 5./4.*eta*eta;
            case 1:
              return - 13./4. + 19./2.*xi - 25./4.*xi*xi
                + 9*eta - 23./2.*xi*eta - 23./4.*eta*eta;
            case 2:
              return - xi + 7./4.*xi*xi + 4*xi*eta;
          }
      case 8:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return -2*eta - 1./2.*xi*eta + 13./4.*eta*eta;
            case 1:
              return 3./4. - 5./2.*xi + 7./4.*xi*xi
                - 4*eta + 7./2.*xi*eta + 17./4.*eta*eta;
            case 2:
              return - 1./4.*xi*xi
                - 2*eta + 3*eta*eta;
          }
      case 9:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return std::sqrt(2.) * (2*xi*eta - eta*eta);
            case 1:
              return std::sqrt(2.) * (- 3 + 10*xi - 7*xi*xi
                + 8*eta - 14*xi*eta - 5*eta*eta);
            case 2:
              return std::sqrt(2.) * (xi*xi);
          }
      case 10:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 4*eta - 4*xi*eta - 8*eta*eta;
            case 1:
              return 4 - 8*xi + 4*xi*xi
                - 12*eta + 12*xi*eta + 8*eta*eta;
            case 2:
              return 4*xi - 4*xi*xi
                - 8*xi*eta;
          }
      case 11:
        switch (subtriangle_lookup(p))
          {
            case 0:
              return 4*xi - 4*xi*xi
                - 4*eta - 8*xi*eta + 10.*eta*eta;
            case 1:
              return 2 - 8*xi + 6*xi*xi
                - 4*eta + 8*xi*eta + 2*eta*eta;
            case 2:
              return - 2*xi*xi;
          }
    }
  }

  libmesh_error();
  return 0.;
}

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

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

  libmesh_error();
  return 0.;
}

  
} // end anonymous namespace



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



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

  clough_compute_coefs(elem);

  const ElemType type = elem->type();
  
  const Order totalorder = static_cast<Order>(order + elem->p_level());

  switch (totalorder)
    {      
      // 2nd-order restricted Clough-Tocher element
    case SECOND:
      {
	switch (type)
	  {
	    // C1 functions on the Clough-Tocher triangle.
	  case TRI6:
	    {
	      libmesh_assert (i<9);
            // FIXME: it would be nice to calculate (and cache)
            // clough_raw_shape(j,p) only once per triangle, not 1-7
            // times
	      switch (i)
		{
	    // Note: these DoF numbers are "scrambled" because my
	    // initial numbering conventions didn't match libMesh
		case 0:
		  return clough_raw_shape(0, p)
                    + d1d2n * clough_raw_shape(10, p)
                    + d1d3n * clough_raw_shape(11, p);
		case 3:
		  return clough_raw_shape(1, p)
                    + d2d3n * clough_raw_shape(11, p)
                    + d2d1n * clough_raw_shape(9, p);
		case 6:
		  return clough_raw_shape(2, p)
                    + d3d1n * clough_raw_shape(9, p)
                    + d3d2n * clough_raw_shape(10, p);
                case 1:
                  return d1xd1x * clough_raw_shape(3, p)
                    + d1xd1y * clough_raw_shape(4, p)
                    + d1xd2n * clough_raw_shape(10, p)
                    + d1xd3n * clough_raw_shape(11, p)
                    + 0.5 * N01x * d3nd3n * clough_raw_shape(11, p)
                    + 0.5 * N02x * d2nd2n * clough_raw_shape(10, p);
		case 2:
                  return d1yd1y * clough_raw_shape(4, p)
                    + d1yd1x * clough_raw_shape(3, p)
                    + d1yd2n * clough_raw_shape(10, p)
                    + d1yd3n * clough_raw_shape(11, p)
                    + 0.5 * N01y * d3nd3n * clough_raw_shape(11, p)
                    + 0.5 * N02y * d2nd2n * clough_raw_shape(10, p);
		case 4:
                  return d2xd2x * clough_raw_shape(5, p)
                    + d2xd2y * clough_raw_shape(6, p)
                    + d2xd3n * clough_raw_shape(11, p)
                    + d2xd1n * clough_raw_shape(9, p)
                    + 0.5 * N10x * d3nd3n * clough_raw_shape(11, p)
                    + 0.5 * N12x * d1nd1n * clough_raw_shape(9, p);
		case 5:
                  return d2yd2y * clough_raw_shape(6, p)
                    + d2yd2x * clough_raw_shape(5, p)
                    + d2yd3n * clough_raw_shape(11, p)
                    + d2yd1n * clough_raw_shape(9, p)
                    + 0.5 * N10y * d3nd3n * clough_raw_shape(11, p)
                    + 0.5 * N12y * d1nd1n * clough_raw_shape(9, p);
		case 7:
                  return d3xd3x * clough_raw_shape(7, p)
                    + d3xd3y * clough_raw_shape(8, p)
                    + d3xd1n * clough_raw_shape(9, p)
                    + d3xd2n * clough_raw_shape(10, p)
                    + 0.5 * N20x * d2nd2n * clough_raw_shape(10, p)
                    + 0.5 * N21x * d1nd1n * clough_raw_shape(9, p);
		case 8:
                  return d3yd3y * clough_raw_shape(8, p)
                    + d3yd3x * clough_raw_shape(7, p)
                    + d3yd1n * clough_raw_shape(9, p)
                    + d3yd2n * clough_raw_shape(10, p)
                    + 0.5 * N20y * d2nd2n * clough_raw_shape(10, p)
                    + 0.5 * N21y * d1nd1n * clough_raw_shape(9, p);