shape20h_ax.f File Reference

Go to the source code of this file.

Functions/Subroutines
subroutine	shape20h_ax (xi, et, ze, xl, xsj, shp, iflag)

Function/Subroutine Documentation

shape20h_ax()

subroutine shape20h_ax	(	real*8, intent(in)	xi,
		real*8, intent(in)	et,
		real*8, intent(in)	ze,
		real*8, dimension(3,20), intent(in)	xl,
		real*8, intent(out)	xsj,
		real*8, dimension(4,20), intent(out)	shp,
		integer, intent(in)	iflag
	)

!
!     shape functions and derivatives for a 20-node quadratic
!     isoparametric brick element. -1<=xi,et,ze<=1
!     special case: axisymmetric elements
!
!     iflag=1: calculate only the value of the shape functions
!     iflag=2: calculate the value of the shape functions and
!              the Jacobian determinant
!     iflag=3: calculate the value of the shape functions, the
!              value of their derivatives w.r.t. the global
!              coordinates and the Jacobian determinant
!
      implicit none
!
      integer j,k,iflag
!
      real*8 shp(4,20),xs(3,3),xsi(3,3),xl(3,20),shpe(4,20),dd,
     &  dd1,dd2,dd3
!
      real*8 xi,et,ze,xsj,omg,omh,omr,opg,oph,opr,
     &  tpgphpr,tmgphpr,tmgmhpr,tpgmhpr,tpgphmr,tmgphmr,tmgmhmr,tpgmhmr,
     &  omgopg,omhoph,omropr,omgmopg,omhmoph,omrmopr
!
      intent(in) xi,et,ze,xl,iflag
!
      intent(out) shp,xsj
!
!     shape functions and their glocal derivatives
!
      omg=1.d0-xi
      omh=1.d0-et
      omr=1.d0-ze
      opg=1.d0+xi
      oph=1.d0+et
      opr=1.d0+ze
      tpgphpr=opg+oph+ze
      tmgphpr=omg+oph+ze
      tmgmhpr=omg+omh+ze
      tpgmhpr=opg+omh+ze
      tpgphmr=opg+oph-ze
      tmgphmr=omg+oph-ze
      tmgmhmr=omg+omh-ze
      tpgmhmr=opg+omh-ze
      omgopg=omg*opg/4.d0
      omhoph=omh*oph/4.d0
      omropr=omr*opr/4.d0
      omgmopg=(omg-opg)/4.d0
      omhmoph=(omh-oph)/4.d0
      omrmopr=(omr-opr)/4.d0
!
!     shape functions
!
      shp(4, 1)=-omg*omh*omr*tpgphpr/8.d0
      shp(4, 2)=-opg*omh*omr*tmgphpr/8.d0
      shp(4, 3)=-opg*oph*omr*tmgmhpr/8.d0
      shp(4, 4)=-omg*oph*omr*tpgmhpr/8.d0
      shp(4, 5)=-omg*omh*opr*tpgphmr/8.d0
      shp(4, 6)=-opg*omh*opr*tmgphmr/8.d0
      shp(4, 7)=-opg*oph*opr*tmgmhmr/8.d0
      shp(4, 8)=-omg*oph*opr*tpgmhmr/8.d0
      shp(4, 9)=omgopg*omh*omr
      shp(4,10)=omhoph*opg*omr
      shp(4,11)=omgopg*oph*omr
      shp(4,12)=omhoph*omg*omr
      shp(4,13)=omgopg*omh*opr
      shp(4,14)=omhoph*opg*opr
      shp(4,15)=omgopg*oph*opr
      shp(4,16)=omhoph*omg*opr
      shp(4,17)=omropr*omg*omh
      shp(4,18)=omropr*opg*omh
      shp(4,19)=omropr*opg*oph
      shp(4,20)=omropr*omg*oph
!
      if(iflag.eq.1) return
!
!     local derivatives of the shape functions: xi-derivative
!
      shpe(1, 1)=omh*omr*(tpgphpr-omg)/8.d0
      shpe(1, 2)=(opg-tmgphpr)*omh*omr/8.d0
      shpe(1, 3)=(opg-tmgmhpr)*oph*omr/8.d0
      shpe(1, 4)=oph*omr*(tpgmhpr-omg)/8.d0
      shpe(1, 5)=omh*opr*(tpgphmr-omg)/8.d0
      shpe(1, 6)=(opg-tmgphmr)*omh*opr/8.d0
      shpe(1, 7)=(opg-tmgmhmr)*oph*opr/8.d0
      shpe(1, 8)=oph*opr*(tpgmhmr-omg)/8.d0
      shpe(1, 9)=omgmopg*omh*omr
      shpe(1,10)=omhoph*omr
      shpe(1,11)=omgmopg*oph*omr
      shpe(1,12)=-omhoph*omr
      shpe(1,13)=omgmopg*omh*opr
      shpe(1,14)=omhoph*opr
      shpe(1,15)=omgmopg*oph*opr
      shpe(1,16)=-omhoph*opr
      shpe(1,17)=-omropr*omh
      shpe(1,18)=omropr*omh
      shpe(1,19)=omropr*oph
      shpe(1,20)=-omropr*oph
!
!     local derivatives of the shape functions: eta-derivative
!
      shpe(2, 1)=omg*omr*(tpgphpr-omh)/8.d0
      shpe(2, 2)=opg*omr*(tmgphpr-omh)/8.d0
      shpe(2, 3)=opg*(oph-tmgmhpr)*omr/8.d0
      shpe(2, 4)=omg*(oph-tpgmhpr)*omr/8.d0
      shpe(2, 5)=omg*opr*(tpgphmr-omh)/8.d0
      shpe(2, 6)=opg*opr*(tmgphmr-omh)/8.d0
      shpe(2, 7)=opg*(oph-tmgmhmr)*opr/8.d0
      shpe(2, 8)=omg*(oph-tpgmhmr)*opr/8.d0
      shpe(2, 9)=-omgopg*omr
      shpe(2,10)=omhmoph*opg*omr
      shpe(2,11)=omgopg*omr
      shpe(2,12)=omhmoph*omg*omr
      shpe(2,13)=-omgopg*opr
      shpe(2,14)=omhmoph*opg*opr
      shpe(2,15)=omgopg*opr
      shpe(2,16)=omhmoph*omg*opr
      shpe(2,17)=-omropr*omg
      shpe(2,18)=-omropr*opg
      shpe(2,19)=omropr*opg
      shpe(2,20)=omropr*omg
!
!     local derivatives of the shape functions: zeta-derivative
!
      shpe(3, 1)=omg*omh*(tpgphpr-omr)/8.d0
      shpe(3, 2)=opg*omh*(tmgphpr-omr)/8.d0
      shpe(3, 3)=opg*oph*(tmgmhpr-omr)/8.d0
      shpe(3, 4)=omg*oph*(tpgmhpr-omr)/8.d0
      shpe(3, 5)=omg*omh*(opr-tpgphmr)/8.d0
      shpe(3, 6)=opg*omh*(opr-tmgphmr)/8.d0
      shpe(3, 7)=opg*oph*(opr-tmgmhmr)/8.d0
      shpe(3, 8)=omg*oph*(opr-tpgmhmr)/8.d0
      shpe(3, 9)=-omgopg*omh
      shpe(3,10)=-omhoph*opg
      shpe(3,11)=-omgopg*oph
      shpe(3,12)=-omhoph*omg
      shpe(3,13)=omgopg*omh
      shpe(3,14)=omhoph*opg
      shpe(3,15)=omgopg*oph
      shpe(3,16)=omhoph*omg
      shpe(3,17)=omrmopr*omg*omh
      shpe(3,18)=omrmopr*opg*omh
      shpe(3,19)=omrmopr*opg*oph
      shpe(3,20)=omrmopr*omg*oph
!
!     computation of the local derivative of the global coordinates
!     (xs)
!
c      do i=1,3
c        do j=1,3
c          xs(i,j)=0.d0
c          do k=1,20
c            xs(i,j)=xs(i,j)+xl(i,k)*shpe(j,k)
c          enddo
c        enddo
c      enddo
      do j=1,3
         xs(1,j)=xl(1,1)*(shpe(j,1)+shpe(j,5))
     &          +xl(1,2)*(shpe(j,2)+shpe(j,6))
     &          +xl(1,3)*(shpe(j,3)+shpe(j,7))
     &          +xl(1,4)*(shpe(j,4)+shpe(j,8))
     &          +xl(1,9)*(shpe(j,9)+shpe(j,13))
     &          +xl(1,10)*(shpe(j,10)+shpe(j,14))
     &          +xl(1,11)*(shpe(j,11)+shpe(j,15))
     &          +xl(1,12)*(shpe(j,12)+shpe(j,16))
     &          +xl(1,17)*shpe(j,17)+xl(1,18)*shpe(j,18)
     &          +xl(1,19)*shpe(j,19)+xl(1,20)*shpe(j,20)
         xs(2,j)=xl(2,1)*(shpe(j,1)+shpe(j,5)+shpe(j,17))
     &          +xl(2,2)*(shpe(j,2)+shpe(j,6)+shpe(j,18))
     &          +xl(2,3)*(shpe(j,3)+shpe(j,7)+shpe(j,19))
     &          +xl(2,4)*(shpe(j,4)+shpe(j,8)+shpe(j,20))
     &          +xl(2,9)*(shpe(j,9)+shpe(j,13))
     &          +xl(2,10)*(shpe(j,10)+shpe(j,14))
     &          +xl(2,11)*(shpe(j,11)+shpe(j,15))
     &          +xl(2,12)*(shpe(j,12)+shpe(j,16))
         xs(3,j)=xl(3,1)*(shpe(j,1)-shpe(j,5))
     &          +xl(3,2)*(shpe(j,2)-shpe(j,6))
     &          +xl(3,3)*(shpe(j,3)-shpe(j,7))
     &          +xl(3,4)*(shpe(j,4)-shpe(j,8))
     &          +xl(3,9)*(shpe(j,9)-shpe(j,13))
     &          +xl(3,10)*(shpe(j,10)-shpe(j,14))
     &          +xl(3,11)*(shpe(j,11)-shpe(j,15))
     &          +xl(3,12)*(shpe(j,12)-shpe(j,16))
      enddo
!
!     computation of the jacobian determinant
!
      dd1=xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2)
      dd2=xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3)
      dd3=xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1)
      xsj=xs(1,1)*dd1+xs(1,2)*dd2+xs(1,3)*dd3
c      xsj=xs(1,1)*(xs(2,2)*xs(3,3)-xs(2,3)*xs(3,2))
c     &   -xs(1,2)*(xs(2,1)*xs(3,3)-xs(2,3)*xs(3,1))
c     &   +xs(1,3)*(xs(2,1)*xs(3,2)-xs(2,2)*xs(3,1))
!
      if(iflag.eq.2) return
!
      dd=1.d0/xsj
!
!     computation of the global derivative of the local coordinates
!     (xsi) (inversion of xs)
!
      xsi(1,1)=dd1*dd
      xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))*dd
      xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))*dd
      xsi(2,1)=dd2*dd
      xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))*dd
      xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))*dd
      xsi(3,1)=dd3*dd
      xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))*dd
      xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))*dd
c      xsi(1,1)=(xs(2,2)*xs(3,3)-xs(3,2)*xs(2,3))*dd
c      xsi(1,2)=(xs(1,3)*xs(3,2)-xs(1,2)*xs(3,3))*dd
c      xsi(1,3)=(xs(1,2)*xs(2,3)-xs(2,2)*xs(1,3))*dd
c      xsi(2,1)=(xs(2,3)*xs(3,1)-xs(2,1)*xs(3,3))*dd
c      xsi(2,2)=(xs(1,1)*xs(3,3)-xs(3,1)*xs(1,3))*dd
c      xsi(2,3)=(xs(1,3)*xs(2,1)-xs(1,1)*xs(2,3))*dd
c      xsi(3,1)=(xs(2,1)*xs(3,2)-xs(3,1)*xs(2,2))*dd
c      xsi(3,2)=(xs(1,2)*xs(3,1)-xs(1,1)*xs(3,2))*dd
c      xsi(3,3)=(xs(1,1)*xs(2,2)-xs(2,1)*xs(1,2))*dd
!
!     computation of the global derivatives of the shape functions
!
      do k=1,20
        do j=1,3
          shp(j,k)=shpe(1,k)*xsi(1,j)+shpe(2,k)*xsi(2,j)
     &          +shpe(3,k)*xsi(3,j)
        enddo
      enddo
!
      return