text2.f90

fortran程序

              case(1)          ! for tetrahedra weights multiplied by 1/6
                 s(1,1)=.25    ; s(1,2)=.25  ;  s(1,3)=.25   ; wt(1)=1./6.
              case(4)
               s(1,1)=.58541020 ; s(1,2)=.13819660  ;  s(1,3)=s(1,2)
               s(2,2)=s(1,1) ; s(2,3)=s(1,2)  ;  s(2,1)=s(1,2)
               s(3,3)=s(1,1) ; s(3,1)=s(1,2)  ;  s(3,2)=s(1,2)
               s(4,1)=s(1,2) ; s(4,2)=s(1,2)  ;  s(4,3)=s(1,2) ; wt(1:4)=.25/6.
              case(5)
                s(1,1)=.25  ;  s(1,2)=.25   ; s(1,3)=.25 ;  s(2,1)=.5
                s(2,2)=1./6. ;  s(2,3)=s(2,2);  s(3,2)=.5
                s(3,3)=1./6.  ;   s(3,1)=s(3,3)   ;   s(4,3)=.5
                s(4,1)=1./6. ;    s(4,2)=s(4,1);    s(5,1)=1./6.
                s(5,2)=s(5,1) ;  s(5,3)=s(5,1) 
                wt(1)=-.8  ;  wt(2)=9./20. ;   wt(3:5)=wt(2)   ; wt =wt/6.  
              case(6)
         wt = 4./3.        ;  s(6,3) = 1.
         s(1,1)=-1. ;s(2,1)=1. ; s(3,2)=-1. ; s(4,2)=1. ;  s(5,3)=-1. 
              case default
               print*,"wrong number of integrating points for a tetrahedron"
            end select
            case('hexahedron')
             select case ( nip )
              case(1)
                     s(1,1) = .0 ; wt(1) = 8.    
              case(8)   
                     s(1,1)= root3;s(1,2)= root3;s(1,3)= root3
                     s(2,1)= root3;s(2,2)= root3;s(2,3)=-root3
                     s(3,1)= root3;s(3,2)=-root3;s(3,3)= root3
                     s(4,1)= root3;s(4,2)=-root3;s(4,3)=-root3
                     s(5,1)=-root3;s(5,2)= root3;s(5,3)= root3
                     s(6,1)=-root3;s(6,2)=-root3;s(6,3)= root3
                     s(7,1)=-root3;s(7,2)= root3;s(7,3)=-root3
                     s(8,1)=-root3;s(8,2)=-root3;s(8,3)=-root3
                     wt = 1.0                                                
              case(14)
          b=0.795822426     ;      c=0.758786911
          wt(1:6)=0.886426593   ; wt(7:) =  0.335180055
          s(1,1)=-b ; s(2,1)=b  ;  s(3,2)=-b ;   s(4,2)=b
          s(5,3)=-b   ;     s(6,3)=b
          s(7:,:) = c
          s(7,1)=-c  ;  s(7,2)=-c  ; s(7,3)=-c ; s(8,2)=-c ;   s(8,3)=-c
          s(9,1)=-c  ;  s(9,3)=-c  ; s(10,3)=-c; s(11,1)=-c
          s(11,2)=-c ;  s(12,2)=-c ; s(13,1)=-c
              case(15)
          b=1.     ;      c=0.674199862
          wt(1)=1.564444444 ;  wt(2:7)=0.355555556  ; wt(8:15)=0.537777778
          s(2,1)=-b  ;    s(3,1)=b  ;    s(4,2)=-b  ;    s(5,2)=b
          s(6,3)=-b  ;    s(7,3)=b  ;    s(8:,:)=c  ;    s(8,1)=-c
          s(8,2)=-c  ;    s(8,3)=-c ;    s(9,2)=-c  ;    s(9,3)=-c
          s(10,1)=-c ;    s(10,3)=-c  ;  s(11,3)=-c ;    s(12,1)=-c
          s(12,2)=-c ;    s(13,2)=-c  ;  s(14,1)=-c                          
              case(27)
                     wt = (/5./9.*v,8./9.*v,5./9.*v/)
                     s(1:7:3,1) = -r15; s(2:8:3,1) = .0
                     s(3:9:3,1) =  r15; s(1:3,3)   = r15
                     s(4:6,3)   =  .0 ; s(7:9,3)   =-r15
                     s(1:9,2)   = -r15
                     s(10:16:3,1) = -r15; s(11:17:3,1) = .0
                     s(12:18:3,1) =  r15; s(10:12,3)   = r15
                     s(13:15,3)   =  .0 ; s(16:18,3)   =-r15
                     s(10:18,2)   = .0   
                     s(19:25:3,1) = -r15; s(20:26:3,1) = .0
                     s(21:27:3,1) =  r15; s(19:21,3)   = r15
                     s(22:24,3)   =  .0 ; s(25:27,3)   =-r15
                     s(19:27,2)   =  r15
               case default
                 print*,"wrong number of integrating points for a hexahedron" 
             end select
            case default
             print*,"not a valid element type" 
     end select
   return
 end subroutine sample
 
 subroutine shape_der(der,points,i)
 implicit none
 integer,intent(in):: i; real,intent(in)::points(:,:)
 real,intent(out)::der(:,:)
  real::eta,xi,zeta,xi0,eta0,zeta0,etam,etap,xim,xip,c1,c2,c3 ! local variables
  real:: t1,t2,t3,t4,t5,t6,t7,t8,t9 ,x2p1,x2m1,e2p1,e2m1,zetam,zetap,x,y,z
  integer :: xii(20), etai(20), zetai(20) ,l,ndim , nod   ! local variables
  ndim = ubound(der , 1); nod = ubound(der , 2)
  select case (ndim)
   case(1) ! one dimensional case
         xi=points(i,1)
     select case (nod)
         case(2)
           der(1,1)=-0.5 ; der(1,2)=0.5
         case(3)
           t1=-1.-xi ; t2=-xi  ; t3=1.-xi
           der(1,1)=-(t3+t2)/2.  ; der(1,2)=(t3+t1)    
           der(1,3)=-(t2+t1)/2.   
         case(4)
           t1=-1.-xi ; t2=-1./3.-xi ; t3=1./3.-xi ; t4=1.-xi
           der(1,1)=-(t3*t4+t2*t4+t2*t3)*9./16.     
           der(1,2)=(t3*t4+t1*t4+t1*t3)*27./16. 
           der(1,3)=-(t2*t4+t1*t4+t1*t2)*27./16. 
           der(1,4)=(t2*t3+t1*t3+t1*t2)*9./16.   
         case(5)
           t1=-1.-xi ; t2=-0.5-xi ; t3=-xi ; t4=0.5-xi ; t5=1.-xi
           der(1,1)=-(t3*t4*t5+t2*t4*t5+t2*t3*t5+t2*t3*t4)*2./3.   
           der(1,2)=(t3*t4*t5+t1*t4*t5+t1*t3*t5+t1*t3*t4)*8./3.
           der(1,3)=-(t2*t4*t5+t1*t4*t5+t1*t2*t5+t1*t2*t4)*4. 
           der(1,4)=(t2*t3*t5+t1*t3*t5+t1*t2*t5+t1*t2*t3)*8./3.
           der(1,5)=-(t2*t3*t4+t1*t3*t4+t1*t2*t4+t1*t2*t3)*2./3.
     case default
       print*,"wrong number of nodes in shape_der"        
     end select
   case(2)      ! two dimensional elements
       xi=points(i,1); eta=points(i,2) ; c1=xi ; c2=eta ; c3=1.-c1-c2
       etam=.25*(1.-eta); etap=.25*(1.+eta); xim=.25*(1.-xi); xip=.25*(1.+xi)
       x2p1=2.*xi+1. ;   x2m1=2.*xi-1. ;  e2p1=2.*eta+1. ;   e2m1=2.*eta-1.
     select case (nod)
      case(3)
       der(1,1)=1.;der(1,3)=0.;der(1,2)=-1.
       der(2,1)=0.;der(2,3)=1.;der(2,2)=-1.
      case(6) 
       der(1,1)=4.*c1-1. ;  der(1,6)=4.*c2;  der(1,5)=0.  ; der(1,4)=-4.*c2
       der(1,3)=-(4.*c3-1.);  der(1,2)=4.*(c3-c1);   der(2,1)=0.
       der(2,6)=4.*c1 ; der(2,5)=4.*c2-1.; der(2,4)=4.*(c3-c2)
       der(2,3)=-(4.*c3-1.)  ; der(2,2)=-4.*c1
      case(15)                          
       t1=c1-.25  ;  t2=c1-.5 ;  t3=c1-.75   ;   t4=c2-.25
       t5=c2-.5   ;  t6=c2-.75 ;  t7=c3-.25  ;   t8=c3-.5 ;  t9=c3-.75
       der(1,1)=32./3.*(t2*t3*(t1+c1)+c1*t1*(t3+t2))
       der(1,12)=128./3.*c2*(t2*(t1+c1)+c1*t1) ;  der(1,11)=64.*c2*t4*(t1+c1)
       der(1,10)=128./3.*c2*t4*t5  ; der(1,9)=0. ; der(1,8)=-128./3.*c2*t4*t5
       der(1,7)=-64.*c2*t4*(t7+c3) ; der(1,6)=-128./3.*c2*(t8*(t7+c3)+c3*t7)
       der(1,5)=-32./3.*(t8*t9*(t7+c3)+c3*t7*(t8+t9))
       der(1,4)=128./3.*(c3*t7*t8-c1*(t8*(t7+c3)+c3*t7))
       der(1,3)=64.*(c3*t7*(t1+c1)-c1*t1*(t7+c3))
       der(1,2)=128./3.*(c3*(t2*(t1+c1)+c1*t1)-c1*t1*t2)
       der(1,13)=128.*c2*(c3*(t1+c1)-c1*t1) ;  der(1,15)=128.*c2*t4*(c3-c1)
       der(1,14)=128.*c2*(c3*t7-c1*(t7+c3))
       der(2,1)=0.0 ;  der(2,12)=128./3.*c1*t1*t2;  der(2,11)=64.*c1*t1*(t4+c2)
       der(2,10)=128./3.*c1*(t5*(t4+c2)+c2*t4)
       der(2,9)=32./3.*(t5*t6*(t4+c2)+c2*t4*(t6+t5))
       der(2,8)=128./3.*((c3*(t5*(t4+c2)+c2*t4))-c2*t4*t5)
       der(2,7)=64.*(c3*t7*(t4+c2)-c2*t4*(t7+c3))
       der(2,6)=128./3.*(c3*t7*t8-c2*(t8*(t7+c3)+c3*t7))
       der(2,5)=-32./3.*(t8*t9*(t7+c3)+c3*t7*(t8+t9))
       der(2,4)=-128./3.*c1*(t8*(t7+c3)+c3*t7)
       der(2,3)=-64.*c1*t1*(t7+c3)  ;  der(2,2)=-128./3.*c1*t1*t2
       der(2,13)=128.*c1*t1*(c3-c2)
       der(2,15)=128.*c1*(c3*(t4+c2)-c2*t4)
       der(2,14)=128.*c1*(c3*t7-c2*(c3+t7))        
      case (4)                                                              
       der(1,1)=-etam; der(1,2)=-etap; der(1,3)=etap; der(1,4)=etam
       der(2,1)=-xim; der(2,2)=xim; der(2,3)=xip; der(2,4)=-xip
      case(8)
       der(1,1)=etam*(2.*xi+eta); der(1,2)=-8.*etam*etap
       der(1,3)=etap*(2.*xi-eta); der(1,4)=-4.*etap*xi
       der(1,5)=etap*(2.*xi+eta); der(1,6)=8.*etap*etam
       der(1,7)=etam*(2.*xi-eta); der(1,8)=-4.*etam*xi
       der(2,1)=xim*(xi+2.*eta); der(2,2)=-4.*xim*eta
       der(2,3)=xim*(2.*eta-xi); der(2,4)=8.*xim*xip
       der(2,5)=xip*(xi+2.*eta); der(2,6)=-4.*xip*eta
       der(2,7)=xip*(2.*eta-xi); der(2,8)=-8.*xim*xip   
     case(9)
       etam = eta - 1.; etap = eta + 1.; xim = xi - 1.; xip = xi + 1.
       der(1,1)=.25*x2m1*eta*etam  ;   der(1,2)=-.5*x2m1*etap*etam
       der(1,3)=.25*x2m1*eta*etap  ;   der(1,4)=-xi*eta*etap
       der(1,5)=.25*x2p1*eta*etap  ;   der(1,6)=-.5*x2p1*etap*etam
       der(1,7)=.25*x2p1*eta*etam  ;   der(1,8)=-xi*eta*etam
       der(1,9)=2.*xi*etap*etam    ;   der(2,1)=.25*xi*xim*e2m1
       der(2,2)=-xi*xim*eta        ;   der(2,3)=.25*xi*xim*e2p1
       der(2,4)=-.5*xip*xim*e2p1   ;   der(2,5)=.25*xi*xip*e2p1
       der(2,6)=-xi*xip*eta        ;   der(2,7)=.25*xi*xip*e2m1
       der(2,8)=-.5*xip*xim*e2m1   ;   der(2,9)=2.*xip*xim*eta
     case default
       print*,"wrong number of nodes in shape_der"        
     end select
   case(3)  ! three dimensional elements
       xi=points(i,1); eta=points(i,2); zeta=points(i,3)
       etam=1.-eta ; xim=1.-xi;  zetam=1.-zeta
       etap=eta+1. ; xip=xi+1. ;  zetap=zeta+1.
    select case (nod)
     case(4)
      der(1:3,1:4) = .0
      der(1,1)=1.;  der(2,2)=1.  ;  der(3,3)=1.
      der(1,4)=-1. ;  der(2,4)=-1. ;  der(3,4)=-1.  
     case(8)
      der(1,1)=-.125*etam*zetam    ;   der(1,2)=-.125*etam*zetap
      der(1,3)=.125*etam*zetap     ;   der(1,4)=.125*etam*zetam
      der(1,5)=-.125*etap*zetam    ;   der(1,6)=-.125*etap*zetap
      der(1,7)=.125*etap*zetap     ;   der(1,8)=.125*etap*zetam
      der(2,1)=-.125*xim*zetam     ;   der(2,2)=-.125*xim*zetap
      der(2,3)=-.125*xip*zetap     ;   der(2,4)=-.125*xip*zetam
      der(2,5)=.125*xim*zetam      ;   der(2,6)=.125*xim*zetap
      der(2,7)=.125*xip*zetap      ;   der(2,8)=.125*xip*zetam
      der(3,1)=-.125*xim*etam      ;   der(3,2)=.125*xim*etam
      der(3,3)=.125*xip*etam       ;   der(3,4)=-.125*xip*etam
      der(3,5)=-.125*xim*etap      ;   der(3,6)=.125*xim*etap
      der(3,7)=.125*xip*etap       ;   der(3,8)=-.125*xip*etap  
 case(14) ! type 6 element
  x= points(i,1)    ;   y= points(i,2)  ;    z= points(i,3) 
  der(1,1)=((2.*x*y+2.*x*z+4.*x+y*z+y+z)*(y-1.)*(z-1.))/8.
  der(1,2)=((2.*x*y-2.*x*z-4.*x+y*z+y-z)*(y+1.)*(z-1.))/8.
  der(1,3)=((2.*x*y+2.*x*z+4.*x-y*z-y-z)*(y-1.)*(z-1.))/8.
  der(1,4)=((2.*x*y-2.*x*z-4.*x-y*z-y+z)*(y+1.)*(z-1.))/8.
  der(1,5)=-((2.*x*y-2.*x*z+4.*x-y*z+y-z)*(y-1.)*(z+1.))/8.
  der(1,6)=-((2.*x*y+2.*x*z-4.*x-y*z+y+z)*(y+1.)*(z+1.))/8.
  der(1,7)=-((2.*x*y-2.*x*z+4.*x+y*z-y+z)*(y-1.)*(z+1.))/8.
  der(1,8)=-((2.*x*y+2.*x*z-4.*x+y*z-y-z)*(y+1.)*(z+1.))/8.
  der(1,9)=-(y+1.)*(y-1.)*(z-1.)*x  ;   der(1,10)=(y+1.)*(y-1.)*(z+1.)*x
  der(1,11)=-(y-1.)*(z+1.)*(z-1.)*x ;   der(1,12)=(y+1.)*(z+1.)*(z-1.)*x
  der(1,13)=-((y+1.)*(y-1.)*(z+1.)*(z-1.))/2.
  der(1,14)=((y+1.)*(y-1.)*(z+1.)*(z-1.))/2.
  der(2,1)=((2.*x*y+x*z+x+2.*y*z+4.*y+z)*(x-1.)*(z-1.))/8.
  der(2,2)=((2.*x*y-x*z-x+2.*y*z+4.*y-z)*(x-1.)*(z-1.))/8.
  der(2,3)=((2.*x*y+x*z+x-2.*y*z-4.*y-z)*(x+1.)*(z-1.))/8.
  der(2,4)=((2.*x*y-x*z-x-2.*y*z-4.*y+z)*(x+1.)*(z-1.))/8.
  der(2,5)=-((2.*x*y-x*z+x-2.*y*z+4.*y-z)*(x-1.)*(z+1.))/8.
  der(2,6)=-((2.*x*y+x*z-x-2.*y*z+4.*y+z)*(x-1.)*(z+1.))/8.
  der(2,7)=-((2.*x*y-x*z+x+2.*y*z-4.*y+z)*(x+1.)*(z+1.))/8.
  der(2,8)=-((2.*x*y+x*z-x+2.*y*z-4.*y-z)*(x+1.)*(z+1.))/8.
  der(2,9)=-(x+1.)*(x-1.)*(z-1.)*y
  der(2,10)=(x+1.)*(x-1.)*(z+1.)*y
  der(2,11)=-((x+1.)*(x-1.)*(z+1.)*(z-1.))/2.
  der(2,12)=((x+1.)*(x-1.)*(z+1.)*(z-1.))/2.
  der(2,13)=-(x-1.)*(z+1.)*(z-1.)*y
  der(2,14)=(x+1.)*(z+1.)*(z-1.)*y
  der(3,1)=((x*y+2.*x*z+x+2.*y*z+y+4.*z)*(x-1.)*(y-1.))/8.
  der(3,2)=((x*y-2.*x*z-x+2.*y*z+y-4.*z)*(x-1.)*(y+1.))/8.
  der(3,3)=((x*y+2.*x*z+x-2.*y*z-y-4.*z)*(x+1.)*(y-1.))/8.
  der(3,4)=((x*y-2.*x*z-x-2.*y*z-y+4.*z)*(x+1.)*(y+1.))/8.
  der(3,5)=-((x*y-2.*x*z+x-2.*y*z+y-4.*z)*(x-1.)*(y-1.))/8.
  der(3,6)=-((x*y+2.*x*z-x-2.*y*z+y+4.*z)*(x-1.)*(y+1.))/8.
  der(3,7)=-((x*y-2.*x*z+x+2.*y*z-y+4.*z)*(x+1.)*(y-1.))/8.
  der(3,8)=-((x*y+2.*x*z-x+2.*y*z-y-4.*z)*(x+1.)*(y+1.))/8.
  der(3,9)=-((x+1.)*(x-1.)*(y+1.)*(y-1.))/2.
  der(3,10)=((x+1.)*(x-1.)*(y+1.)*(y-1.))/2.
  der(3,11)=-(x+1.)*(x-1.)*(y-1.)*z  ; der(3,12)=(x+1.)*(x-1.)*(y+1.)*z
  der(3,13)=-(x-1.)*(y+1.)*(y-1.)*z  ; der(3,14)=(x+1.)*(y+1.)*(y-1.)*z
    case(20)
      xii=(/-1,-1,-1,0,1,1,1,0,-1,-1,1,1,-1,-1,-1,0,1,1,1,0/)
      etai=(/-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,1,1,1,1,1,1,1,1/)
      zetai=(/-1,0,1,1,1,0,-1,-1,-1,1,1,-1,-1,0,1,1,1,0,-1,-1/)
      do l=1,20
         xi0=xi*xii(l); eta0=eta*etai(l); zeta0=zeta*zetai(l)
         if(l==4.or.l==8.or.l==16.or.l==20) then
          der(1,l)=-.5*xi*(1.+eta0)*(1.+zeta0)
          der(2,l)=.25*etai(l)*(1.-xi*xi)*(1.+zeta0)
          der(3,l)=.25*zetai(l)*(1.-xi*xi)*(1.+eta0)
         else if(l>=9.and.l<=12)then
          der(1,l)=.25*xii(l)*(1.-eta*eta)*(1.+zeta0)
          der(2,l)=-.5*eta*(1.+xi0)*(1.+zeta0)
          der(3,l)=.25*zetai(l)*(1.+xi0)*(1.-eta*eta)
         else if(l==2.or.l==6.or.l==14.or.l==18) then
          der(1,l)=.25*xii(l)*(1.+eta0)*(1.-zeta*zeta)
          der(2,l)=.25*etai(l)*(1.+xi0)*(1.-zeta*zeta)
          der(3,l)=-.5*zeta*(1.+xi0)*(1.+eta0)
         else
          der(1,l)=.125*xii(l)*(1.+eta0)*(1.+zeta0)*(2.*xi0+eta0+zeta0-1.)
          der(2,l)=.125*etai(l)*(1.+xi0)*(1.+zeta0)*(xi0+2.*eta0+zeta0-1.)
          der(3,l)=.125*zetai(l)*(1.+xi0)*(1.+eta0)*(xi0+eta0+2.*zeta0-1.)
         end if
      end do 
     case default
       print*,"wrong number of nodes in shape_der"        
   end select
  case default
   print*,"wrong number of dimensions in shape_der"
  end select
 return
 end subroutine shape_der
 
 subroutine shape_fun(fun,points,i)
  implicit none  
  integer,intent(in):: i; real,intent(in)::points(:,:)
  real,intent(out)::fun(:)
  real :: eta,xi,etam,etap,xim,xip,zetam,zetap,c1,c2,c3     !local variables
  real :: t1,t2,t3,t4,t5,t6,t7,t8,t9,x,y,z
  real :: zeta,xi0,eta0,zeta0; integer::xii(20),etai(20),zetai(20),l,ndim,nod
        ndim = ubound(points , 2 ); nod = ubound(fun , 1 )  
    select case (ndim)
      case(1) ! one dimensional cases
           xi=points(i,1)
        select case(nod)
         case(2)
           t1=-1.-xi ; t2=1.-xi
           fun(1)=t2/2. ; fun(2)=-t1/2.
         case(3)
           t1=-1.-xi ; t2=-xi ; t3=1.-xi
           fun(1)=t2*t3/2. ; fun(2)=-t1*t3 ; fun(3)=t1*t2/2.
         case(4)
           t1=-1.-xi ; t2=-1./3.-xi ; t3=1./3.-xi ; t4=1.-xi
           fun(1)=t2*t3*t4*9./16.  ; fun(2)=-t1*t3*t4*27./16.
           fun(3)=t1*t2*t4*27./16. ; fun(4)=-t1*t2*t3*9./16.
         case(5)
           t1=-1.-xi ; t2=-0.5-xi ; t3=-xi ; t4=0.5-xi ; t5=1.-xi
           fun(1)=t2*t3*t4*t5*2./3. ; fun(2)=-t1*t3*t4*t5*8./3.
           fun(3)=t1*t2*t4*t5*4. ; fun(4)=-t1*t2*t3*t5*8./3.
           fun(5)=t1*t2*t3*t4*2./3.
          case default
             print*,"wrong number of nodes in shape_fun"
        end select
      case(2) ! two dimensional cases
           c1=points(i,1); c2=points(i,2); c3=1.-c1-c2 
           xi=points(i,1);  eta=points(i,2)
           etam=.25*(1.-eta); etap=.25*(1.+eta)
           xim=.25*(1.-xi); xip=.25*(1.+xi)
        select case(nod)
          case(3)
            fun = (/c1,c3,c2/)  
          case(6)
            fun(1)=(2.*c1-1.)*c1 ;  fun(6)=4.*c1*c2 ;  fun(5)=(2.*c2-1.)*c2
            fun(4)=4.*c2*c3      ;  fun(3)=(2.*c3-1.)*c3 ; fun(2)=4.*c3*c1
          case(15)
            t1=c1-.25  ;  t2=c1-.5 ;  t3=c1-.75   ;   t4=c2-.25
            t5=c2-.5   ;  t6=c2-.75 ;  t7=c3-.25  ;   t8=c3-.5 ;  t9=c3-.75
            fun(1)=32./3.*c1*t1*t2*t3   ;  fun(12)=128./3.*c1*c2*t1*t2
            fun(11)=64.*c1*c2*t1*t4     ;  fun(10)=128./3.*c1*c2*t4*t5
            fun(9)=32./3.*c2*t4*t5*t6   ;  fun(8)=128./3.*c2*c3*t4*t5
            fun(7)=64.*c2*c3*t4*t7      ;  fun(6)=128./3.*c2*c3*t7*t8
            fun(5)=32./3.*c3*t7*t8*t9   ;  fun(4)=128./3.*c3*c1*t7*t8
            fun(3)=64.*c3*c1*t1*t7      ;  fun(2)=128./3.*c3*c1*t1*t2
            fun(13)=128.*c1*c2*t1*c3    ;  fun(15)=128.*c1*c2*c3*t4
            fun(14)=128.*c1*c2*c3*t7      
          case(4)
            fun=(/4.*xim*etam,4.*xim*etap,4.*xip*etap,4.*xip*etam/)
          case(8)
            fun=(/4.*etam*xim*(-xi-eta-1.),32.*etam*xim*etap,&
                  4.*etap*xim*(-xi+eta-1.),32.*xim*xip*etap, &
                  4.*etap*xip*(xi+eta-1.), 32.*etap*xip*etam,&
                  4.*xip*etam*(xi-eta-1.), 32.*xim*xip*etam/)
          case(9)
            etam = eta - 1.; etap= eta + 1.; xim = xi - 1.; xip = xi + 1.
            fun=(/.25*xi*xim*eta*etam,-.5*xi*xim*etap*etam,&
                  .25*xi*xim*eta*etap,-.5*xip*xim*eta*etap,&
                  .25*xi*xip*eta*etap,-.5*xi*xip*etap*etam,&
                  .25*xi*xip*eta*etam,-.5*xip*xim*eta*etam,xip*xim*etap*etam/)
          case default
             print*,"wrong number of nodes in shape_fun"
        end select