fe_bernstein_shape_3d.c

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

// $Id: fe_bernstein_shape_3D.C 2858 2008-06-13 18:59:07Z benkirk $

// The Next Great Finite Element Library.
// Copyright (C) 2002-2007  Benjamin S. Kirk, John W. Peterson
  
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
  
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA


// C++ includes


// Local includes
#include "libmesh_config.h"
#ifdef ENABLE_HIGHER_ORDER_SHAPES

#include "fe.h"
#include "elem.h"




template <>
Real FE<3,BERNSTEIN>::shape(const ElemType,
			    const Order,
			    const unsigned int,
			    const Point&)
{
  std::cerr << "Bernstein polynomials require the element type\n"
	    << "because edge and face orientation is needed."
	    << std::endl;
  
  libmesh_error();
  return 0.;
}



template <>
Real FE<3,BERNSTEIN>::shape(const Elem* elem,
			    const Order order,
			    const unsigned int i,
			    const Point& p)
{
  
#if DIM == 3
  
  libmesh_assert (elem != NULL);
  const ElemType type = elem->type();

  const Order totalorder = static_cast<Order>(order + elem->p_level());
  
  switch (totalorder)
    {
      
      // 1st order Bernstein.
    case FIRST:
      {
	switch (type)
	  {
	    
	    // Bernstein shape functions on the tetrahedron.
	  case TET4:
	  case TET10:
	    {
	      libmesh_assert(i<4);
	      
	      // Area coordinates
	      const Real zeta1 = p(0);
	      const Real zeta2 = p(1);
	      const Real zeta3 = p(2);
	      const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
	      
	      switch(i)
		{
		case  0:  return zeta0;	   
		case  1:  return zeta1;	   
		case  2:  return zeta2;	   
		case  3:  return zeta3;
		  
		default:
		  libmesh_error();
		}
	    }
	    
	    // Bernstein shape functions on the hexahedral.
	  case HEX20:
	  case HEX27:
	    {
	      libmesh_assert (i<8);
	      
	      // Compute hex shape functions as a tensor-product
	      const Real xi   = p(0);
	      const Real eta  = p(1);
	      const Real zeta = 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!
	      //                                0  1  2  3  4  5  6  7
	      static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1};
	      static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1};
	      
	      return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
	    }

	    
	  default:
	    libmesh_error();
	  }
	}
      
      
      
      
    case SECOND:
      {
	switch (type)
	  {
	
	    // Bernstein shape functions on the tetrahedron.
	  case TET10:
	    {
	      libmesh_assert(i<10);
	
	      // Area coordinates
	      const Real zeta1 = p(0);
	      const Real zeta2 = p(1);
	      const Real zeta3 = p(2);
	      const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
	      
	      switch(i)
		{
		case  0:  return zeta0*zeta0;			
		case  1:  return zeta1*zeta1;			
		case  2:  return zeta2*zeta2;        	
		case  3:  return zeta3*zeta3;
		case  4:  return 2.*zeta0*zeta1;
		case  5:  return 2.*zeta1*zeta2;
		case  6:  return 2.*zeta0*zeta2;
		case  7:  return 2.*zeta3*zeta0;
		case  8:  return 2.*zeta1*zeta3;
		case  9:  return 2.*zeta2*zeta3;     
		  
		default:
		  libmesh_error();
		}
	    }
	    
	    // Bernstein shape functions on the 20-noded hexahedral.
	  case HEX20:
	    {
	      libmesh_assert (i<20);
	      
	      // Compute hex shape functions as a tensor-product
	      const Real xi   = p(0);
	      const Real eta  = p(1);
	      const Real zeta = 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!
	      //                                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
	      static const unsigned int i0[] = {0,     1,     1,     0,     0,     1,     1,     0,     2,     1,     2,     0,     0,     1,     1,     0,     2,     1,     2,     0,  2, 2, 1, 2, 0, 2, 2};
	      static const unsigned int i1[] = {0,     0,     1,     1,     0,     0,     1,     1,     0,     2,     1,     2,     0,     0,     1,     1,     0,     2,     1,     2,  2, 0, 2, 1, 2, 2, 2};
	      static const unsigned int i2[] = {0,     0,     0,     0,     1,     1,     1,     1,     0,     0,     0,     0,     2,     2,     2,     2,     1,     1,     1,     1,  0, 2, 2, 2, 2, 1, 2};
	      //To compute the hex20 shape functions the original shape functions for hex27 are used.
	      //scalx[i] tells how often the original x-th shape function has to be added to the original i-th shape function
	      //to compute the new i-th shape function for hex20
	      //example: B_0^HEX20 = B_0^HEX27 - 0.25*B_20^HEX27 - 0.25*B_21^HEX27 + 0*B_22^HEX27 + 0*B_23^HEX27 - 0.25*B_24^HEX27 + 0*B_25^HEX27 - 0.25*B_26^HEX27
	      //         B_0^HEX20 = B_0^HEX27 + scal20[0]*B_20^HEX27 + scal21[0]*B_21^HEX27 + ...
	      static const Real scal20[] =     {-0.25, -0.25, -0.25, -0.25, 0,     0,     0,     0,     0.5,   0.5,   0.5,   0.5,   0,     0,     0,     0,     0,     0,     0,     0};
	      static const Real scal21[] =     {-0.25, -0.25, 0,     0,     -0.25, -0.25, 0,     0,     0.5,   0,     0,     0,     0.5,   0.5,   0,     0,     0.5,   0,     0,     0};
	      static const Real scal22[] =     {0,     -0.25, -0.25, 0,     0,     -0.25, -0.25, 0,     0,     0.5,   0,     0,     0,     0.5,   0.5,   0,     0,     0.5,   0,     0};
	      static const Real scal23[] =     {0,     0,     -0.25, -0.25, 0,     0,     -0.25, -0.25, 0,     0,     0.5,   0,     0,     0,     0.5,   0.5,   0,     0,     0.5,   0};
	      static const Real scal24[] =     {-0.25, 0,     0,     -0.25, -0.25, 0,     0,     -0.25, 0,     0,     0,     0.5,   0.5,   0,     0,     0.5,   0,     0,     0,     0.5};
	      static const Real scal25[] =     {0,     0,     0,     0,     -0.25, -0.25, -0.25, -0.25, 0,     0,     0,     0,     0,     0,     0,     0,     0.5,   0.5,   0.5,   0.5};
	      static const Real scal26[] =     {-0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, -0.25, 0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25,  0.25};
	      
	      return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta)
		      +scal20[i]*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[20], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[20], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[20], zeta)
		      +scal21[i]*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[21], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[21], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[21], zeta)
		      +scal22[i]*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[22], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[22], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[22], zeta)
		      +scal23[i]*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[23], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[23], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[23], zeta)
		      +scal24[i]*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[24], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[24], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[24], zeta)
		      +scal25[i]*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[25], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[25], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[25], zeta)
		      +scal26[i]*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[26], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[26], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[26], zeta));
	    }
	    
	    // Bernstein shape functions on the hexahedral.
	  case HEX27:
	    {
	      libmesh_assert (i<27);
	      
	      // Compute hex shape functions as a tensor-product
	      const Real xi   = p(0);
	      const Real eta  = p(1);
	      const Real zeta = 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!
	      //                                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
	      static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 0, 2, 2, 1, 2, 0, 2, 2};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 2, 1, 2, 0, 0, 1, 1, 0, 2, 1, 2, 2, 0, 2, 1, 2, 2, 2};
	      static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 2, 2, 1, 1, 1, 1, 0, 2, 2, 2, 2, 1, 2};
	      
	      return (FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i0[i], xi)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i1[i], eta)*
		      FE<1,BERNSTEIN>::shape(EDGE3, totalorder, i2[i], zeta));
	    }

	    
	  default:
	    libmesh_error();
	  }
	
      }

      
      
      // 3rd-order Bernstein.
    case THIRD:
      {
	switch (type)
	  {
	  
// 	    // Bernstein shape functions on the tetrahedron.
// 	  case TET10:
// 	    {
// 	      libmesh_assert(i<20);
	      
// 	      // Area coordinates
// 	      const Real zeta1 = p(0);
// 	      const Real zeta2 = p(1);
// 	      const Real zeta3 = p(2);
// 	      const Real zeta0 = 1. - zeta1 - zeta2 - zeta3;
	      
	      
// 	      unsigned int shape=i;
	      
// 	      // handle the edge orientation
	      
// 	      if ((i== 4||i== 5) && elem->node(0) > elem->node(1))shape= 9-i;   //Edge 0
// 	      if ((i== 6||i== 7) && elem->node(1) > elem->node(2))shape=13-i;   //Edge 1
// 	      if ((i== 8||i== 9) && elem->node(0) > elem->node(2))shape=17-i;   //Edge 2
// 	      if ((i==10||i==11) && elem->node(0) > elem->node(3))shape=21-i;   //Edge 3
// 	      if ((i==12||i==13) && elem->node(1) > elem->node(3))shape=25-i;   //Edge 4
// 	      if ((i==14||i==15) && elem->node(2) > elem->node(3))shape=29-i;   //Edge 5
	      
// 	      // No need to handle face orientation in 3rd order.
	      
	      
// 	      switch(shape)
// 		{
// 		  //point function
// 		case  0:  return zeta0*zeta0*zeta0;
// 		case  1:  return zeta1*zeta1*zeta1;
// 		case  2:  return zeta2*zeta2*zeta2;
// 		case  3:  return zeta3*zeta3*zeta3;
		  
// 		  //edge functions
// 		case  4:  return 3.*zeta0*zeta0*zeta1;
// 		case  5:  return 3.*zeta1*zeta1*zeta0;
			  
// 		case  6:  return 3.*zeta1*zeta1*zeta2;
// 		case  7:  return 3.*zeta2*zeta2*zeta1;
		  
// 		case  8:  return 3.*zeta0*zeta0*zeta2;
// 		case  9:  return 3.*zeta2*zeta2*zeta0;
		  
// 		case 10:  return 3.*zeta0*zeta0*zeta3;
// 		case 11:  return 3.*zeta3*zeta3*zeta0;
		  
// 		case 12:  return 3.*zeta1*zeta1*zeta3;
// 		case 13:  return 3.*zeta3*zeta3*zeta1;
		  
// 		case 14:  return 3.*zeta2*zeta2*zeta3;
// 		case 15:  return 3.*zeta3*zeta3*zeta2;
			  
// 		  //face functions
// 		case 16:  return 6.*zeta0*zeta1*zeta2;
// 		case 17:  return 6.*zeta0*zeta1*zeta3;
// 		case 18:  return 6.*zeta1*zeta2*zeta3;
// 		case 19:  return 6.*zeta2*zeta0*zeta3;
		  
// 		default:
// 		  libmesh_error();
// 		}
// 	    }
	    
	    
	    // Bernstein shape functions on the hexahedral.
	  case HEX27:
	    {
	      libmesh_assert (i<64);
	      
	      // 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
	      static const unsigned int i0[] = {0, 1, 1, 0, 0, 1, 1, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3, 2, 3};
	      static const unsigned int i1[] = {0, 0, 1, 1, 0, 0, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 2, 3, 1, 1, 2, 3, 2, 2, 3, 3, 0, 0, 0, 0, 2, 3, 2, 3, 1, 1, 1, 1, 2, 3, 2, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3};
	      static const unsigned int i2[] = {0, 0, 0, 0, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 2, 3, 2, 3, 2, 3, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 2, 2, 3, 3, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3};

	      

	      // 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 ))
		  {
		    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;
			}