xref: /phasta/phSolver/incompressible/e3bvar.f (revision 595995161822a203c8467e0e4a253d7bd7d6df32)

      subroutine e3bvar (yl,      acl,     ul,
     &                   shpb,    shglb,
     &                   xlb,     lnode,
     &                   WdetJb,  bnorm,   pres,
     &                   u1,      u2,      u3,      rmu,
     &                   unm,     tau1n,   tau2n,   tau3n,
     &                   vdot,    rlKwall,
     &                   xKebe,   rKwall_glob)
c
c----------------------------------------------------------------------
c
c   This routine computes the variables at integration points for
c the boundary element routine.
c
c input:
c  yl     (npro,nshl,ndof)      : primitive variables (local)
c          ndof: 5[p,v1,v2,v3,T]+number of scalars solved
c  acl    (npro,nshl,ndof)      : acceleration (local)
c  ul     (npro,nshlb,nsd)       : displacement (local)
c  shpb   (nen)                 : boundary element shape-functions
c  shglb  (nsd,nen)             : boundary element grad-shape-functions
c  xlb    (npro,nenl,nsd)       : nodal coordinates at current step
c  lnode  (nenb)                : local nodes on the boundary
c
c output:
c  g1yi   (npro,ndof)           : grad-v in direction 1
c  g2yi   (npro,ndof)           : grad-v in direction 2
c  g3yi   (npro,ndof)           : grad-v in direction 3
c  WdetJb (npro)                : weighted Jacobian
c  bnorm  (npro,nsd)            : outward normal
c  pres   (npro)                : pressure
c  u1     (npro)                : x1-velocity component
c  u2     (npro)                : x2-velocity component
c  u3     (npro)                : x3-velocity component
c  unm    (npro)                : BC u dot n
c  p      (npro)                : BC pressure
c  tau1n  (npro)                : BC viscous flux 1
c  tau2n  (npro)                : BC viscous flux 2
c  tau3n  (npro)                : BC viscous flux 3
c  vdot   (npro,nsd)            : acceleration at quadrature points
c  rlKwall(npro,nshlb,nsd)      : wall stiffness contribution to the local residual
c
c Zdenek Johan, Summer 1990.  (Modified from e2bvar.f)
c Zdenek Johan, Winter 1991.  (Fortran 90)
c Alberto Figueroa, Winter 2004.  CMM-FSI
c----------------------------------------------------------------------
c
      use        turbsa
      include "common.h"
c
      dimension yl(npro,nshl,ndof),        rmu(npro),
     &            shpb(npro,nshl),           shglb(npro,nsd,nshl),
     &            xlb(npro,nenl,nsd),
     &            lnode(27),                 g1yi(npro,ndof),
     &            g2yi(npro,ndof),           g3yi(npro,ndof),
     &            WdetJb(npro),              bnorm(npro,nsd),
     &            pres(npro),
     &            u1(npro),                  u2(npro),
     &            u3(npro),
     &            unm(npro),
     &            tau1n(npro),               tau2n(npro),
     &            tau3n(npro),
     &            acl(npro,nshl,ndof),       ul(npro,nshl,nsd),
     &            vdot(npro,nsd),            rlKwall(npro,nshlb,nsd)
c
      dimension gl1yi(npro,ndof),          gl2yi(npro,ndof),
     &            gl3yi(npro,ndof),          dxdxib(npro,nsd,nsd),
     &            dxidxb(npro,nsd,nsd),      temp(npro),
     &            temp1(npro),               temp2(npro),
     &            temp3(npro),
     &            v1(npro,nsd),              v2(npro,nsd),
     &            v3(npro,nsd),
     &            rotnodallocal(npro,nsd,nsd),
     &            x1rot(npro,nsd),           x2rot(npro,nsd),
     &            x3rot(npro,nsd),           detJacrot(npro),
     &            B1(npro,5,3),              B2(npro,5,3),
     &            B3(npro,5,3),              Dmatrix(npro,5,5),
     &            DtimesB1(npro,5,3),        DtimesB2(npro,5,3),
     &            DtimesB3(npro,5,3),
     &            rKwall_local11(npro,nsd,nsd),
     &            rKwall_local12(npro,nsd,nsd),
     &            rKwall_local13(npro,nsd,nsd),
     &            rKwall_local21(npro,nsd,nsd),
     &            rKwall_local22(npro,nsd,nsd),
     &            rKwall_local23(npro,nsd,nsd),
     &            rKwall_local31(npro,nsd,nsd),
     &            rKwall_local32(npro,nsd,nsd),
     &            rKwall_local33(npro,nsd,nsd),
     &            rKwall_glob11(npro,nsd,nsd),
     &            rKwall_glob12(npro,nsd,nsd),
     &            rKwall_glob13(npro,nsd,nsd),
     &            rKwall_glob21(npro,nsd,nsd),
     &            rKwall_glob22(npro,nsd,nsd),
     &            rKwall_glob23(npro,nsd,nsd),
     &            rKwall_glob31(npro,nsd,nsd),
     &            rKwall_glob32(npro,nsd,nsd),
     &            rKwall_glob33(npro,nsd,nsd)
c
      dimension   rKwall_glob(npro,9,nshl,nshl),
     &            xKebe(npro,9,nshl,nshl)
c
      real*8      lhmFctvw, tsFctvw(npro)

      dimension   tmp1(npro)
c
      real*8    Turb(npro),                xki,
     &            xki3,                      fv1
c
      integer   e, i, j
c
      integer   aa, b


c
c.... ------------------->  integration variables  <--------------------
c
c.... compute the primitive variables at the integration point
c
      pres = zero
      u1   = zero
      u2   = zero
      u3   = zero
c

      do n = 1, nshlb
         nodlcl = lnode(n)
c
         pres = pres + shpb(:,nodlcl) * yl(:,nodlcl,1)
         u1   = u1   + shpb(:,nodlcl) * yl(:,nodlcl,2)
         u2   = u2   + shpb(:,nodlcl) * yl(:,nodlcl,3)
         u3   = u3   + shpb(:,nodlcl) * yl(:,nodlcl,4)

      enddo
c
c.... ---------------------->  Element Metrics  <-----------------------
c
c.... compute the deformation gradient
c
      dxdxib = zero
c
      do n = 1, nenl
         dxdxib(:,1,1) = dxdxib(:,1,1) + xlb(:,n,1) * shglb(:,1,n)
         dxdxib(:,1,2) = dxdxib(:,1,2) + xlb(:,n,1) * shglb(:,2,n)
         dxdxib(:,1,3) = dxdxib(:,1,3) + xlb(:,n,1) * shglb(:,3,n)
         dxdxib(:,2,1) = dxdxib(:,2,1) + xlb(:,n,2) * shglb(:,1,n)
         dxdxib(:,2,2) = dxdxib(:,2,2) + xlb(:,n,2) * shglb(:,2,n)
         dxdxib(:,2,3) = dxdxib(:,2,3) + xlb(:,n,2) * shglb(:,3,n)
         dxdxib(:,3,1) = dxdxib(:,3,1) + xlb(:,n,3) * shglb(:,1,n)
         dxdxib(:,3,2) = dxdxib(:,3,2) + xlb(:,n,3) * shglb(:,2,n)
         dxdxib(:,3,3) = dxdxib(:,3,3) + xlb(:,n,3) * shglb(:,3,n)
      enddo
c
c.... compute the normal to the boundary
c
c$$$      if (lcsyst .eq. 4) then   ! wedge-quad
c$$$         temp1 =  dxdxib(:,2,1) * dxdxib(:,3,3) -
c$$$     &            dxdxib(:,2,3) * dxdxib(:,3,1)
c$$$         temp2 =  dxdxib(:,3,1) * dxdxib(:,1,3) -
c$$$     &            dxdxib(:,3,3) * dxdxib(:,1,1)
c$$$         temp3 =  dxdxib(:,1,1) * dxdxib(:,2,3) -
c$$$     &            dxdxib(:,1,3) * dxdxib(:,2,1)
c$$$      elseif( lcyst .eq. 6) then  ! pyr-tri face
c$$$         temp1 =  dxdxib(:,2,1) * dxdxib(:,3,3) -
c$$$     &            dxdxib(:,2,3) * dxdxib(:,3,1)
c$$$         temp2 =  dxdxib(:,3,1) * dxdxib(:,1,3) -
c$$$     &            dxdxib(:,3,3) * dxdxib(:,1,1)
c$$$         temp3 =  dxdxib(:,1,1) * dxdxib(:,2,3) -
c$$$     &            dxdxib(:,1,3) * dxdxib(:,2,1)
c$$$      elseif( lcyst .eq. 1) then !usual wrong way tets
c$$$         temp1 = -dxdxib(:,2,2) * dxdxib(:,3,1) +
c$$$     &            dxdxib(:,2,1) * dxdxib(:,3,2)
c$$$         temp2 = -dxdxib(:,3,2) * dxdxib(:,1,1) +
c$$$     &            dxdxib(:,3,1) * dxdxib(:,1,2)
c$$$         temp3 = -dxdxib(:,1,2) * dxdxib(:,2,1) +
c$$$     &            dxdxib(:,1,1) * dxdxib(:,2,2)
c$$$      else
c$$$         temp1 =  dxdxib(:,2,2) * dxdxib(:,3,1) -
c$$$     &            dxdxib(:,2,1) * dxdxib(:,3,2)
c$$$         temp2 =  dxdxib(:,3,2) * dxdxib(:,1,1) -
c$$$     &            dxdxib(:,3,1) * dxdxib(:,1,2)
c$$$         temp3 =  dxdxib(:,1,2) * dxdxib(:,2,1) -
c$$$     &            dxdxib(:,1,1) * dxdxib(:,2,2)
c$$$      endif
c$$$c
c$$$      temp       = one / sqrt ( temp1**2 + temp2**2 + temp3**2 )
c$$$      bnorm(:,1) = temp1 * temp
c$$$      bnorm(:,2) = temp2 * temp
c$$$      bnorm(:,3) = temp3 * temp
c$$$c
c$$$      WdetJb     = Qwtb(lcsyst,intp) / temp
c$$$      if(lcsyst .eq. 3) WdetJb=WdetJb*two
c
c.... compute the normal to the boundary. This is achieved by taking
c     the cross product of two vectors in the plane of the 2-d
c     boundary face.
c
      if(lcsyst.eq.1) then      ! set to curl into element all others out
         ipt2=2
         ipt3=3
      elseif(lcsyst.eq.2) then
         ipt2=4
         ipt3=2
      elseif(lcsyst.eq.3) then
         ipt2=3
         ipt3=2
      elseif(lcsyst.eq.4) then
         ipt2=2
         ipt3=4
      elseif(lcsyst.eq.5) then
         ipt2=4
         ipt3=2
      elseif(lcsyst.eq.6) then
         ipt2=2
         ipt3=5
      endif
      v1 = xlb(:,ipt2,:) - xlb(:,1,:)
      v2 = xlb(:,ipt3,:) - xlb(:,1,:)
c
c compute cross product
c
      temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3)
      temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3)
      temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2)
c
c mag is area for quads, twice area for tris
c
      temp       = one / sqrt ( temp1**2 + temp2**2 + temp3**2 )
      bnorm(:,1) = temp1 * temp
      bnorm(:,2) = temp2 * temp
      bnorm(:,3) = temp3 * temp
c

      if (lcsyst .eq. 1) then
         WdetJb     = Qwtb(lcsyst,intp) / (four*temp)
      elseif (lcsyst .eq. 2) then
         WdetJb     = Qwtb(lcsyst,intp) / (four*temp)
      elseif (lcsyst .eq. 3) then
         WdetJb     = Qwtb(lcsyst,intp) / (two*temp)
      elseif (lcsyst .eq. 4) then
         WdetJb     = Qwtb(lcsyst,intp) / (four*temp)
      elseif (lcsyst .eq. 5) then
         WdetJb     = Qwtb(lcsyst,intp) / (four*temp)
      elseif (lcsyst .eq. 6) then
         WdetJb     = Qwtb(lcsyst,intp) / (two*temp)
      endif
c
c.... -------------------------->  Grad-V  <----------------------------
c
c.... compute grad-v for Navier-Stokes terms
c
      if (Navier .eq. 1) then
c
c.... compute the inverse of deformation gradient
c
         dxidxb(:,1,1) =   dxdxib(:,2,2) * dxdxib(:,3,3)
     &        - dxdxib(:,3,2) * dxdxib(:,2,3)
         dxidxb(:,1,2) =   dxdxib(:,3,2) * dxdxib(:,1,3)
     &        - dxdxib(:,1,2) * dxdxib(:,3,3)
         dxidxb(:,1,3) =   dxdxib(:,1,2) * dxdxib(:,2,3)
     &        - dxdxib(:,1,3) * dxdxib(:,2,2)
         temp          = one / ( dxidxb(:,1,1) * dxdxib(:,1,1)
     &        + dxidxb(:,1,2) * dxdxib(:,2,1)
     &        + dxidxb(:,1,3) * dxdxib(:,3,1) )
         dxidxb(:,1,1) =  dxidxb(:,1,1) * temp
         dxidxb(:,1,2) =  dxidxb(:,1,2) * temp
         dxidxb(:,1,3) =  dxidxb(:,1,3) * temp
         dxidxb(:,2,1) = (dxdxib(:,2,3) * dxdxib(:,3,1)
     &        - dxdxib(:,2,1) * dxdxib(:,3,3)) * temp
         dxidxb(:,2,2) = (dxdxib(:,1,1) * dxdxib(:,3,3)
     &        - dxdxib(:,3,1) * dxdxib(:,1,3)) * temp
         dxidxb(:,2,3) = (dxdxib(:,2,1) * dxdxib(:,1,3)
     &        - dxdxib(:,1,1) * dxdxib(:,2,3)) * temp
         dxidxb(:,3,1) = (dxdxib(:,2,1) * dxdxib(:,3,2)
     &        - dxdxib(:,2,2) * dxdxib(:,3,1)) * temp
         dxidxb(:,3,2) = (dxdxib(:,3,1) * dxdxib(:,1,2)
     &        - dxdxib(:,1,1) * dxdxib(:,3,2)) * temp
         dxidxb(:,3,3) = (dxdxib(:,1,1) * dxdxib(:,2,2)
     &        - dxdxib(:,1,2) * dxdxib(:,2,1)) * temp
c
c.... compute local-grad-Y
c
         gl1yi = zero
         gl2yi = zero
         gl3yi = zero
c
         do n = 1, nshl
            gl1yi(:,1) = gl1yi(:,1) + shglb(:,1,n) * yl(:,n,1)
            gl1yi(:,2) = gl1yi(:,2) + shglb(:,1,n) * yl(:,n,2)
            gl1yi(:,3) = gl1yi(:,3) + shglb(:,1,n) * yl(:,n,3)
            gl1yi(:,4) = gl1yi(:,4) + shglb(:,1,n) * yl(:,n,4)
c
            gl2yi(:,1) = gl2yi(:,1) + shglb(:,2,n) * yl(:,n,1)
            gl2yi(:,2) = gl2yi(:,2) + shglb(:,2,n) * yl(:,n,2)
            gl2yi(:,3) = gl2yi(:,3) + shglb(:,2,n) * yl(:,n,3)
            gl2yi(:,4) = gl2yi(:,4) + shglb(:,2,n) * yl(:,n,4)
c
            gl3yi(:,1) = gl3yi(:,1) + shglb(:,3,n) * yl(:,n,1)
            gl3yi(:,2) = gl3yi(:,2) + shglb(:,3,n) * yl(:,n,2)
            gl3yi(:,3) = gl3yi(:,3) + shglb(:,3,n) * yl(:,n,3)
            gl3yi(:,4) = gl3yi(:,4) + shglb(:,3,n) * yl(:,n,4)
         enddo
c
c.... convert local-grads to global-grads
c
         g1yi(:,2) = dxidxb(:,1,1) * gl1yi(:,2) +
     &        dxidxb(:,2,1) * gl2yi(:,2) +
     &        dxidxb(:,3,1) * gl3yi(:,2)
         g2yi(:,2) = dxidxb(:,1,2) * gl1yi(:,2) +
     &        dxidxb(:,2,2) * gl2yi(:,2) +
     &        dxidxb(:,3,2) * gl3yi(:,2)
         g3yi(:,2) = dxidxb(:,1,3) * gl1yi(:,2) +
     &        dxidxb(:,2,3) * gl2yi(:,2) +
     &        dxidxb(:,3,3) * gl3yi(:,2)
c
         g1yi(:,3) = dxidxb(:,1,1) * gl1yi(:,3) +
     &        dxidxb(:,2,1) * gl2yi(:,3) +
     &        dxidxb(:,3,1) * gl3yi(:,3)
         g2yi(:,3) = dxidxb(:,1,2) * gl1yi(:,3) +
     &        dxidxb(:,2,2) * gl2yi(:,3) +
     &        dxidxb(:,3,2) * gl3yi(:,3)
         g3yi(:,3) = dxidxb(:,1,3) * gl1yi(:,3) +
     &        dxidxb(:,2,3) * gl2yi(:,3) +
     &        dxidxb(:,3,3) * gl3yi(:,3)
c
         g1yi(:,4) = dxidxb(:,1,1) * gl1yi(:,4) +
     &        dxidxb(:,2,1) * gl2yi(:,4) +
     &        dxidxb(:,3,1) * gl3yi(:,4)
         g2yi(:,4) = dxidxb(:,1,2) * gl1yi(:,4) +
     &        dxidxb(:,2,2) * gl2yi(:,4) +
     &        dxidxb(:,3,2) * gl3yi(:,4)
         g3yi(:,4) = dxidxb(:,1,3) * gl1yi(:,4) +
     &        dxidxb(:,2,3) * gl2yi(:,4) +
     &        dxidxb(:,3,3) * gl3yi(:,4)
c
c.... end grad-v
c
      endif

c
c.... mass flux
c
      unm = bnorm(:,1) * u1 +bnorm(:,2) * u2  +bnorm(:,3) * u3
! no rho in continuity eq.


c
c.... viscous flux
c
      tau1n = bnorm(:,1) * two * rmu *  g1yi(:,2)
     &     + bnorm(:,2) *      (rmu * (g2yi(:,2) + g1yi(:,3)))
     &     + bnorm(:,3) *      (rmu * (g3yi(:,2) + g1yi(:,4)))
      tau2n = bnorm(:,1) *      (rmu * (g2yi(:,2) + g1yi(:,3)))
     &     + bnorm(:,2) * two * rmu *  g2yi(:,3)
     &     + bnorm(:,3) *      (rmu * (g3yi(:,3) + g2yi(:,4)))
      tau3n = bnorm(:,1) *      (rmu * (g3yi(:,2) + g1yi(:,4)))
     &     + bnorm(:,2) *      (rmu * (g3yi(:,3) + g2yi(:,4)))
     &     + bnorm(:,3) * two * rmu *  g3yi(:,4)
c
      temp1 = bnorm(:,1) * tau1n
     &     + bnorm(:,2) * tau2n
     &     + bnorm(:,3) * tau3n

      pres  = pres - temp1

      tau1n = tau1n - bnorm(:,1) * temp1
      tau2n = tau2n - bnorm(:,2) * temp1
      tau3n = tau3n - bnorm(:,3) * temp1

c
c.... viscous flux control
c
c     if iviscflux = 1, we consider the viscous flux on the RHS
c     otherwise, if iviscflux = 0, we eliminate this term: stability in
c     pressure-flow coupled boundaries for cardiovascular applications in
c     situations of flow reversal
      tau1n = tau1n * iviscflux
      tau2n = tau2n * iviscflux
      tau3n = tau3n * iviscflux

      vdot = zero
      rlKwall = zero
      if (intp.eq.ngaussb)   then    ! do this only for the last gauss point
        rKwall_glob = zero
      endif

      if(ideformwall.eq.1) then
      do n = 1, nshlb
         nodlcl = lnode(n)
c
         vdot(:,1) = vdot(:,1) + shpb(:,nodlcl) * acl(:,nodlcl,2)
         vdot(:,2) = vdot(:,2) + shpb(:,nodlcl) * acl(:,nodlcl,3)
         vdot(:,3) = vdot(:,3) + shpb(:,nodlcl) * acl(:,nodlcl,4)

      enddo
      vdot = vdot * thicknessvw * rhovw
c
c.... --------------------->  Stiffness matrix & residual  <-----------------
c
c.... B^t * D * B formulation for plane stress enhanced membrane
c
c
c.... rotation matrix
c
      v1 = xlb(:,ipt2,:) - xlb(:,1,:)
      temp       = one / sqrt ( v1(:,1)**2 + v1(:,2)**2 + v1(:,3)**2 )
      v1(:,1) = v1(:,1) * temp
      v1(:,2) = v1(:,2) * temp
      v1(:,3) = v1(:,3) * temp

      v2 = xlb(:,ipt3,:) - xlb(:,1,:)

c     compute cross product
      temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3)
      temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3)
      temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2)

      temp       = one / sqrt ( temp1**2 + temp2**2 + temp3**2 )
      v3(:,1) = temp1 * temp
      v3(:,2) = temp2 * temp
      v3(:,3) = temp3 * temp

c     cross product again for v2
      temp1 = v3(:,2) * v1(:,3) - v1(:,2) * v3(:,3)
      temp2 = v1(:,1) * v3(:,3) - v3(:,1) * v1(:,3)
      temp3 = v3(:,1) * v1(:,2) - v1(:,1) * v3(:,2)

      temp       = one / sqrt ( temp1**2 + temp2**2 + temp3**2 )
      v2(:,1) = temp1 * temp
      v2(:,2) = temp2 * temp
      v2(:,3) = temp3 * temp

      do j = 1, nsd
         rotnodallocal(:,1,j) = v1(:,j)
         rotnodallocal(:,2,j) = v2(:,j)
         rotnodallocal(:,3,j) = v3(:,j)
      enddo

c
c.... rotated coordinates
c
      x1rot = zero
      x2rot = zero
      x3rot = zero

      do i = 1, nsd
         do j = 1, nsd
            x1rot(:,i) = x1rot(:,i)+rotnodallocal(:,i,j)*xlb(:,1,j)
            x2rot(:,i) = x2rot(:,i)+rotnodallocal(:,i,j)*xlb(:,ipt2,j)
            x3rot(:,i) = x3rot(:,i)+rotnodallocal(:,i,j)*xlb(:,ipt3,j)
         enddo
      enddo

c
c.... B matrices
c
      B1 = zero
      B2 = zero
      B3 = zero
      detJacrot = (x2rot(:,1)-x1rot(:,1)) * (x3rot(:,2)-x1rot(:,2)) -
     &     (x3rot(:,1)-x1rot(:,1)) * (x2rot(:,2)-x1rot(:,2))

      B1(:,1,1) = (x2rot(:,2)-x3rot(:,2))/detJacrot(:)
      B1(:,2,2) = (x3rot(:,1)-x2rot(:,1))/detJacrot(:)
      B1(:,3,1) = (x3rot(:,1)-x2rot(:,1))/detJacrot(:)
      B1(:,3,2) = (x2rot(:,2)-x3rot(:,2))/detJacrot(:)
      B1(:,4,3) = (x2rot(:,2)-x3rot(:,2))/detJacrot(:)
      B1(:,5,3) = (x3rot(:,1)-x2rot(:,1))/detJacrot(:)

      B2(:,1,1) = (x3rot(:,2)-x1rot(:,2))/detJacrot(:)
      B2(:,2,2) = (x1rot(:,1)-x3rot(:,1))/detJacrot(:)
      B2(:,3,1) = (x1rot(:,1)-x3rot(:,1))/detJacrot(:)
      B2(:,3,2) = (x3rot(:,2)-x1rot(:,2))/detJacrot(:)
      B2(:,4,3) = (x3rot(:,2)-x1rot(:,2))/detJacrot(:)
      B2(:,5,3) = (x1rot(:,1)-x3rot(:,1))/detJacrot(:)

      B3(:,1,1) = (x1rot(:,2)-x2rot(:,2))/detJacrot(:)
      B3(:,2,2) = (x2rot(:,1)-x1rot(:,1))/detJacrot(:)
      B3(:,3,1) = (x2rot(:,1)-x1rot(:,1))/detJacrot(:)
      B3(:,3,2) = (x1rot(:,2)-x2rot(:,2))/detJacrot(:)
      B3(:,4,3) = (x1rot(:,2)-x2rot(:,2))/detJacrot(:)
      B3(:,5,3) = (x2rot(:,1)-x1rot(:,1))/detJacrot(:)

C      B1 = B1 / detJacrot
C      B2 = B2 / detJacrot
C      B3 = B3 / detJacrot

c
c.... D matrix
c
      Dmatrix = zero
      temp1 = evw / (1.0d0 - rnuvw*rnuvw)
      temp2 = rnuvw * temp1
      temp3 = pt5 * (1.0d0 - rnuvw) * temp1
      Dmatrix(:,1,1) = temp1
      Dmatrix(:,1,2) = temp2
      Dmatrix(:,2,1) = temp2
      Dmatrix(:,2,2) = temp1
      Dmatrix(:,3,3) = temp3
      Dmatrix(:,4,4) = temp3*rshearconstantvw
      Dmatrix(:,5,5) = temp3*rshearconstantvw
c
c.... D * [B1|B2|B3]
c
      DtimesB1 = zero
      DtimesB2 = zero
      DtimesB3 = zero
      do i = 1, 5
         do j = 1, 3
            do k = 1, 5
               DtimesB1(:,i,j) = DtimesB1(:,i,j)
     &              + Dmatrix(:,i,k) * B1(:,k,j)
               DtimesB2(:,i,j) = DtimesB2(:,i,j)
     &              + Dmatrix(:,i,k) * B2(:,k,j)
               DtimesB3(:,i,j) = DtimesB3(:,i,j)
     &              + Dmatrix(:,i,k) * B3(:,k,j)
            enddo
         enddo
      enddo
c
c.... [B1|B2|B3]^T * D * [B1|B2|B3]
c
      rKwall_local11 = zero
      rKwall_local12 = zero
      rKwall_local13 = zero
      rKwall_local21 = zero
      rKwall_local22 = zero
      rKwall_local23 = zero
      rKwall_local31 = zero
      rKwall_local32 = zero
      rKwall_local33 = zero

      do i = 1, 3               ! i is a node index: i=1, nenbl=3
         do j = 1, 3            ! same is true for j
            do k = 1, 5
               rKwall_local11(:,i,j)  = rKwall_local11(:,i,j)
     &              + B1(:,k,i) * DtimesB1(:,k,j)
               rKwall_local12(:,i,j)  = rKwall_local12(:,i,j)
     &              + B1(:,k,i) * DtimesB2(:,k,j)
               rKwall_local13(:,i,j)  = rKwall_local13(:,i,j)
     &              + B1(:,k,i) * DtimesB3(:,k,j)
               rKwall_local21(:,i,j)  = rKwall_local21(:,i,j)
     &              + B2(:,k,i) * DtimesB1(:,k,j)
               rKwall_local22(:,i,j)  = rKwall_local22(:,i,j)
     &              + B2(:,k,i) * DtimesB2(:,k,j)
               rKwall_local23(:,i,j)  = rKwall_local23(:,i,j)
     &              + B2(:,k,i) * DtimesB3(:,k,j)
               rKwall_local31(:,i,j)  = rKwall_local31(:,i,j)
     &              + B3(:,k,i) * DtimesB1(:,k,j)
               rKwall_local32(:,i,j)  = rKwall_local32(:,i,j)
     &              + B3(:,k,i) * DtimesB2(:,k,j)
               rKwall_local33(:,i,j)  = rKwall_local33(:,i,j)
     &              + B3(:,k,i) * DtimesB3(:,k,j)
            enddo
         enddo
      enddo

c
c.... Now we need to rotate each of these submatrices to the global frame
c
      call rotatestiff(rKwall_local11, rotnodallocal, rKwall_glob11)
      call rotatestiff(rKwall_local12, rotnodallocal, rKwall_glob12)
      call rotatestiff(rKwall_local13, rotnodallocal, rKwall_glob13)
      call rotatestiff(rKwall_local21, rotnodallocal, rKwall_glob21)
      call rotatestiff(rKwall_local22, rotnodallocal, rKwall_glob22)
      call rotatestiff(rKwall_local23, rotnodallocal, rKwall_glob23)
      call rotatestiff(rKwall_local31, rotnodallocal, rKwall_glob31)
      call rotatestiff(rKwall_local32, rotnodallocal, rKwall_glob32)
      call rotatestiff(rKwall_local33, rotnodallocal, rKwall_glob33)

c     multiply the nodal matrices by the area and the thickness
      do i =1, nsd
         do j = 1, nsd
            rKwall_glob11(:,i,j) = rKwall_glob11(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
            rKwall_glob12(:,i,j) = rKwall_glob12(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
            rKwall_glob13(:,i,j) = rKwall_glob13(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
            rKwall_glob21(:,i,j) = rKwall_glob21(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
            rKwall_glob22(:,i,j) = rKwall_glob22(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
            rKwall_glob23(:,i,j) = rKwall_glob23(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
            rKwall_glob31(:,i,j) = rKwall_glob31(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
            rKwall_glob32(:,i,j) = rKwall_glob32(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
            rKwall_glob33(:,i,j) = rKwall_glob33(:,i,j) * detJacrot(:)
     &                           * pt5 * thicknessvw
         enddo
      enddo

c
c.... Final K * u product (in global coordinates) to get the residual
c
      do i = 1, 3               ! now i is a spatial index: i=1, nsd=3
         rlKwall(:,1,1) = rlKwall(:,1,1)
     &                  + rKwall_glob11(:,1,i) * ul(:,1,i)
     &                  + rKwall_glob12(:,1,i) * ul(:,2,i)
     &                  + rKwall_glob13(:,1,i) * ul(:,3,i)
         rlKwall(:,1,2) = rlKwall(:,1,2)
     &                  + rKwall_glob11(:,2,i) * ul(:,1,i)
     &                  + rKwall_glob12(:,2,i) * ul(:,2,i)
     &                  + rKwall_glob13(:,2,i) * ul(:,3,i)
         rlKwall(:,1,3) = rlKwall(:,1,3)
     &                  + rKwall_glob11(:,3,i) * ul(:,1,i)
     &                  + rKwall_glob12(:,3,i) * ul(:,2,i)
     &                  + rKwall_glob13(:,3,i) * ul(:,3,i)
         rlKwall(:,2,1) = rlKwall(:,2,1)
     &                  + rKwall_glob21(:,1,i) * ul(:,1,i)
     &                  + rKwall_glob22(:,1,i) * ul(:,2,i)
     &                  + rKwall_glob23(:,1,i) * ul(:,3,i)
         rlKwall(:,2,2) = rlKwall(:,2,2)
     &                  + rKwall_glob21(:,2,i) * ul(:,1,i)
     &                  + rKwall_glob22(:,2,i) * ul(:,2,i)
     &                  + rKwall_glob23(:,2,i) * ul(:,3,i)
         rlKwall(:,2,3) = rlKwall(:,2,3)
     &                  + rKwall_glob21(:,3,i) * ul(:,1,i)
     &                  + rKwall_glob22(:,3,i) * ul(:,2,i)
     &                  + rKwall_glob23(:,3,i) * ul(:,3,i)
         rlKwall(:,3,1) = rlKwall(:,3,1)
     &                  + rKwall_glob31(:,1,i) * ul(:,1,i)
     &                  + rKwall_glob32(:,1,i) * ul(:,2,i)
     &                  + rKwall_glob33(:,1,i) * ul(:,3,i)
         rlKwall(:,3,2) = rlKwall(:,3,2)
     &                  + rKwall_glob31(:,2,i) * ul(:,1,i)
     &                  + rKwall_glob32(:,2,i) * ul(:,2,i)
     &                  + rKwall_glob33(:,2,i) * ul(:,3,i)
         rlKwall(:,3,3) = rlKwall(:,3,3)
     &                  + rKwall_glob31(:,3,i) * ul(:,1,i)
     &                  + rKwall_glob32(:,3,i) * ul(:,2,i)
     &                  + rKwall_glob33(:,3,i) * ul(:,3,i)
      enddo
c
c.... --------------> End of Stiffness matrix & residual  <-----------------
c

c
c.... -----> Wall Stiffness and Mass matrices for implicit LHS  <-----------
c

c....  Here we just add the mass matrix contribution.  The stiffness contribution
c....  is added in e3b

c      lhmFct = almi * (one - flmpl)      Maybe we have to define flmplW: lumped
                                        ! mass parameter for the wall
      lhmFctvw = almi * (one - flmpl)
c
c.... scale variables for efficiency
c
      tsFctvw     = lhmFctvw * WdetJb * rhovw * thicknessvw
c
c.... compute mass and convection terms
c
c.... NOTE:  the wall mass contributions should only have 3 nodal components
c.... since the fourth node is an interior node... therefore, the loops should
c.... be done from 1 to nshlb=3...

      do b = 1, nshlb
         do aa = 1, nshlb
            tmp1 = tsFctvw * shpb(:,aa) * shpb(:,b)
c
c           tmp1=alpha_m*(1-lmp)*WdetJ*N^aN^b*rho*thickness   the time term
c
            xKebe(:,1,aa,b) = xKebe(:,1,aa,b) + tmp1
            xKebe(:,5,aa,b) = xKebe(:,5,aa,b) + tmp1
            xKebe(:,9,aa,b) = xKebe(:,9,aa,b) + tmp1
         enddo
      enddo

c
c.... assemble the nodal stiffness into the element stiffness matrix rKwall_glob
c
c.... We have passed the integer intp to make this operation only once: we are
c.... not using the gauss points structure to compute the stiffness of the wall
c.... elements, so we don't want to be redundant and calculate ngaussb times the
c.... stiffness matrix which is constant for linear triangles...

c.... This is ugly, but I will fix it later...

      if (intp.eq.ngaussb)   then    ! do this only for the last gauss point
        rKwall_glob(:,1,1,1) = rKwall_glob11(:,1,1)
        rKwall_glob(:,2,1,1) = rKwall_glob11(:,1,2)
        rKwall_glob(:,3,1,1) = rKwall_glob11(:,1,3)
        rKwall_glob(:,4,1,1) = rKwall_glob11(:,2,1)
        rKwall_glob(:,5,1,1) = rKwall_glob11(:,2,2)
        rKwall_glob(:,6,1,1) = rKwall_glob11(:,2,3)
        rKwall_glob(:,7,1,1) = rKwall_glob11(:,3,1)
        rKwall_glob(:,8,1,1) = rKwall_glob11(:,3,2)
        rKwall_glob(:,9,1,1) = rKwall_glob11(:,3,3)

        rKwall_glob(:,1,1,2) = rKwall_glob12(:,1,1)
        rKwall_glob(:,2,1,2) = rKwall_glob12(:,1,2)
        rKwall_glob(:,3,1,2) = rKwall_glob12(:,1,3)
        rKwall_glob(:,4,1,2) = rKwall_glob12(:,2,1)
        rKwall_glob(:,5,1,2) = rKwall_glob12(:,2,2)
        rKwall_glob(:,6,1,2) = rKwall_glob12(:,2,3)
        rKwall_glob(:,7,1,2) = rKwall_glob12(:,3,1)
        rKwall_glob(:,8,1,2) = rKwall_glob12(:,3,2)
        rKwall_glob(:,9,1,2) = rKwall_glob12(:,3,3)

        rKwall_glob(:,1,1,3) = rKwall_glob13(:,1,1)
        rKwall_glob(:,2,1,3) = rKwall_glob13(:,1,2)
        rKwall_glob(:,3,1,3) = rKwall_glob13(:,1,3)
        rKwall_glob(:,4,1,3) = rKwall_glob13(:,2,1)
        rKwall_glob(:,5,1,3) = rKwall_glob13(:,2,2)
        rKwall_glob(:,6,1,3) = rKwall_glob13(:,2,3)
        rKwall_glob(:,7,1,3) = rKwall_glob13(:,3,1)
        rKwall_glob(:,8,1,3) = rKwall_glob13(:,3,2)
        rKwall_glob(:,9,1,3) = rKwall_glob13(:,3,3)

        rKwall_glob(:,1,2,1) = rKwall_glob21(:,1,1)
        rKwall_glob(:,2,2,1) = rKwall_glob21(:,1,2)
        rKwall_glob(:,3,2,1) = rKwall_glob21(:,1,3)
        rKwall_glob(:,4,2,1) = rKwall_glob21(:,2,1)
        rKwall_glob(:,5,2,1) = rKwall_glob21(:,2,2)
        rKwall_glob(:,6,2,1) = rKwall_glob21(:,2,3)
        rKwall_glob(:,7,2,1) = rKwall_glob21(:,3,1)
        rKwall_glob(:,8,2,1) = rKwall_glob21(:,3,2)
        rKwall_glob(:,9,2,1) = rKwall_glob21(:,3,3)

        rKwall_glob(:,1,2,2) = rKwall_glob22(:,1,1)
        rKwall_glob(:,2,2,2) = rKwall_glob22(:,1,2)
        rKwall_glob(:,3,2,2) = rKwall_glob22(:,1,3)
        rKwall_glob(:,4,2,2) = rKwall_glob22(:,2,1)
        rKwall_glob(:,5,2,2) = rKwall_glob22(:,2,2)
        rKwall_glob(:,6,2,2) = rKwall_glob22(:,2,3)
        rKwall_glob(:,7,2,2) = rKwall_glob22(:,3,1)
        rKwall_glob(:,8,2,2) = rKwall_glob22(:,3,2)
        rKwall_glob(:,9,2,2) = rKwall_glob22(:,3,3)

        rKwall_glob(:,1,2,3) = rKwall_glob23(:,1,1)
        rKwall_glob(:,2,2,3) = rKwall_glob23(:,1,2)
        rKwall_glob(:,3,2,3) = rKwall_glob23(:,1,3)
        rKwall_glob(:,4,2,3) = rKwall_glob23(:,2,1)
        rKwall_glob(:,5,2,3) = rKwall_glob23(:,2,2)
        rKwall_glob(:,6,2,3) = rKwall_glob23(:,2,3)
        rKwall_glob(:,7,2,3) = rKwall_glob23(:,3,1)
        rKwall_glob(:,8,2,3) = rKwall_glob23(:,3,2)
        rKwall_glob(:,9,2,3) = rKwall_glob23(:,3,3)

        rKwall_glob(:,1,3,1) = rKwall_glob31(:,1,1)
        rKwall_glob(:,2,3,1) = rKwall_glob31(:,1,2)
        rKwall_glob(:,3,3,1) = rKwall_glob31(:,1,3)
        rKwall_glob(:,4,3,1) = rKwall_glob31(:,2,1)
        rKwall_glob(:,5,3,1) = rKwall_glob31(:,2,2)
        rKwall_glob(:,6,3,1) = rKwall_glob31(:,2,3)
        rKwall_glob(:,7,3,1) = rKwall_glob31(:,3,1)
        rKwall_glob(:,8,3,1) = rKwall_glob31(:,3,2)
        rKwall_glob(:,9,3,1) = rKwall_glob31(:,3,3)

        rKwall_glob(:,1,3,2) = rKwall_glob32(:,1,1)
        rKwall_glob(:,2,3,2) = rKwall_glob32(:,1,2)
        rKwall_glob(:,3,3,2) = rKwall_glob32(:,1,3)
        rKwall_glob(:,4,3,2) = rKwall_glob32(:,2,1)
        rKwall_glob(:,5,3,2) = rKwall_glob32(:,2,2)
        rKwall_glob(:,6,3,2) = rKwall_glob32(:,2,3)
        rKwall_glob(:,7,3,2) = rKwall_glob32(:,3,1)
        rKwall_glob(:,8,3,2) = rKwall_glob32(:,3,2)
        rKwall_glob(:,9,3,2) = rKwall_glob32(:,3,3)

        rKwall_glob(:,1,3,3) = rKwall_glob33(:,1,1)
        rKwall_glob(:,2,3,3) = rKwall_glob33(:,1,2)
        rKwall_glob(:,3,3,3) = rKwall_glob33(:,1,3)
        rKwall_glob(:,4,3,3) = rKwall_glob33(:,2,1)
        rKwall_glob(:,5,3,3) = rKwall_glob33(:,2,2)
        rKwall_glob(:,6,3,3) = rKwall_glob33(:,2,3)
        rKwall_glob(:,7,3,3) = rKwall_glob33(:,3,1)
        rKwall_glob(:,8,3,3) = rKwall_glob33(:,3,2)
        rKwall_glob(:,9,3,3) = rKwall_glob33(:,3,3)

        rKwall_glob = rKwall_glob*betai*Delt(itseq)*Delt(itseq)*alfi

      else
c....   nothing happens
        goto 123
      endif

123   continue

      endif
c
c.... return
c
      return
      end

c---------------------------------------------------------------------
c
c     variables for boundary elements
c
c---------------------------------------------------------------------
        subroutine e3bvarSclr (yl,        shdrv,    xlb,
     &                         shape,     WdetJb,   bnorm,
     &                         flux,      dwl )

        include "common.h"
c
        dimension yl(npro,nshl,ndof),        shdrv(npro,nsd,nshl),
     &            xlb(npro,nenl,nsd),        shape(npro,nshl),
     &            WdetJb(npro),              bnorm(npro,nsd),
     &            flux(npro)
c
        dimension dxdxib(npro,nsd,nsd),
     &            dxidxb(npro,nsd,nsd),      temp(npro),
     &            temp1(npro),               temp2(npro),
     &            temp3(npro),
     &            v1(npro,nsd),              v2(npro,nsd),
     &            gradSl(npro,nsd),          gradS(npro,nsd)

        real*8    diffus(npro),              dwl(npro,nshl)

        call getdiffsclr(shape,dwl,yl,diffus)
c
c.... ---------------------->  Element Metrics  <-----------------------
c
c.... compute the deformation gradient
c
        dxdxib = zero
c
        do n = 1, nenl
           dxdxib(:,1,1) = dxdxib(:,1,1) + xlb(:,n,1) * shdrv(:,1,n)
           dxdxib(:,1,2) = dxdxib(:,1,2) + xlb(:,n,1) * shdrv(:,2,n)
           dxdxib(:,1,3) = dxdxib(:,1,3) + xlb(:,n,1) * shdrv(:,3,n)
           dxdxib(:,2,1) = dxdxib(:,2,1) + xlb(:,n,2) * shdrv(:,1,n)
           dxdxib(:,2,2) = dxdxib(:,2,2) + xlb(:,n,2) * shdrv(:,2,n)
           dxdxib(:,2,3) = dxdxib(:,2,3) + xlb(:,n,2) * shdrv(:,3,n)
           dxdxib(:,3,1) = dxdxib(:,3,1) + xlb(:,n,3) * shdrv(:,1,n)
           dxdxib(:,3,2) = dxdxib(:,3,2) + xlb(:,n,3) * shdrv(:,2,n)
           dxdxib(:,3,3) = dxdxib(:,3,3) + xlb(:,n,3) * shdrv(:,3,n)
        enddo
c
c.... compute the normal to the boundary. This is achieved by taking
c     the cross product of two vectors in the plane of the 2-d
c     boundary face.
c
        v1 = xlb(:,2,:) - xlb(:,1,:)
        v2 = xlb(:,3,:) - xlb(:,1,:)

c
c.....The following are done in order to correct temp1..3
c     based on the results from compressible code.  This is done only
c     for wedges, depending on the bounary face.(tri or quad)
c
        if (lcsyst .eq. 4) then
           temp1 = dxdxib(:,2,1) * dxdxib(:,3,3) -
     &             dxdxib(:,2,3) * dxdxib(:,3,1)
           temp2 = dxdxib(:,3,1) * dxdxib(:,1,3) -
     &             dxdxib(:,3,3) * dxdxib(:,1,1)
           temp3 = dxdxib(:,1,1) * dxdxib(:,2,3) -
     &             dxdxib(:,1,3) * dxdxib(:,2,1)

        elseif (lcsyst .eq. 1) then
           temp1 = v1(:,2) * v2(:,3) - v2(:,2) * v1(:,3)
           temp2 = v2(:,1) * v1(:,3) - v1(:,1) * v2(:,3)
           temp3 = v1(:,1) * v2(:,2) - v2(:,1) * v1(:,2)
        else
           temp1 = - v1(:,2) * v2(:,3) + v2(:,2) * v1(:,3)
           temp2 = - v2(:,1) * v1(:,3) + v1(:,1) * v2(:,3)
           temp3 = - v1(:,1) * v2(:,2) + v2(:,1) * v1(:,2)
        endif
c
        temp       = one / sqrt ( temp1**2 + temp2**2 + temp3**2 )
        bnorm(:,1) = temp1 * temp
        bnorm(:,2) = temp2 * temp
        bnorm(:,3) = temp3 * temp
c

        if (lcsyst .eq. 3) then
           WdetJb     = (1 - Qwtb(lcsyst,intp)) / (four*temp)
        elseif (lcsyst .eq. 4) then
           WdetJb     = Qwtb(lcsyst,intp) / temp
        else
           WdetJb     = Qwtb(lcsyst,intp) / (four*temp)
        endif
c
c.... -------------------------->  Grad-V  <----------------------------
c
c.... compute grad-v for Navier-Stokes terms
c
        if (Navier .eq. 1) then
c
c.... compute the inverse of deformation gradient
c
          dxidxb(:,1,1) =   dxdxib(:,2,2) * dxdxib(:,3,3)
     &                    - dxdxib(:,3,2) * dxdxib(:,2,3)
          dxidxb(:,1,2) =   dxdxib(:,3,2) * dxdxib(:,1,3)
     &                    - dxdxib(:,1,2) * dxdxib(:,3,3)
          dxidxb(:,1,3) =   dxdxib(:,1,2) * dxdxib(:,2,3)
     &                    - dxdxib(:,1,3) * dxdxib(:,2,2)
          temp          = one / ( dxidxb(:,1,1) * dxdxib(:,1,1)
     &                          + dxidxb(:,1,2) * dxdxib(:,2,1)
     &                          + dxidxb(:,1,3) * dxdxib(:,3,1) )
          dxidxb(:,1,1) =  dxidxb(:,1,1) * temp
          dxidxb(:,1,2) =  dxidxb(:,1,2) * temp
          dxidxb(:,1,3) =  dxidxb(:,1,3) * temp
          dxidxb(:,2,1) = (dxdxib(:,2,3) * dxdxib(:,3,1)
     &                   - dxdxib(:,2,1) * dxdxib(:,3,3)) * temp
          dxidxb(:,2,2) = (dxdxib(:,1,1) * dxdxib(:,3,3)
     &                   - dxdxib(:,3,1) * dxdxib(:,1,3)) * temp
          dxidxb(:,2,3) = (dxdxib(:,2,1) * dxdxib(:,1,3)
     &                   - dxdxib(:,1,1) * dxdxib(:,2,3)) * temp
          dxidxb(:,3,1) = (dxdxib(:,2,1) * dxdxib(:,3,2)
     &                   - dxdxib(:,2,2) * dxdxib(:,3,1)) * temp
          dxidxb(:,3,2) = (dxdxib(:,3,1) * dxdxib(:,1,2)
     &                   - dxdxib(:,1,1) * dxdxib(:,3,2)) * temp
          dxidxb(:,3,3) = (dxdxib(:,1,1) * dxdxib(:,2,2)
     &                   - dxdxib(:,1,2) * dxdxib(:,2,1)) * temp
c
c.... compute local-grad-Y
c
c
          gradSl = zero
          isc=5+isclr
          do n = 1, nshl
            gradSl(:,1) = gradSl(:,1) + shdrv(:,1,n) * yl(:,n,isc)
            gradSl(:,2) = gradSl(:,2) + shdrv(:,2,n) * yl(:,n,isc)
            gradSl(:,3) = gradSl(:,3) + shdrv(:,3,n) * yl(:,n,isc)
          enddo
c
c.... convert local-grads to global-grads
c
          gradS(:,1) = dxidxb(:,1,1) * gradSl(:,1) +
     &                 dxidxb(:,2,1) * gradSl(:,2) +
     &                 dxidxb(:,3,1) * gradSl(:,3)

c
          gradS(:,2) = dxidxb(:,1,2) * gradSl(:,1) +
     &                 dxidxb(:,2,2) * gradSl(:,2) +
     &                 dxidxb(:,3,2) * gradSl(:,3)

          gradS(:,3) = dxidxb(:,1,3) * gradSl(:,1) +
     &                 dxidxb(:,2,3) * gradSl(:,2) +
     &                 dxidxb(:,3,3) * gradSl(:,3)
c
c.... end grad-T
c
        endif

        flux = diffus * ( gradS(:,1) * bnorm(:,1)
     &                  + gradS(:,2) * bnorm(:,2)
     &                  + gradS(:,3) * bnorm(:,3) )
c
c.... return
c
        return
        end


c---------------------------------------------------------------------
c
c     rotates the local nodal stiffnesses to the goblal frame
c
c---------------------------------------------------------------------
      subroutine rotatestiff(rKlocal, rotation,
     &                       rKglobal)

      include "common.h"

      dimension rKlocal(npro,nsd,nsd), rotation(npro,nsd,nsd),
     &          rKglobal(npro,nsd,nsd)

      dimension tempm(npro,nsd,nsd)

      tempm = zero
      do i = 1, 3
         do j = 1, 3
            do k = 1, 3
               tempm(:,i,j) = tempm(:,i,j)
     &              + rKlocal(:,i,k) * rotation(:,k,j)
            enddo
         enddo
      enddo

      rKglobal = zero
      do i = 1, 3
         do j = 1, 3
            do k = 1, 3
               rKglobal(:,i,j) = rKglobal(:,i,j)
     &              + rotation(:,k,i) * tempm(:,k,j)
            enddo
         enddo
      enddo

      return
      end