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

      subroutine e3StsLhs( xl,  lStsVec )
c-----------------------------------------------------------------------
c
c  compute the terms needed for the left hand side matrices
c  needed for the conservative projection
c
c-----------------------------------------------------------------------
      use     stats
      include "common.h"

      integer i
      real*8  lDir(npro,nshl,3), lStsVec(npro,nshl,nResDims),
     &        xl(npro,nenl,3)

      call e3StsDir( xl,  lDir )

      do i = 1, nshl
         lStsVec(:,i,1) = lDir(:,i,1) * lDir(:,i,1)
         lStsVec(:,i,2) = lDir(:,i,2) * lDir(:,i,2)
         lStsVec(:,i,3) = lDir(:,i,3) * lDir(:,i,3)

         lStsVec(:,i,4) = lDir(:,i,1) * lDir(:,i,2)
         lStsVec(:,i,5) = lDir(:,i,2) * lDir(:,i,3)
         lStsVec(:,i,6) = lDir(:,i,3) * lDir(:,i,1)

         lStsVec(:,i,7) = 0.0
         lStsVec(:,i,8) = 0.0
         lStsVec(:,i,9) = 0.0
         lStsVec(:,i,10) = 0.0
         lStsVec(:,i,11) = 0.0
      enddo

      return
      end


      subroutine e3StsRes( xl, rl, lStsVec )
c-----------------------------------------------------------------------
c
c  compute the residual terms for the consistent projection
c
c-----------------------------------------------------------------------
      use     stats
      include "common.h"

      real*8  xl(npro,nenl,3),  rl(npro,nshl,ndof)
      real*8  lDir(npro,nshl,3), lStsVec(npro,nshl,nResDims)

      call e3StsDir( xl,  lDir )

      do i = 1, nshl
         lStsVec(:,i,1) = lDir(:,i,1) * rl(:,i,4)
         lStsVec(:,i,2) = lDir(:,i,2) * rl(:,i,4)
         lStsVec(:,i,3) = lDir(:,i,3) * rl(:,i,4)

         lStsVec(:,i,4) = lDir(:,i,1) * rl(:,i,1)
         lStsVec(:,i,5) = lDir(:,i,2) * rl(:,i,2)
         lStsVec(:,i,6) = lDir(:,i,3) * rl(:,i,3)

         lStsVec(:,i,7) = lDir(:,i,1) * rl(:,i,2)
     &                  + lDir(:,i,2) * rl(:,i,1)
         lStsVec(:,i,8) = lDir(:,i,2) * rl(:,i,3)
     &                  + lDir(:,i,3) * rl(:,i,2)
         lStsVec(:,i,9) = lDir(:,i,3) * rl(:,i,1)
     &        + lDir(:,i,1) * rl(:,i,3)
         lStsVec(:,i,10) = 0
         lStsVec(:,i,11) = 0
      enddo


      return
      end

      subroutine e3StsDir( xl,  lDir )
c-----------------------------------------------------------------------
c
c  compute the normal to each of the nodes
c
c-----------------------------------------------------------------------
      include "common.h"

      real*8  xl(npro,nenl,3), lDir(npro,nshl,3)
      integer e

c
c.... linear tets
c
      if (nshl .eq. 4 ) then
         fct = 1.d0 / 6.d0
         do e = 1, npro

            x12         = xl(e,2,1) - xl(e,1,1)
            x13         = xl(e,3,1) - xl(e,1,1)
            x14         = xl(e,4,1) - xl(e,1,1)
            x23         = xl(e,3,1) - xl(e,2,1)
            x24         = xl(e,4,1) - xl(e,2,1)
            x34         = xl(e,4,1) - xl(e,3,1)

            y12         = xl(e,2,2) - xl(e,1,2)
            y13         = xl(e,3,2) - xl(e,1,2)
            y14         = xl(e,4,2) - xl(e,1,2)
            y23         = xl(e,3,2) - xl(e,2,2)
            y24         = xl(e,4,2) - xl(e,2,2)
            y34         = xl(e,4,2) - xl(e,3,2)

            z12         = xl(e,2,3) - xl(e,1,3)
            z13         = xl(e,3,3) - xl(e,1,3)
            z14         = xl(e,4,3) - xl(e,1,3)
            z23         = xl(e,3,3) - xl(e,2,3)
            z24         = xl(e,4,3) - xl(e,2,3)
            z34         = xl(e,4,3) - xl(e,3,3)

c
c.. The calculation of the direction of a vertex is based on the average of
c.. the normals of the neighbor faces(3); And the calculation of the direction
c.. of the edge is based on the neighbor faces(2);
c
	    lDir(e,1,1) = fct * (y14*(z12 - z13) + y12*(z13 - z14)
     &                         + y13*(-z12 + z14))
	    lDir(e,1,2) = fct * ( x14*(-z12 + z13) + x13*(z12 - z14)
     &                          + x12*(-z13 + z14))
	    lDir(e,1,3) = fct * ( x14*(y12 - y13)
     &                          + x12*(y13 - y14) + x13*(-y12 + y14))
c
	    lDir(e,2,1) = fct * (-(y13*z12) + y14*z12 + y12*z13
     &                          - y12*z14 + y24*z23 - y23*z24)
            lDir(e,2,2) = fct * (x13*z12 - x14*z12-x12*z13 + x12*z14
     &                          - x24*z23 + x23*z24 )
            lDir(e,2,3) = fct * (-(x13*y12) + x14*y12 + x12*y13
     &                         - x12*y14 + x24*y23 - x23*y24)
c
            lDir(e,3,1) = fct * (y12*z13 - y14*z13 + y13*(-z12 + z14)
     &                         + y24*z23 - y23*z24)
            lDir(e,3,2) = fct * (-(x12*z13) + x14*z13 + x13*(z12 - z14)
     &                         - x24*z23 + x23*z24)
            lDir(e,3,3) = fct * (x12*y13 - x14*y13 + x13*(-y12 + y14)
     &                         + x24*y23 - x23*y24)
c
            lDir(e,4,1) = fct * (y14*(z12 - z13) - y12*z14 + y13*z14
     &                         + y24*z23 - y23*z24)
            lDir(e,4,2) = fct * (x14*(-z12 + z13) + x12*z14 - x13*z14
     &                         - x24*z23 + x23*z24)
            lDir(e,4,3) = fct * (x14*(y12 - y13) - x12*y14 + x13*y14
     &                         + x24*y23 - x23*y24)
c
         enddo
c
c.... quadratic tets
c
      else if (nshl .eq. 10 ) then
         fct = 1.d0 / 6.d0
         do e = 1, npro

            x12         = xl(e,2,1) - xl(e,1,1)
            x13         = xl(e,3,1) - xl(e,1,1)
            x14         = xl(e,4,1) - xl(e,1,1)
            x23         = xl(e,3,1) - xl(e,2,1)
            x24         = xl(e,4,1) - xl(e,2,1)
            x34         = xl(e,4,1) - xl(e,3,1)

            y12         = xl(e,2,2) - xl(e,1,2)
            y13         = xl(e,3,2) - xl(e,1,2)
            y14         = xl(e,4,2) - xl(e,1,2)
            y23         = xl(e,3,2) - xl(e,2,2)
            y24         = xl(e,4,2) - xl(e,2,2)
            y34         = xl(e,4,2) - xl(e,3,2)

            z12         = xl(e,2,3) - xl(e,1,3)
            z13         = xl(e,3,3) - xl(e,1,3)
            z14         = xl(e,4,3) - xl(e,1,3)
            z23         = xl(e,3,3) - xl(e,2,3)
            z24         = xl(e,4,3) - xl(e,2,3)
            z34         = xl(e,4,3) - xl(e,3,3)

c
c.... vertex modes
c
	    lDir(e,1,1) = fct * (y14*(z12 - z13) + y12*(z13 - z14)
     &                         + y13*(-z12 + z14))
	    lDir(e,1,2) = fct * ( x14*(-z12 + z13) + x13*(z12 - z14)
     &                          + x12*(-z13 + z14))
	    lDir(e,1,3) = fct * ( x14*(y12 - y13)
     &                          + x12*(y13 - y14) + x13*(-y12 + y14))
c
	    lDir(e,2,1) = fct * (-(y13*z12) + y14*z12 + y12*z13
     &                          - y12*z14 + y24*z23 - y23*z24)
            lDir(e,2,2) = fct * (x13*z12 - x14*z12-x12*z13 + x12*z14
     &                          - x24*z23 + x23*z24 )
            lDir(e,2,3) = fct * (-(x13*y12) + x14*y12 + x12*y13
     &                         - x12*y14 + x24*y23 - x23*y24)
c
            lDir(e,3,1) = fct * (y12*z13 - y14*z13 + y13*(-z12 + z14)
     &                         + y24*z23 - y23*z24)
            lDir(e,3,2) = fct * (-(x12*z13) + x14*z13 + x13*(z12 - z14)
     &                         - x24*z23 + x23*z24)
            lDir(e,3,3) = fct * (x12*y13 - x14*y13 + x13*(-y12 + y14)
     &                         + x24*y23 - x23*y24)
c
            lDir(e,4,1) = fct * (y14*(z12 - z13) - y12*z14 + y13*z14
     &                         + y24*z23 - y23*z24)
            lDir(e,4,2) = fct * (x14*(-z12 + z13) + x12*z14 - x13*z14
     &                         - x24*z23 + x23*z24)
            lDir(e,4,3) = fct * (x14*(y12 - y13) - x12*y14 + x13*y14
     &                         + x24*y23 - x23*y24)
c
c.... edge modes (quadratic)
c
            lDir(e,5,1) = pt25*(-(y13*z12) + y14*z12 + y12*(z13-z14))
            lDir(e,5,2) = pt25*(x13*z12 - x14*z12 + x12*(-z13 + z14))
            lDir(e,5,3) = pt25*(-(x13*y12) + x14*y12 + x12*(y13 - y14))

            lDir(e,6,1) = pt25*(-(y13*z12) + y12*z13 + y24*z23-y23*z24)
            lDir(e,6,2) = pt25*(x13*z12 - x12*z13 - x24*z23 + x23*z24)
            lDir(e,6,3) = pt25*(-(x13*y12) + x12*y13 + x24*y23-x23*y24)

            lDir(e,7,1) = pt25*((y12 - y14)*z13 + y13*(-z12 + z14))
            lDir(e,7,2) = pt25*((-x12 + x14)*z13 + x13*(z12 - z14))
            lDir(e,7,3) = pt25*((x12 - x14)*y13 + x13*(-y12 + y14))

            lDir(e,8,1) = pt25*(y14*(z12 - z13) + (-y12 + y13)*z14)
            lDir(e,8,2) = pt25*(x14*(-z12 + z13) + (x12 - x13)*z14)
            lDir(e,8,3) = pt25*(x14*(y12 - y13) + (-x12 + x13)*y14)

            lDir(e,9,1) = pt25*(y14*z12 - y12*z14 + y24*z23 - y23*z24)
            lDir(e,9,2) = pt25*(-(x14*z12) + x12*z14 - x24*z23+x23*z24)
            lDir(e,9,3) = pt25*(x14*y12 - x12*y14 + x24*y23 - x23*y24)

            lDir(e,10,1) = pt25*(-(y14*z13) + y13*z14+y24*z23-y23*z24)
            lDir(e,10,2) = pt25*(x14*z13 - x13*z14-x24*z23 + x23*z24)
            lDir(e,10,3) = pt25*(-(x14*y13) + x13*y14+x24*y23-x23*y24)


         enddo
c
c.... cubic tets
c
      else if (nshl .eq. 20 ) then
         fct = 1.d0 / 6.d0
         do e = 1, npro

            x12         = xl(e,2,1) - xl(e,1,1)
            x13         = xl(e,3,1) - xl(e,1,1)
            x14         = xl(e,4,1) - xl(e,1,1)
            x23         = xl(e,3,1) - xl(e,2,1)
            x24         = xl(e,4,1) - xl(e,2,1)
c$$$            x34         = xl(e,4,1) - xl(e,3,1)

            y12         = xl(e,2,2) - xl(e,1,2)
            y13         = xl(e,3,2) - xl(e,1,2)
            y14         = xl(e,4,2) - xl(e,1,2)
            y23         = xl(e,3,2) - xl(e,2,2)
            y24         = xl(e,4,2) - xl(e,2,2)
c$$$            y34         = xl(e,4,2) - xl(e,3,2)

            z12         = xl(e,2,3) - xl(e,1,3)
            z13         = xl(e,3,3) - xl(e,1,3)
            z14         = xl(e,4,3) - xl(e,1,3)
            z23         = xl(e,3,3) - xl(e,2,3)
            z24         = xl(e,4,3) - xl(e,2,3)
c$$$            z34         = xl(e,4,3) - xl(e,3,3)

c
c.... vertex modes
c
	    lDir(e,1,1) = fct * (y14*(z12 - z13) + y12*(z13 - z14)
     &                         + y13*(-z12 + z14))
	    lDir(e,1,2) = fct * ( x14*(-z12 + z13) + x13*(z12 - z14)
     &                          + x12*(-z13 + z14))
	    lDir(e,1,3) = fct * ( x14*(y12 - y13)
     &                          + x12*(y13 - y14) + x13*(-y12 + y14))
c
	    lDir(e,2,1) = fct * (-(y13*z12) + y14*z12 + y12*z13
     &                          - y12*z14 + y24*z23 - y23*z24)
            lDir(e,2,2) = fct * (x13*z12 - x14*z12-x12*z13 + x12*z14
     &                          - x24*z23 + x23*z24 )
            lDir(e,2,3) = fct * (-(x13*y12) + x14*y12 + x12*y13
     &                         - x12*y14 + x24*y23 - x23*y24)
c
            lDir(e,3,1) = fct * (y12*z13 - y14*z13 + y13*(-z12 + z14)
     &                         + y24*z23 - y23*z24)
            lDir(e,3,2) = fct * (-(x12*z13) + x14*z13 + x13*(z12 - z14)
     &                         - x24*z23 + x23*z24)
            lDir(e,3,3) = fct * (x12*y13 - x14*y13 + x13*(-y12 + y14)
     &                         + x24*y23 - x23*y24)
c
            lDir(e,4,1) = fct * (y14*(z12 - z13) - y12*z14 + y13*z14
     &                         + y24*z23 - y23*z24)
            lDir(e,4,2) = fct * (x14*(-z12 + z13) + x12*z14 - x13*z14
     &                         - x24*z23 + x23*z24)
            lDir(e,4,3) = fct * (x14*(y12 - y13) - x12*y14 + x13*y14
     &                         + x24*y23 - x23*y24)
c
c.... edge modes (quadratic and cubic)
c
            lDir(e,5,1) = pt25*(-(y13*z12) + y14*z12 + y12*(z13-z14))
            lDir(e,5,2) = pt25*(x13*z12 - x14*z12 + x12*(-z13 + z14))
            lDir(e,5,3) = pt25*(-(x13*y12) + x14*y12 + x12*(y13 - y14))
            lDir(e,6,1) = lDir(e,5,1)
            lDir(e,6,2) = lDir(e,5,2)
            lDir(e,6,3) = lDir(e,5,3)

            lDir(e,7,1) = pt25*(-(y13*z12) + y12*z13 + y24*z23-y23*z24)
            lDir(e,7,2) = pt25*(x13*z12 - x12*z13 - x24*z23 + x23*z24)
            lDir(e,7,3) = pt25*(-(x13*y12) + x12*y13 + x24*y23-x23*y24)
            lDir(e,8,1) = lDir(e,7,1)
            lDir(e,8,2) = lDir(e,7,2)
            lDir(e,8,3) = lDir(e,7,3)

            lDir(e,9,1) = pt25*((y12 - y14)*z13 + y13*(-z12 + z14))
            lDir(e,9,2) = pt25*((-x12 + x14)*z13 + x13*(z12 - z14))
            lDir(e,9,3) = pt25*((x12 - x14)*y13 + x13*(-y12 + y14))
            lDir(e,10,1) = lDir(e,9,1)
            lDir(e,10,2) = lDir(e,9,2)
            lDir(e,10,3) = lDir(e,9,3)

            lDir(e,11,1) = pt25*(y14*(z12 - z13) + (-y12 + y13)*z14)
            lDir(e,11,2) = pt25*(x14*(-z12 + z13) + (x12 - x13)*z14)
            lDir(e,11,3) = pt25*(x14*(y12 - y13) + (-x12 + x13)*y14)
            lDir(e,12,1) = lDir(e,11,1)
            lDir(e,12,2) = lDir(e,11,2)
            lDir(e,12,3) = lDir(e,11,3)


            lDir(e,13,1) = pt25*(y14*z12 - y12*z14 + y24*z23 - y23*z24)
            lDir(e,13,2) = pt25*(-(x14*z12) + x12*z14-x24*z23+x23*z24)
            lDir(e,13,3) = pt25*(x14*y12 - x12*y14 + x24*y23 - x23*y24)
            lDir(e,14,1) = lDir(e,13,1)
            lDir(e,14,2) = lDir(e,13,2)
            lDir(e,14,3) = lDir(e,13,3)

            lDir(e,15,1) = pt25*(-(y14*z13) + y13*z14+y24*z23-y23*z24)
            lDir(e,15,2) = pt25*(x14*z13 - x13*z14-x24*z23 + x23*z24)
            lDir(e,15,3) = pt25*(-(x14*y13) + x13*y14+x24*y23-x23*y24)
            lDir(e,16,1) = lDir(e,15,1)
            lDir(e,16,2) = lDir(e,15,2)
            lDir(e,16,3) = lDir(e,15,3)
c
c.... face modes (cubic)
c
            lDir(e,17,1) = pt5*(-(y13*z12) + y12*z13)
            lDir(e,17,2) = pt5*(x13*z12 - x12*z13)
            lDir(e,17,3) = pt5*(-(x13*y12) + x12*y13)

            lDir(e,18,1) = pt5*(y14*z12 - y12*z14)
            lDir(e,18,2) = pt5*(-(x14*z12) + x12*z14)
            lDir(e,18,3) = pt5*(x14*y12 - x12*y14)

            lDir(e,19,1) = pt5*(y24*z23 - y23*z24)
            lDir(e,19,2) = pt5*(-(x24*z23) + x23*z24)
            lDir(e,19,3) = pt5*(x24*y23 - x23*y24)

            lDir(e,20,1) = pt5*(-(y14*z13) + y13*z14)
            lDir(e,20,2) = pt5*(x14*z13 - x13*z14)
            lDir(e,20,3) = pt5*(-(x14*y13) + x13*y14)

         enddo
c
c.... hexes
c
      else if (nenl .eq. 8) then
         fct = 1.d0 / 12.d0
         do e = 1, npro
	    x13		= xl(e,1,1) - xl(e,3,1)
	    x16		= xl(e,1,1) - xl(e,6,1)
	    x18		= xl(e,1,1) - xl(e,8,1)
	    x24		= xl(e,2,1) - xl(e,4,1)
	    x25		= xl(e,2,1) - xl(e,5,1)
	    x27		= xl(e,2,1) - xl(e,7,1)
	    x36		= xl(e,3,1) - xl(e,6,1)
	    x38		= xl(e,3,1) - xl(e,8,1)
	    x45		= xl(e,4,1) - xl(e,5,1)
	    x47		= xl(e,4,1) - xl(e,7,1)
	    x57		= xl(e,5,1) - xl(e,7,1)
	    x68		= xl(e,6,1) - xl(e,8,1)
c
	    y13		= xl(e,1,2) - xl(e,3,2)
	    y16		= xl(e,1,2) - xl(e,6,2)
	    y18		= xl(e,1,2) - xl(e,8,2)
	    y24		= xl(e,2,2) - xl(e,4,2)
	    y25		= xl(e,2,2) - xl(e,5,2)
	    y27		= xl(e,2,2) - xl(e,7,2)
	    y36		= xl(e,3,2) - xl(e,6,2)
	    y38		= xl(e,3,2) - xl(e,8,2)
	    y45		= xl(e,4,2) - xl(e,5,2)
	    y47		= xl(e,4,2) - xl(e,7,2)
	    y57		= xl(e,5,2) - xl(e,7,2)
	    y68		= xl(e,6,2) - xl(e,8,2)
c
	    z13		= xl(e,1,3) - xl(e,3,3)
	    z16		= xl(e,1,3) - xl(e,6,3)
	    z18		= xl(e,1,3) - xl(e,8,3)
	    z24		= xl(e,2,3) - xl(e,4,3)
	    z25		= xl(e,2,3) - xl(e,5,3)
	    z27		= xl(e,2,3) - xl(e,7,3)
	    z36		= xl(e,3,3) - xl(e,6,3)
	    z38		= xl(e,3,3) - xl(e,8,3)
	    z45		= xl(e,4,3) - xl(e,5,3)
	    z47		= xl(e,4,3) - xl(e,7,3)
	    z57		= xl(e,5,3) - xl(e,7,3)
	    z68		= xl(e,6,3) - xl(e,8,3)
c
	x31= -x13
	x61= -x16
	x81= -x18
	x42= -x24
	x52= -x25
	x72= -x27
	x63= -x36
	x83= -x38
	x54= -x45
	x74= -x47
	x75= -x57
	x86= -x68
	y31= -y13
	y61= -y16
	y81= -y18
	y42= -y24
	y52= -y25
	y72= -y27
	y63= -y36
	y83= -y38
	y54= -y45
	y74= -y47
	y75= -y57
	y86= -y68
	z31= -z13
	z61= -z16
	z81= -z18
	z42= -z24
	z52= -z25
	z72= -z27
	z63= -z36
	z83= -z38
	z54= -z45
	z74= -z47
	z75= -z57
	z86= -z68

	    lDir(e,1,1) = fct * (-y24 * z45 + y36 * z24 - y68 * z45
     1				+ z24 * y45 - z36 * y24 + z68 * y45 )
	    lDir(e,2,1) = fct * (-y16 * z63 + y54 * z16 - y47 * z63
     1				+ z16 * y63 - z54 * y16 + z47 * y63 )
	    lDir(e,3,1) = fct * (-y42 * z27 + y18 * z42 - y86 * z27
     1				+ z42 * y27 - z18 * y42 + z86 * y27 )
	    lDir(e,4,1) = fct * (-y38 * z81 + y72 * z38 - y25 * z81
     1				+ z38 * y81 - z72 * y38 + z25 * y81 )
	    lDir(e,5,1) = fct * (-y61 * z18 + y27 * z61 - y74 * z18
     1				+ z61 * y18 - z27 * y61 + z74 * y18 )
	    lDir(e,6,1) = fct * (-y57 * z72 + y81 * z57 - y13 * z72
     1				+ z57 * y72 - z81 * y57 + z13 * y72 )
	    lDir(e,7,1) = fct * (-y83 * z36 + y45 * z83 - y52 * z36
     1				+ z83 * y36 - z45 * y83 + z52 * y36 )
	    lDir(e,8,1) = fct * (-y75 * z54 + y63 * z75 - y31 * z54
     1				+ z75 * y54 - z63 * y75 + z31 * y54 )
c
	    lDir(e,1,2) = fct * (-z24 * x45 + z36 * x24 - z68 * x45
     1				+ x24 * z45 - x36 * z24 + x68 * z45 )
	    lDir(e,2,2) = fct * (-z16 * x63 + z54 * x16 - z47 * x63
     1				+ x16 * z63 - x54 * z16 + x47 * z63 )
	    lDir(e,3,2) = fct * (-z42 * x27 + z18 * x42 - z86 * x27
     1				+ x42 * z27 - x18 * z42 + x86 * z27 )
	    lDir(e,4,2) = fct * (-z38 * x81 + z72 * x38 - z25 * x81
     1				+ x38 * z81 - x72 * z38 + x25 * z81 )
	    lDir(e,5,2) = fct * (-z61 * x18 + z27 * x61 - z74 * x18
     1				+ x61 * z18 - x27 * z61 + x74 * z18 )
	    lDir(e,6,2) = fct * (-z57 * x72 + z81 * x57 - z13 * x72
     1				+ x57 * z72 - x81 * z57 + x13 * z72 )
	    lDir(e,7,2) = fct * (-z83 * x36 + z45 * x83 - z52 * x36
     1				+ x83 * z36 - x45 * z83 + x52 * z36 )
	    lDir(e,8,2) = fct * (-z75 * x54 + z63 * x75 - z31 * x54
     1				+ x75 * z54 - x63 * z75 + x31 * z54 )
c
	    lDir(e,1,3) = fct * (-x24 * y45 + x36 * y24 - x68 * y45
     1				+ y24 * x45 - y36 * x24 + y68 * x45 )
	    lDir(e,2,3) = fct * (-x16 * y63 + x54 * y16 - x47 * y63
     1				+ y16 * x63 - y54 * x16 + y47 * x63 )
	    lDir(e,3,3) = fct * (-x42 * y27 + x18 * y42 - x86 * y27
     1				+ y42 * x27 - y18 * x42 + y86 * x27 )
	    lDir(e,4,3) = fct * (-x38 * y81 + x72 * y38 - x25 * y81
     1				+ y38 * x81 - y72 * x38 + y25 * x81 )
	    lDir(e,5,3) = fct * (-x61 * y18 + x27 * y61 - x74 * y18
     1				+ y61 * x18 - y27 * x61 + y74 * x18 )
	    lDir(e,6,3) = fct * (-x57 * y72 + x81 * y57 - x13 * y72
     1				+ y57 * x72 - y81 * x57 + y13 * x72 )
	    lDir(e,7,3) = fct * (-x83 * y36 + x45 * y83 - x52 * y36
     1				+ y83 * x36 - y45 * x83 + y52 * x36 )
	    lDir(e,8,3) = fct * (-x75 * y54 + x63 * y75 - x31 * y54
     1				+ y75 * x54 - y63 * x75 + y31 * x54 )
c
         enddo
      else
         write(*,*) 'Error in e3sts: elt type not impl.'
         stop
      endif

      return
      end