fe_bernstein_shape_3d.c

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

		    else if (elem->point(3) == min_point)
		      if (elem->point(0) == std::min(elem->point(0), elem->point(2)))
			{
			  // Case 3
			  xi_mapped  = -eta;
			  eta_mapped = xi;
			}
		      else
			{
			  // Case 4
			  xi_mapped  = xi;
			  eta_mapped = -eta;
			}
		    
		    else if (elem->point(2) == min_point)
		      if (elem->point(3) == std::min(elem->point(3), elem->point(1)))
			{
			  // Case 5
			  xi_mapped  = -xi;
			  eta_mapped = -eta;
			}
		      else
			{
			  // Case 6
			  xi_mapped  = -eta;
			  eta_mapped = -xi;
			}
		    
		    else if (elem->point(1) == min_point)
		      {
			if (elem->point(2) == std::min(elem->point(2), elem->point(0)))
			  {
			    // Case 7
			    xi_mapped  = eta;
			    eta_mapped = -xi;
			  }
			else
			  {
			    // Case 8
			    xi_mapped  = -xi;
			    eta_mapped = eta;
			  }
		      }
		  }
		
		
		// Face 1
		else if ((i1[i] == 0) && (i0[i] >= 2) && (i2[i] >= 2))
		  {
		    const Point min_point = std::min(elem->point(0),
							   std::min(elem->point(1),
								    std::min(elem->point(5),
									     elem->point(4))));
		    if (elem->point(0) == min_point)
		      if (elem->point(1) == std::min(elem->point(1), elem->point(4)))
			{
			  // Case 1
			  xi_mapped   = xi;
			  zeta_mapped = zeta;
			}
		      else
			{
			  // Case 2
			  xi_mapped   = zeta;
			  zeta_mapped = xi;
			}

		    else if (elem->point(1) == min_point)
		      if (elem->point(5) == std::min(elem->point(5), elem->point(0)))
			{
			  // Case 3
			  xi_mapped   = zeta;
			  zeta_mapped = -xi;
			}
		      else
			{
			  // Case 4
			  xi_mapped   = -xi;
			  zeta_mapped = zeta;
			}
		    
		    else if (elem->point(5) == min_point)
		      if (elem->point(4) == std::min(elem->point(4), elem->point(1)))
			{
			  // Case 5
			  xi_mapped   = -xi;
			  zeta_mapped = -zeta;
			}
		      else
			{
			  // Case 6
			  xi_mapped   = -zeta;
			  zeta_mapped = -xi;
			}
		    
		    else if (elem->point(4) == min_point)
		      {
			if (elem->point(0) == std::min(elem->point(0), elem->point(5)))
			  {
			    // Case 7
			    xi_mapped   = -xi;
			    zeta_mapped = zeta;
			  }
			else
			  {
			    // Case 8
			    xi_mapped   = xi;
			    zeta_mapped = -zeta;
			  }
		      }
		  }
		
		
		// Face 2
		else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] >= 2))
		  {
		    const Point min_point = std::min(elem->point(1),
							   std::min(elem->point(2),
								    std::min(elem->point(6),
									     elem->point(5))));
		    if (elem->point(1) == min_point)
		      if (elem->point(2) == std::min(elem->point(2), elem->point(5)))
			{
			  // Case 1
			  eta_mapped  = eta;
			  zeta_mapped = zeta;
			}
		      else
			{
			  // Case 2
			  eta_mapped  = zeta;
			  zeta_mapped = eta;
			}
		    
		    else if (elem->point(2) == min_point)
		      if (elem->point(6) == std::min(elem->point(6), elem->point(1)))
			{
			  // Case 3
			  eta_mapped  = zeta;
			  zeta_mapped = -eta;
			}
		      else
			{
			  // Case 4
			  eta_mapped  = -eta;
			  zeta_mapped = zeta;
			}
		    
		    else if (elem->point(6) == min_point)
		      if (elem->point(5) == std::min(elem->point(5), elem->point(2)))
			{
			  // Case 5
			  eta_mapped  = -eta;
			  zeta_mapped = -zeta;
			}
		      else
			{
			  // Case 6
			  eta_mapped  = -zeta;
			  zeta_mapped = -eta;
			}

		    else if (elem->point(5) == min_point)
		      {
			if (elem->point(1) == std::min(elem->point(1), elem->point(6)))
			  {
			    // Case 7
			    eta_mapped  = -zeta;
			    zeta_mapped = eta;
			  }
			else
			  {
			    // Case 8
			    eta_mapped   = eta;
			    zeta_mapped = -zeta;
			  }
		      }
		  }

		
		// Face 3
		else if ((i1[i] == 1) && (i0[i] >= 2) && (i2[i] >= 2))
		  {
 		    const Point min_point = std::min(elem->point(2),
							   std::min(elem->point(3),
								    std::min(elem->point(7),
									     elem->point(6))));
		    if (elem->point(3) == min_point)
		      if (elem->point(2) == std::min(elem->point(2), elem->point(7)))
			{
			  // Case 1
			  xi_mapped   = xi;
			  zeta_mapped = zeta;
			}
		      else
			{
			  // Case 2
			  xi_mapped   = zeta;
			  zeta_mapped = xi;
			}
		    
		    else if (elem->point(7) == min_point)
		      if (elem->point(3) == std::min(elem->point(3), elem->point(6)))
			{
			  // Case 3
			  xi_mapped   = -zeta;
			  zeta_mapped = xi;
			}
		      else
			{
			  // Case 4
			  xi_mapped   = xi;
			  zeta_mapped = -zeta;
			}
		    
		    else if (elem->point(6) == min_point)
		      if (elem->point(7) == std::min(elem->point(7), elem->point(2)))
			{
			  // Case 5
			  xi_mapped   = -xi;
			  zeta_mapped = -zeta;
			}
		      else
			{
			  // Case 6
			  xi_mapped   = -zeta;
			  zeta_mapped = -xi;
			}

		    else if (elem->point(2) == min_point)
		      {
			if (elem->point(6) == std::min(elem->point(3), elem->point(6)))
			  {
			    // Case 7
			    xi_mapped   = zeta;
			    zeta_mapped = -xi;
			  }
			else
			  {
			    // Case 8
			    xi_mapped   = -xi;
			    zeta_mapped = zeta;
			  }
		      }
		  }
		
		
		// Face 4
		else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] >= 2))
		  {
 		    const Point min_point = std::min(elem->point(3),
							   std::min(elem->point(0),
								    std::min(elem->point(4),
									     elem->point(7))));
		    if (elem->point(0) == min_point)
		      if (elem->point(3) == std::min(elem->point(3), elem->point(4)))
			{
			  // Case 1
			  eta_mapped  = eta;
			  zeta_mapped = zeta;
			}
		      else
			{
			  // Case 2
			  eta_mapped  = zeta;
			  zeta_mapped = eta;
			}

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

		    else if (elem->point(3) == min_point)
		      {
			if (elem->point(7) == std::min(elem->point(7), elem->point(0)))
			  {
			    // Case 7
			    eta_mapped   = zeta;
			    zeta_mapped = -eta;
			  }
			else
			  {
			    // Case 8
			    eta_mapped  = -eta;
			    zeta_mapped = zeta;
			  }
		      }
		  }
		
		
		// Face 5
		else if ((i2[i] == 1) && (i0[i] >= 2) && (i1[i] >= 2))
		  {
 		    const Point min_point = std::min(elem->point(4),
							   std::min(elem->point(5),
								    std::min(elem->point(6),
									     elem->point(7))));
		    if (elem->point(4) == min_point)
		      if (elem->point(5) == std::min(elem->point(5), elem->point(7)))
			{
			  // Case 1
			  xi_mapped  = xi;
			  eta_mapped = eta;
			}
		      else
			{
			  // Case 2
			  xi_mapped  = eta;
			  eta_mapped = xi;
			}
		    
		    else if (elem->point(5) == min_point)
		      if (elem->point(6) == std::min(elem->point(6), elem->point(4)))
			{
			  // Case 3
			  xi_mapped  = eta;
			  eta_mapped = -xi;
			}
		      else
			{
			  // Case 4
			  xi_mapped  = -xi;
			  eta_mapped = eta;
			}
		    
		    else if (elem->point(6) == min_point)
		      if (elem->point(7) == std::min(elem->point(7), elem->point(5)))
			{
			  // Case 5
			  xi_mapped  = -xi;
			  eta_mapped = -eta;
			}
		      else
			{
			  // Case 6
			  xi_mapped  = -eta;
			  eta_mapped = -xi;
			}
		    
		    else if (elem->point(7) == min_point)
		      {
			if (elem->point(4) == std::min(elem->point(4), elem->point(6)))
			  {
			    // Case 7
			    xi_mapped  = -eta;
			    eta_mapped = xi;
			  }
			else
			  {
			    // Case 8
			    xi_mapped  = xi;
			    eta_mapped = eta;
			  }
		      }
		  }
	      }
	      
	      return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi_mapped)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta_mapped)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta_mapped));
	    }
	    
	    
	  default:
	    libmesh_error();
	    
	  } //case HEX27
 
      }	//case THIRD
      
      
      // 4th-order Bernstein.
    case FOURTH:
      {
	switch (type)
	  {  
	    
	    // Bernstein shape functions on the hexahedral.
	  case HEX27:
	    {
	      libmesh_assert (i<125);
	      
	      // Compute hex shape functions as a tensor-product
	      const Real xi    = p(0);
	      const Real eta   = p(1);
	      const Real zeta  = p(2);
	      Real xi_mapped   = p(0);
	      Real eta_mapped  = p(1);
	      Real zeta_mapped = p(2);
	      
	      // The only way to make any sense of this
	      // is to look at the mgflo/mg2/mgf documentation
	      // and make the cut-out cube!
	      //  Nodes                         0  1  2  3  4  5  6  7  8  8  9  9 10 10 11 11 12 12 13 13 14 14 15 15 16 16 17 17 18 18 19 19 20 20 20 20 21 21 21 21 22 22 22 22 23 23 23 23 24 24 24 24 25 25 25 25 26 26 26 26 26 26 26 26 
	      //  DOFS                          0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 18 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 60 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
	      static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0, 0, 2, 3, 4, 1, 1, 1, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4};
	      static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 2, 3, 4, 2, 3, 4, 2, 3, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 2, 2, 2, 3, 3, 3, 4, 4, 4, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4, 4};
	      
	      
	      
	      // handle the edge orientation
	      {
		// Edge 0
		if ( (i0[i] >= 2) && (i1[i] == 0) && (i2[i] == 0))
		  {
		    if (elem->point(0) != std::min(elem->point(0), elem->point(1)))
		      xi_mapped = -xi;
		  }
		// Edge 1
		else if ((i0[i] == 1) && (i1[i] >= 2) && (i2[i] == 0))
		  {
		    if (elem->point(1) != std::min(elem->point(1), elem->point(2)))
		      eta_mapped = -eta;
		  }
		// Edge 2
		else if ((i0[i] >= 2) && (i1[i] == 1) && (i2[i] == 0))
		  {
		    if (elem->point(3) != std::min(elem->point(3), elem->point(2)))
		      xi_mapped = -xi;
		  }
		// Edge 3
		else if ((i0[i] == 0) && (i1[i] >= 2) && (i2[i] == 0))
		  {
		    if (elem->point(0) != std::min(elem->point(0), elem->point(3)))
		      eta_mapped = -eta;
		  }
		// Edge 4
		else if ((i0[i] == 0) && (i1[i] == 0) && (i2[i] >=2 ))