fe_bernstein_shape_3d.c

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

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

		    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;
			  }
		      }
		  }