3DRWALK/source/trinvs.for

c----------------------------------------------------------------trinvs3
       subroutine trinvs3(nel,cno,itype,ebox,egrp,nblk,egps,nelb,nbb,
     *                        nee,nxyz,nen,it,x,y,z,xi,eta,zeta,num,nn)
c----------------------------------------------------------------------c
c purpose:                                                             c
c         To locate the element the coordinate belongs to and to       c
c         calculate the shape functions.                               c
c----------------------------------------------------------------------c
ccc    implicit real*8 (a-h,o-z)
       common /etype/inode(4)
       common /shape/shap(20),shpx(20),shpy(20),shpz(20)
       integer*2 nel(20,1)
       dimension pos(3),cno(3,1),nbb(1),nee(1),itype(1),ebox(6,1),
     *           nelb(1),egps(6,1),nblk(2,1),egrp(6,1)
       real shap,xi,eta,zeta
c
       pos(1)=x
       pos(2)=y
       pos(3)=z
       nn=0
       if(num.gt.0) then
c
c----------box test
c
           if(pos(1).lt.ebox(1,num).or.pos(1).gt.ebox(4,num)) goto 1
           if(pos(2).lt.ebox(2,num).or.pos(2).gt.ebox(5,num)) goto 1
           if(pos(3).lt.ebox(3,num).or.pos(3).gt.ebox(6,num)) goto 1
           it=itype(num)
           nen=inode(it)
           call invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii)
           if(key.eq.1) goto 30
       endif
c      call trinvs1(nel,cno,egrp,nblk,egps,nelb,nbb,nee,nxyz,nen,
c    *                                     it,x,y,z,xi,eta,zeta,num,nn)
c      if(nn.eq.1) goto 30
c      goto 35
    1  do n=1,nxyz
          if(pos(1).gt.egrp(4,n).or.pos(1).lt.egrp(1,n)) goto 10
          if(pos(2).gt.egrp(5,n).or.pos(2).lt.egrp(2,n)) goto 10
          if(pos(3).gt.egrp(6,n).or.pos(3).lt.egrp(3,n)) goto 10
          if(nblk(1,n).gt.1) then
            do j=nblk(2,n),nblk(2,n)+nblk(1,n)-1
               if(pos(1).gt.egps(4,j).or.pos(1).lt.egps(1,j)) goto 5
               if(pos(2).gt.egps(5,j).or.pos(2).lt.egps(2,j)) goto 5
               if(pos(3).gt.egps(6,j).or.pos(3).lt.egps(3,j)) goto 5
               ng=j
               goto 20
    5       enddo
          else
            ng=nblk(2,n)
            goto 20
          endif
   10  enddo
c
c------out of the finite element mesh
c
       goto 900

   20  do 35 i=nbb(ng),nee(ng)
           num=nelb(i)
c
c----------box test
c
           if(pos(1).lt.ebox(1,num).or.pos(1).gt.ebox(4,num)) goto 35
           if(pos(2).lt.ebox(2,num).or.pos(2).gt.ebox(5,num)) goto 35
           if(pos(3).lt.ebox(3,num).or.pos(3).gt.ebox(6,num)) goto 35
           it=itype(num)
           nen=inode(it)
           call invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii)
           if(key.eq.1) goto 30
   35  continue
c      write(*,'(/a/)') 'The input coordinate is out of the FE mesh.'
       goto 900
   30  nn=1
  900  return
 1000  format(2i10)
 1100  format(i10,3f10.0)
 1200  format(21i5)
 1300  format(i6,6e12.4,i6)
       end
c--------------------------------------------------------------invs_el3
       subroutine invs_el3(nel,cno,pos,xi,eta,zeta,nen,it,num,key,ii)
c---------------------------------------------------------------------c
c Purpose:                                                            c
c         To perform inverse mapping for element i.                   c
c---------------------------------------------------------------------c
       common /shape/shap(20),shpx(20),shpy(20),shpz(20)
       integer*2 nel(20,1)
       dimension nele(20),cnod(3,20),pos(3),cno(3,1)
       real shap,xi,eta,zeta
       real a(3,3),b(3),a1,a2,a3
       real*8 detj,ass
c
       do 40 j=1,nen
             nele(j)=nel(j,num)
             cnod(1,j)=cno(1,nele(j))
             cnod(2,j)=cno(2,nele(j))
             cnod(3,j)=cno(3,nele(j))
   40  continue
c
c------Newton-Raphson iteration method
c
       ii=0
       xi0=0.0
       eta0=0.0
       zeta0=0.0
c
c------calculate the values of shape funcations and their derivatives
c
    1  ii=ii+1
       call xn3x(it,nen,xi0,eta0,zeta0)
c
c------calculate the Jacobian matrix
c
       do 20 j=1,3
       a(j,1)=0.0
       a(j,2)=0.0
       a(j,3)=0.0
       b(j)=-pos(j)
       do 20 i=1,nen
             a(j,1)=a(j,1)+shpx(i)*cnod(j,i)
             a(j,2)=a(j,2)+shpy(i)*cnod(j,i)
             a(j,3)=a(j,3)+shpz(i)*cnod(j,i)
             b(j)=b(j)+shap(i)*cnod(j,i)
   20  continue

       if(abs(b(1)).lt.1.0e-5.and.abs(b(2)).lt.1.0e-5.and.
     &                            abs(b(3)).lt.1.0e-5) goto 30

       do i=1,3
          ass=max(abs(a(i,1)),abs(a(i,2)),abs(a(i,3)))
          if(dabs(ass).lt.1.0e-10) goto 25
          do j=1,3
             a(i,j)=a(i,j)/ass
          enddo
          b(i)=b(i)/ass
   25  enddo
       a1=a(2,2)*a(3,3)-a(2,3)*a(3,2)
       a2=a(3,2)*a(1,3)-a(3,3)*a(1,2)
       a3=a(1,2)*a(2,3)-a(1,3)*a(2,2)
       detj=a(1,1)*a1+a(2,1)*a2+a(3,1)*a3
       if(dabs(detj).lt.1.0e-10) goto 900
       xi=xi0-(a1*b(1)+a2*b(2)+a3*b(3))/detj
       a1=a(3,1)*a(2,3)-a(3,3)*a(2,1)
       a2=a(1,1)*a(3,3)-a(1,3)*a(3,1)
       a3=a(2,1)*a(1,3)-a(2,3)*a(1,1)
       eta=eta0-(a1*b(1)+a2*b(2)+a3*b(3))/detj
       a1=a(2,1)*a(3,2)-a(2,2)*a(3,1)
       a2=a(3,1)*a(1,2)-a(3,2)*a(1,1)
       a3=a(1,1)*a(2,2)-a(1,2)*a(2,1)
       zeta=zeta0-(a1*b(1)+a2*b(2)+a3*b(3))/detj
       if(abs(eta-eta0).le.0.5e-4.and.abs(xi-xi0).le.0.5e-4.and.
     *    abs(zeta-zeta0).le.0.5e-4) goto 30
       if(ii.gt.11) goto 30
       xi0=xi
       eta0=eta
       zeta0=zeta
       goto 1
c 900  write(6,*) 'Jacobian Determinant = 0',num,nen,it,xi0,eta0,zeta0,
c    *             detj
c      stop
  900  goto 35
c      goto 35
c
c------determine whether it is inside the element
c
   30  goto (50,60,70,80), it
   50  if(xi.lt.-5.0e-3.or.eta.lt.-5.0e-3.or.zeta.lt.-5.0e-3)
     *                                                goto 35
       if((xi+eta+zeta).gt.1.005) goto 35
       goto 90
   60  if(xi.lt.-5.0e-3.or.xi.gt.1.005) goto 35
       if(eta.lt.-1.005.or.eta.gt.1.005) goto 35
       if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35
       if((xi+abs(eta)).gt.1.005) goto 35
       goto 90
   70  if(xi.lt.-1.005.or.xi.gt.1.005) goto 35
       if(eta.lt.-1.005.or.eta.gt.1.005) goto 35
       if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35
       goto 90
   80  if(xi.lt.-5.0e-3.or.xi.gt.1.005) goto 35
       if(eta.lt.-1.005.or.eta.gt.1.005) goto 35
       if(zeta.lt.-1.005.or.zeta.gt.1.005) goto 35
       if((xi+abs(eta)).gt.1.005) goto 35
       if((xi+abs(zeta)).gt.1.005) goto 35
   90  key=1
       return
   35  key=0
       return
       end

c-------------------------------------------------------------------xn3
      subroutine xn3(it,nen,x,y,z)
c---------------------------------------------------------------------c
c
c
c     subroutine to evaluate shape function
c
c      IT is element type 1 10 point tetrahedron
c                    type 2 15 point prism
c                    type 3 20 point parallelipiped
c                    type 4 13 point rectangular base pyramid
c      nen  is number of shape functions
c      X, Y, Z  are cordinates of point to be evaluated in local coord
c
c---------------------------------------------------------------------c
ccc   implicit real*8 (a-h,o-z)
      save
      common /shape/shap(20),shpx(20),shpy(20),shpz(20)
      dimension irf(10,2)
      dimension a(5),b(5),c(5),xt(3),yt(3),jx(3),kx(3),xm(5)
     1         ,sn(2),shpp(15,4),ilokup(13)
c
      data ilokup/5,6,1,2,3,4,9,7,14,15,10,11,13/
      data irf/1,1,2,2,3,3,1,2,3,4
     1        ,1,2,2,3,3,1,4,4,4,4/
      data xt/0.0,1.0,0.0/,yt/-1.0,0.0,1.0/,jx/2,3,1/,kx/3,1,2/
      data shpp/0.0,0.0,1.0,12*0.0,
     +      0.25,-1.0,3.0,-1.0,0.25,4*0.0,0.25,-1.0,0.0,-1.0,0.25,0.0,
     +      0.25,-1.0,0.0,1.0,-0.25,4*0.0,0.25,-1.0,0.0,1.0,-0.25,0.0,
     +      0.25,-1.0,0.0,-1.0,0.25,4*0.0,-0.25,1.0,0.0,1.0,-0.25,0.0/
      data ncall/0/
c
      if(it .eq. 4) go to 80
      if( it - 2 ) 500,80,300
c-
c-----shape functions for right prism and pyramid.....
c-
   80 do j0=1,nen
         if(it.eq.2) then
            i=j0
         else
            i=ilokup(j0)
         endif
         ncall = ncall + 1
         if( ncall .gt. 1 ) go to 125
c-
c-----calculate invarient triangular functions.....
c-
         n = 0
         do 100 j = 1,5,2
         n = n + 1
         jj = jx(n)
         kk = kx(n)
         a(j) = xt(jj)*yt(kk) - xt(kk)*yt(jj)
         b(j) = yt(jj) - yt(kk)
         c(j) = xt(kk) - xt(jj)
  100    continue
c-
c-.....shape function calculations.....
c-
  125    do 130 k = 1, 5, 2
            xm(k) = 0.5*(a(k)+b(k)*x+c(k)*y)
  130    continue
c-
c-.....set indexes and normalize coordinates.....
c-
         if( i .le. 6 ) then
            zooo = -z
         else
            zooo = z
         endif
         if(it .eq. 2) go to 132
         IF(X .EQ. 1.0) GO TO 250
c        if(x .eq. 1.0) x=0.999
         if(abs(1.0-x) .le. 1.0e-6) then
c          print *,'x=',x
           goto 250
         endif
         zooo=zooo/(1.-x)
c        if(x .eq. 1.) zooo=sign(x,zooo)
         if(i .eq. 3) go to 200
  132    if( i .gt. 6 .and. i .lt. 10 ) go to 170
         l = i
         if( i .gt. 9 ) l = i - 9
         if( mod(l,2) .eq. 0 ) go to 150
c-
c-.....corner nodes.....
c-
         shap(j0)=xm(l)*((2.0*xm(l)-1.0)*(1.0+zooo)-(1.0-zooo**2))/2.0
         if(it .eq. 4) shap(j0)=shap(j0)+xm(l)/2.0*x*(1.-zooo**2)
         goto 240
c-
c-.....mid side node.....
c-
  150    n1 = l - 1
         n2 = mod(l+1,6)
         shap(j0) = 2.0*xm(n1)*xm(n2)*(1.0+zooo)
         goto 240
c-
c-.....mid side rectangle.....
c-
  170    n1 = 1
         if( i .eq. 8 ) n1 = 3
         if( i .eq. 9 ) n1 = 5
         shap(j0) = xm(n1)*(1.0-zooo**2)
         if(it .eq. 4) shap(j0)=shap(j0)*(1.-x)
  180    goto 240
c-
c......special case for no 3 shape function of pyramid
c-
  200    l=i
         shap(j0)=xm(l)*(2.*xm(l)-1.0)
         goto 240
c
c-----Special case when x exactly equals 1.0
c
  250    continue
         shap(j0)=shpp(i,1)
  240 enddo
      return
c-
c-.....shape functions for cube.....
c-
  300 zm1=0.125*(1.0-z)
      zp1=0.125*(1.0+z)
c-
c......corner nodes
c-
      shap(1)=(1.-x)*(1.-y)*zm1*(-x-y-z-2.)
      shap(3)=(1.+x)*(1.-y)*zm1*(x-y-z-2.)
      shap(5)=(1.+x)*(1.+y)*zm1*(x+y-z-2.)
      shap(7)=(1.-x)*(1.+y)*zm1*(-x+y-z-2.)
      shap(13)=(1.-x)*(1.-y)*zp1*(-x-y+z-2.)
      shap(15)=(1.+x)*(1.-y)*zp1*(x-y+z-2.)
      shap(17)=(1.+x)*(1.+y)*zp1*(x+y+z-2.)
      shap(19)=(1.-x)*(1.+y)*zp1*(-x+y+z-2.)
c-
c......mid-side nodes
c-
      x21=0.25*(1.0-x*x)
      y21=0.25*(1.0-y*y)
      z21=0.25*(1.0-z*z)
      shap(2)=x21*(1.-y)*(1.-z)
      shap(4)=(1.+x)*y21*(1.-z)
      shap(6)=x21*(1.+y)*(1.-z)
      shap(8)=(1.-x)*y21*(1.-z)
      shap(9)=(1.-x)*(1.-y)*z21
      shap(10)=(1.+x)*(1.-y)*z21
      shap(11)=(1.+x)*(1.+y)*z21
      shap(12)=(1.-x)*(1.+y)*z21
      shap(14)=x21*(1.-y)*(1.+z)
      shap(16)=(1.+x)*y21*(1.+z)
      shap(18)=x21*(1.+y)*(1.+z)
      shap(20)=(1.-x)*y21*(1.+z)
      return
  500 continue
c-
c......section for tetrahedron
c-
      do i=1,nen
         i1=irf(i,1)
         i2=irf(i,2)
         ia=i1
         do 550 k=1,2
         go to (505,515,525,535),ia
  505    sn(k)=1.-x-y-z
         go to 550
  515    sn(k)=z
         go to 550
  525    sn(k)=x
         go to 550
  535    sn(k)=y
  550    ia=i2
c-
c......test for corner nodes
c-
         if(i1 .ne. i2) then
c-
c......for midsides evaluate function etc
c-
            shap(i)=4.*sn(1)*sn(2)
         else
c-
c......corner node evaluation
c-
            shap(i)=sn(1)*(2.*sn(1)-1.)
         endif
      enddo
c-
c......final step
c-
      return
      end

c------------------------------------------------------------------xn3x
      subroutine xn3x(it,nen,x,y,z)
c---------------------------------------------------------------------c
c
c
c     subroutine to evaluate shape function and its derivatives
c                   for three dimensions
c
c      IT is element type 1 10 point tetrahedron
c                    type 2 15 point prism
c                    type 3 20 point parallelipiped
c                    type 4 13 point rectangular base pyramid
c      I  is shape function number
c      X, Y, Z  are cordinates of point to be evaluated in local coord
c
c---------------------------------------------------------------------c
ccc   implicit real*8 (a-h,o-z)
      save
      common /shape/shap(20),shpx(20),shpy(20),shpz(20)
      dimension irf(10,2)
      dimension a(5),b(5),c(5),xt(3),yt(3),jx(3),kx(3),xm(5)
     1         ,sn(2),snx(2),sny(2),snz(2),shpp(15,4),ilokup(13)
c
      data ilokup/5,6,1,2,3,4,9,7,14,15,10,11,13/
      data irf/1,1,2,2,3,3,1,2,3,4
     1        ,1,2,2,3,3,1,4,4,4,4/
      data xt/0.0,1.0,0.0/,yt/-1.0,0.0,1.0/,jx/2,3,1/,kx/3,1,2/
      data shpp/0.0,0.0,1.0,12*0.0,
     +      0.25,-1.0,3.0,-1.0,0.25,4*0.0,0.25,-1.0,0.0,-1.0,0.25,0.0,
     +      0.25,-1.0,0.0,1.0,-0.25,4*0.0,0.25,-1.0,0.0,1.0,-0.25,0.0,
     +      0.25,-1.0,0.0,-1.0,0.25,4*0.0,-0.25,1.0,0.0,1.0,-0.25,0.0/
      data ncall/0/
c
      if(it .eq. 4) go to 80
      if( it - 2 ) 500,80,300
c-
c-----shape functions for right prism and pyramid.....
c-
   80 do j0=1,nen
         if(it.eq.2) then
            i=j0
         else
            i=ilokup(j0)
         endif
         ncall = ncall + 1
         if( ncall .gt. 1 ) go to 125
c-
c-----calculate invarient triangular functions.....
c-
         n = 0
         do 100 j = 1,5,2
         n = n + 1
         jj = jx(n)
         kk = kx(n)
         a(j) = xt(jj)*yt(kk) - xt(kk)*yt(jj)
         b(j) = yt(jj) - yt(kk)
         c(j) = xt(kk) - xt(jj)
  100    continue
c-
c-.....shape function calculations.....
c-
  125    do 130 k = 1, 5, 2
            xm(k) = 0.5*(a(k)+b(k)*x+c(k)*y)
  130    continue
c-
c-.....set indexes and normalize coordinates.....
c-
         if( i .le. 6 ) then
            zooo = -z
         else
            zooo = z
         endif
         if(it .eq. 2) go to 132
         IF(X .EQ. 1.0) GO TO 250
c      if(x .eq. 1.0) x=0.999
         if(abs(1.0-x) .le. 1.0e-6) then
c           print *,'x=',x
            goto 250
         endif
         zooo=zooo/(1.-x)
c        if(x .eq. 1.) zooo=sign(x,zooo)
         if(i .eq. 3) go to 200
  132    if( i .gt. 6 .and. i .lt. 10 ) go to 170
         l = i
         if( i .gt. 9 ) l = i - 9
         if( mod(l,2) .eq. 0 ) go to 150
c-
c-.....corner nodes.....
c-
         shap(j0)=xm(l)*((2.0*xm(l)-1.0)*(1.0+zooo)-(1.0-zooo**2))/2.0
         if(it .eq. 4) shap(j0)=shap(j0)+xm(l)/2.0*x*(1.-zooo**2)
         shpx(j0)=0.5*b(l)*((1.+zooo)*(2.*xm(l)-0.5)-0.5*(1.0-zooo**2))
         if(it .eq. 4) then
           shpx(j0)=shpx(j0)+xm(l)/2.*((2.*xm(l)-1.)*zooo
     *            /(1.-x)+zooo**2+1.)+b(l)/4.*(1.-zooo**2)*x
         endif
         shpy(j0)=0.5*c(l)*((1.+zooo)*(2.*xm(l)-0.5)-0.5*(1.0-zooo**2))
         if(it .eq. 4) shpy(j0)=shpy(j0)+c(l)/4.*(1-zooo**2)*x
         shpz(j0) = xm(l)*(xm(l)+zooo-0.5)
         if(it .eq. 4) then
           shpz(j0)=(shpz(j0)-zooo*xm(l))/(1.-x)+xm(l)*zooo
         endif
         if(i .le. 6) shpz(j0) = -shpz(j0)
         goto 240
c-
c-.....mid side node.....
c-
  150    n1 = l - 1
         n2 = mod(l+1,6)
         shap(j0) = 2.0*xm(n1)*xm(n2)*(1.0+zooo)
         shpx(j0) = (1.0+zooo)*(xm(n1)*b(n2) + xm(n2)*b(n1))
         if(it .eq. 4) then
           shpx(j0)=shpx(j0)+2.*xm(n1)*xm(n2)*zooo/(1.-x)
         endif
         shpy(j0) = (1.0+zooo)*( xm(n1)*c(n2)+xm(n2)*c(n1))
         shpz(j0) = 2.0*xm(n1)*xm(n2)
         if(i .le. 6) shpz(j0) = -shpz(j0)
         if(it .eq. 2) goto 240
         shpz(j0)=shpz(j0)/(1.-x)
         goto 240
c-
c-.....mid side rectangle.....
c-
  170    n1 = 1
         if( i .eq. 8 ) n1 = 3
         if( i .eq. 9 ) n1 = 5
         shap(j0) = xm(n1)*(1.0-zooo**2)
         if(it .eq. 4) shap(j0)=shap(j0)*(1.-x)
         shpx(j0) = 0.5*(1.0-zooo**2)*b(n1)
         if(it .eq. 4) shpx(j0)=shpx(j0)*(1.-x)-xm(n1)*(zooo**2+1.)
         shpy(j0) = 0.5*(1.0-zooo**2)*c(n1)
         if(it .eq. 4) shpy(j0)=shpy(j0)*(1.-x)
         shpz(j0) = -2.0*xm(n1)*zooo
  180    goto 240
c-
c......special case for no 3 shape function of pyramid
c-
  200    l=i
         shap(j0)=xm(l)*(2.*xm(l)-1.0)
         shpx(j0)=b(l)/2.*(4.*xm(l)-1.0)
         shpy(j0)=c(l)/2.*(4.*xm(l)-1.0)
         shpz(j0)=0.0
         goto 240
c
c--------Special case when x exactly equals 1.0
c
  250    continue
         shap(j0)=shpp(i,1)
         shpx(j0)=shpp(i,2)
         shpy(j0)=shpp(i,3)
         shpy(j0)=shpp(i,4)
  240 enddo
      return
c-
c-.....shape functions for cube.....
c-
  300 zm1=0.125*(1.0-z)
      zp1=0.125*(1.0+z)
      ym1=0.125*(1.0-y)
      yp1=0.125*(1.0+y)
c-
c......corner nodes
c-
      shap(1)=(1.-x)*(1.-y)*zm1*(-x-y-z-2.)
      shpx(1)=(1.-y)*zm1*(2.*x+y+z+1.)
      shpy(1)=(1.-x)*zm1*(x+2.*y+z+1.)
      shpz(1)=(1.-x)*ym1*(x+y+2.*z+1.)

      shap(3)=(1.+x)*(1.-y)*zm1*(x-y-z-2.)
      shpx(3)=(1.-y)*zm1*(2.*x-y-z-1.)
      shpy(3)=(1.+x)*zm1*(-x+2.*y+z+1.)
      shpz(3)=(1.+x)*ym1*(-x+y+2.*z+1.)

      shap(5)=(1.+x)*(1.+y)*zm1*(x+y-z-2.)
      shpx(5)=(1.+y)*zm1*(2.*x+y-z-1.)
      shpy(5)=(1.+x)*zm1*(x+2.*y-z-1.)
      shpz(5)=(1.+x)*yp1*(-x-y+2.*z+1.)

      shap(7)=(1.-x)*(1.+y)*zm1*(-x+y-z-2.)
      shpx(7)=(1.+y)*zm1*(2.*x-y+z+1.)
      shpy(7)=(1.-x)*zm1*(-x+2.*y-z-1.)
      shpz(7)=(1.-x)*yp1*(x-y+2.*z+1.)

      shap(13)=(1.-x)*(1.-y)*zp1*(-x-y+z-2.)
      shpx(13)=(1.-y)*zp1*(2.*x+y-z+1.)
      shpy(13)=(1.-x)*zp1*(x+2.*y-z+1.)
      shpz(13)=(1.-x)*ym1*(-x-y+2.*z-1.)

      shap(15)=(1.+x)*(1.-y)*zp1*(x-y+z-2.)
      shpx(15)=(1.-y)*zp1*(2.*x-y+z-1.)
      shpy(15)=(1.+x)*zp1*(-x+2.*y-z+1.)
      shpz(15)=(1.+x)*ym1*(x-y+2.*z-1.)

      shap(17)=(1.+x)*(1.+y)*zp1*(x+y+z-2.)
      shpx(17)=(1.+y)*zp1*(2.*x+y+z-1.)
      shpy(17)=(1.+x)*zp1*(x+2.*y+z-1.)
      shpz(17)=(1.+x)*yp1*(x+y+2.*z-1.)

      shap(19)=(1.-x)*(1.+y)*zp1*(-x+y+z-2.)
      shpx(19)=(1.+y)*zp1*(2.*x-y-z+1.)
      shpy(19)=(1.-x)*zp1*(-x+2.*y+z-1.)
      shpz(19)=(1.-x)*yp1*(-x+y+2.*z-1.)
c-
c......mid-side nodes
c-
      x21=0.25*(1.0-x*x)
      y21=0.25*(1.0-y*y)
      z21=0.25*(1.0-z*z)
      x2=0.5*x
      y2=0.5*y
      z2=0.5*z
      shap(2)=x21*(1.-y)*(1.-z)
      shpx(2)=-x2*(1.-y)*(1.-z)
      shpy(2)=-x21*(1.-z)
      shpz(2)=-x21*(1.-y)

      shap(4)=(1.+x)*y21*(1.-z)
      shpx(4)=y21*(1.-z)
      shpy(4)=-y2*(1.+x)*(1.-z)
      shpz(4)=-y21*(1.+x)

      shap(6)=x21*(1.+y)*(1.-z)
      shpx(6)=-x2*(1.+y)*(1.-z)
      shpy(6)=x21*(1.-z)
      shpz(6)=-x21*(1.+y)

      shap(8)=(1.-x)*y21*(1.-z)
      shpx(8)=-y21*(1.-z)
      shpy(8)=-y2*(1.-x)*(1.-z)
      shpz(8)=-y21*(1.-x)

      shap(9)=(1.-x)*(1.-y)*z21
      shpx(9)=-(1.-y)*z21
      shpy(9)=-(1.-x)*z21
      shpz(9)=-z2*(1.-x)*(1.-y)

      shap(10)=(1.+x)*(1.-y)*z21
      shpx(10)=(1.-y)*z21
      shpy(10)=-(1.+x)*z21
      shpz(10)=-z2*(1.+x)*(1.-y)

      shap(11)=(1.+x)*(1.+y)*z21
      shpx(11)=(1.+y)*z21
      shpy(11)=(1.+x)*z21
      shpz(11)=-z2*(1.+x)*(1.+y)

      shap(12)=(1.-x)*(1.+y)*z21
      shpx(12)=-(1.+y)*z21
      shpy(12)=(1.-x)*z21
      shpz(12)=-z2*(1.-x)*(1.+y)

      shap(14)=x21*(1.-y)*(1.+z)
      shpx(14)=-x2*(1.-y)*(1.+z)
      shpy(14)=-x21*(1.+z)
      shpz(14)=x21*(1.-y)

      shap(16)=(1.+x)*y21*(1.+z)
      shpx(16)=y21*(1.+z)
      shpy(16)=-y2*(1.+x)*(1.+z)
      shpz(16)=y21*(1.+x)

      shap(18)=x21*(1.+y)*(1.+z)
      shpx(18)=-x2*(1.+y)*(1.+z)
      shpy(18)=x21*(1.+z)
      shpz(18)=x21*(1.+y)

      shap(20)=(1.-x)*y21*(1.+z)
      shpx(20)=-y21*(1.+z)
      shpy(20)=-y2*(1.-x)*(1.+z)
      shpz(20)=y21*(1.-x)

      return
  500 continue
c-
c......section for tetrahedron
c-
      do i=1,nen
         i1=irf(i,1)
         i2=irf(i,2)
         ia=i1
         do 550 k=1,2
         go to (505,515,525,535),ia
  505    sn(k)=1.-x-y-z
         snx(k)=-1.
         sny(k)=-1.
         snz(k)=-1.
         go to 550
  515    sn(k)=z
         snx(k)=0.
         sny(k)=0.
         snz(k)=1.
         go to 550
  525    sn(k)=x
         snx(k)=1.
         sny(k)=0.
         snz(k)=0.
         go to 550
  535    sn(k)=y
         snx(k)=0.
         sny(k)=1.
         snz(k)=0.
  550    ia=i2
c-
c......test for corner nodes
c-
         if(i1 .ne. i2) then
c-
c......for midsides evaluate function etc
c-
            shap(i)=4.*sn(1)*sn(2)
            shpx(i)=4.*(snx(1)*sn(2)+sn(1)*snx(2))
            shpy(i)=4.*(sny(1)*sn(2)+sn(1)*sny(2))
            shpz(i)=4.*(snz(1)*sn(2)+sn(1)*snz(2))
         else
c-
c......corner node evaluation
c-
            shap(i)=sn(1)*(2.*sn(1)-1.)
            shpx(i)=snx(1)*(4.*sn(1)-1.)
            shpy(i)=sny(1)*(4.*sn(1)-1.)
            shpz(i)=snz(1)*(4.*sn(1)-1.)
         endif
      enddo
c-
c......final step
c-
      return
      end