s4deri3.F File Reference

#include "implicit_f.inc"
#include "comlock.inc"
#include "mvsiz_p.inc"
#include "param_c.inc"
#include "scr17_c.inc"
#include "vectorize.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	s4deri3 (off, det, ngl, deltax, mxt, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, px1, px2, px3, px4, py1, py2, py3, py4, pz1, pz2, pz3, pz4, rx, ry, rz, sx, sy, sz, tx, ty, tz, sav, offg, nel, pm, voldp, ismstr, iformdt, jlag)

Function/Subroutine Documentation

s4deri3()

subroutine s4deri3	(		off,
			det,
		integer, dimension(*)	ngl,
			deltax,
		integer, dimension(mvsiz)	mxt,
		double precision, dimension(*)	x1,
		double precision, dimension(*)	x2,
		double precision, dimension(*)	x3,
		double precision, dimension(*)	x4,
		double precision, dimension(*)	y1,
		double precision, dimension(*)	y2,
		double precision, dimension(*)	y3,
		double precision, dimension(*)	y4,
		double precision, dimension(*)	z1,
		double precision, dimension(*)	z2,
		double precision, dimension(*)	z3,
		double precision, dimension(*)	z4,
			px1,
			px2,
			px3,
			px4,
			py1,
			py2,
			py3,
			py4,
			pz1,
			pz2,
			pz3,
			pz4,
			rx,
			ry,
			rz,
			sx,
			sy,
			sz,
			tx,
			ty,
			tz,
		double precision, dimension(nel,9)	sav,
			offg,
		integer	nel,
			pm,
		double precision, dimension(nel)	voldp,
		integer, intent(in)	ismstr,
		integer, intent(in)	iformdt,
		integer, intent(in)	jlag )

Definition at line 30 of file s4deri3.F.

C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.inc"
#include      "comlock.inc"
C-----------------------------------------------
C   G l o b a l   P a r a m e t e r s
C-----------------------------------------------
#include      "mvsiz_p.inc"
C-----------------------------------------------
C   C o m m o n   B l o c k s
C-----------------------------------------------
#include      "param_c.inc"
#include      "scr17_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER, INTENT(IN) :: JLAG
      INTEGER, INTENT(IN) :: ISMSTR
      INTEGER, INTENT(IN) :: IFORMDT
      INTEGER NEL, MXT(MVSIZ)
      double precision
     .   x1(*), x2(*), x3(*), x4(*),
     .   y1(*), y2(*), y3(*), y4(*),
     .   z1(*), z2(*), z3(*), z4(*), sav(nel,9),voldp(nel)
     
      my_real
     .   off(*), det(*), deltax(*), pm(npropm,*),
     .   px1(*), px2(*), px3(*), px4(*),  
     .   py1(*), py2(*), py3(*), py4(*),  
     .   pz1(*), pz2(*), pz3(*), pz4(*),offg(*),
     .   rx(*), ry(*), rz(*), sx(*), sy(*), sz(*),tx(*), ty(*), tz(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER NGL(*), I,J
      INTEGER NNEGA,INDEX(MVSIZ)
      double precision
     .   x41, y41, z41, x42, y42, z42, x43, y43, z43,b1dp,c1dp,d1dp
      my_real
     .   a1, a2, a3, a4, d, areamx2, 
     .   b1(mvsiz), b2(mvsiz), b3(mvsiz), b4(mvsiz), 
     .   c1(mvsiz), c2(mvsiz), c3(mvsiz), c4(mvsiz),
     .   d1(mvsiz), d2(mvsiz), d3(mvsiz), d4(mvsiz),
     .   pxx, pyy, pzz, pxy, pyz, pxz, gfac, aa, bb, p, ld
C-----------------------------------------------
C
      DO i=1,nel
       x43 = x4(i) - x3(i)
       y43 = y4(i) - y3(i)
       z43 = z4(i) - z3(i)
       x41 = x4(i) - x1(i)
       y41 = y4(i) - y1(i)
       z41 = z4(i) - z1(i)
       x42 = x4(i) - x2(i)
       y42 = y4(i) - y2(i)
       z42 = z4(i) - z2(i)
C
       rx(i) =  -x41
       ry(i) =  -y41
       rz(i) =  -z41
       sx(i) =  -x42
       sy(i) =  -y42
       sz(i) =  -z42
C
       tx(i) =  -x43
       ty(i) =  -y43
       tz(i) =  -z43
C
       b1dp  =  y43*z42 - y42*z43
       b1(i) =  b1dp
       b2(i) =  y41*z43 - y43*z41
       b3(i) =  y42*z41 - y41*z42
       b4(i) =  -(b1(i) + b2(i) + b3(i))
C
       c1dp  =  z43*x42 - z42*x43
       c1(i) =  c1dp
       c2(i) =  z41*x43 - z43*x41
       c3(i) =  z42*x41 - z41*x42
       c4(i) =  -(c1(i) + c2(i) + c3(i))
C
       d1dp  =  x43*y42 - x42*y43
       d1(i) =  d1dp
       d2(i) =  x41*y43 - x43*y41
       d3(i) =  x42*y41 - x41*y42
       d4(i) =  -(d1(i) + d2(i) + d3(i))
C
       voldp(i) = (x41*b1dp + y41*c1dp + z41*d1dp)*one_over_6
       det(i) = voldp(i)
C
      ENDDO
C
C
      CALL schkjabt3(
     1   off,     det,     ngl,     offg,
     2   nnega,   index,   nel,     ismstr,
     3   jlag)
      IF (nnega>0) THEN
        IF (ismstr==10.OR.ismstr==12) THEN
#include "vectorize.inc"
         DO j=1,nnega
          i = index(j)
            x1(i)=sav(i,1)
            y1(i)=sav(i,4)
            z1(i)=sav(i,7)
            x2(i)=sav(i,2)
            y2(i)=sav(i,5)
            z2(i)=sav(i,8)
            x3(i)=sav(i,3)
            y3(i)=sav(i,6)
            z3(i)=sav(i,9)
            x4(i)=zero
            y4(i)=zero
            z4(i)=zero
         ENDDO
        ELSE
#include "vectorize.inc"
         DO j=1,nnega
          i = index(j)
          x1(i)=sav(i,1)
          y1(i)=sav(i,2)
          z1(i)=sav(i,3)
          x2(i)=sav(i,4)
          y2(i)=sav(i,5)
          z2(i)=sav(i,6)
          x3(i)=sav(i,7)
          y3(i)=sav(i,8)
          z3(i)=sav(i,9)
          x4(i)=zero
          y4(i)=zero
          z4(i)=zero
         ENDDO
          END IF
#include "vectorize.inc"
         DO j=1,nnega
          i = index(j)
          x43 = x4(i) - x3(i)
          y43 = y4(i) - y3(i)
          z43 = z4(i) - z3(i)
          x41 = x4(i) - x1(i)
          y41 = y4(i) - y1(i)
          z41 = z4(i) - z1(i)
          x42 = x4(i) - x2(i)
          y42 = y4(i) - y2(i)
          z42 = z4(i) - z2(i)
C
          rx(i) =  -x41
          ry(i) =  -y41
          rz(i) =  -z41
          sx(i) =  -x42
          sy(i) =  -y42
          sz(i) =  -z42
          tx(i) =  -x43
          ty(i) =  -y43
          tz(i) =  -z43
C
C
          b1dp  =  y43*z42 - y42*z43
          b1(i) =  b1dp
          b2(i) =  y41*z43 - y43*z41
          b3(i) =  y42*z41 - y41*z42
          b4(i) =  -(b1(i) + b2(i) + b3(i))
C
          c1dp  =  z43*x42 - z42*x43
          c1(i) =  c1dp
          c2(i) =  z41*x43 - z43*x41
          c3(i) =  z42*x41 - z41*x42
          c4(i) =  -(c1(i) + c2(i) + c3(i))
C
          d1dp  =  x43*y42 - x42*y43
          d1(i) =  d1dp
          d2(i) =  x41*y43 - x43*y41
          d3(i) =  x42*y41 - x41*y42
          d4(i) =  -(d1(i) + d2(i) + d3(i))
C
          voldp(i) = (x41*b1dp + y41*c1dp + z41*d1dp)*one_over_6
          det(i) = voldp(i)
          offg(i) = two
         ENDDO
       END IF
C
      DO i=1,nel
        d = one/det(i)/six
        px1(i)=-b1(i)*d
        py1(i)=-c1(i)*d
        pz1(i)=-d1(i)*d
        px2(i)=-b2(i)*d
        py2(i)=-c2(i)*d
        pz2(i)=-d2(i)*d
        px3(i)=-b3(i)*d
        py3(i)=-c3(i)*d
        pz3(i)=-d3(i)*d
        px4(i)=-b4(i)*d
        py4(i)=-c4(i)*d
        pz4(i)=-d4(i)*d
       END DO
 
       IF(idt1sol==0)THEN
 
         DO i=1,nel
             d = max(px1(i)*px1(i)+py1(i)*py1(i)+pz1(i)*pz1(i),
     .         px2(i)*px2(i)+py2(i)*py2(i)+pz2(i)*pz2(i),
     .       px3(i)*px3(i)+py3(i)*py3(i)+pz3(i)*pz3(i),
     .       px4(i)*px4(i)+py4(i)*py4(i)+pz4(i)*pz4(i))
           deltax(i) = one / sqrt(d)
         END DO
 
       ELSEIF(iformdt==0)THEN
         DO i=1,nel
             d = px1(i)*px1(i)+py1(i)*py1(i)+pz1(i)*pz1(i)
     .       + px2(i)*px2(i)+py2(i)*py2(i)+pz2(i)*pz2(i)
     .       + px3(i)*px3(i)+py3(i)*py3(i)+pz3(i)*pz3(i)
     .       + px4(i)*px4(i)+py4(i)*py4(i)+pz4(i)*pz4(i)
           deltax(i) = one / sqrt(d)
         END DO
 
       ELSEIF(iformdt==1)THEN
 
         gfac=pm(105,mxt(1))
         ld  =two*sqrt(max(one-gfac,zero))+one
         DO i=1,nel
           pxx=px1(i)*px1(i)+px2(i)*px2(i)+px3(i)*px3(i)+px4(i)*px4(i)
           pyy=py1(i)*py1(i)+py2(i)*py2(i)+py3(i)*py3(i)+py4(i)*py4(i)
           pzz=pz1(i)*pz1(i)+pz2(i)*pz2(i)+pz3(i)*pz3(i)+pz4(i)*pz4(i)
           pxy=px1(i)*py1(i)+px2(i)*py2(i)+px3(i)*py3(i)+px4(i)*py4(i)
           pxz=px1(i)*pz1(i)+px2(i)*pz2(i)+px3(i)*pz3(i)+px4(i)*pz4(i)
           pyz=py1(i)*pz1(i)+py2(i)*pz2(i)+py3(i)*pz3(i)+py4(i)*pz4(i)
C
           aa = -(pxx+pyy+pzz)
           bb =  (pxx*pyy+pxx*pzz+pyy*pzz-pxy**2-pxz**2-pyz**2) 
           p  = bb-third*aa*aa
           d  = two*sqrt(third*max(-p,zero))-third*aa
C
             d=ld*d
C
           deltax(i) = one / sqrt(d)
         END DO
 
       ELSEIF(iformdt==2)THEN
 
         gfac=pm(105,mxt(1))
         DO i=1,nel
           pxx=px1(i)*px1(i)+px2(i)*px2(i)+px3(i)*px3(i)+px4(i)*px4(i)
           pyy=py1(i)*py1(i)+py2(i)*py2(i)+py3(i)*py3(i)+py4(i)*py4(i)
           pzz=pz1(i)*pz1(i)+pz2(i)*pz2(i)+pz3(i)*pz3(i)+pz4(i)*pz4(i)
           pxy=px1(i)*py1(i)+px2(i)*py2(i)+px3(i)*py3(i)+px4(i)*py4(i)
           pxz=px1(i)*pz1(i)+px2(i)*pz2(i)+px3(i)*pz3(i)+px4(i)*pz4(i)
           pyz=py1(i)*pz1(i)+py2(i)*pz2(i)+py3(i)*pz3(i)+py4(i)*pz4(i)
C
           aa = -(pxx+pyy+pzz)
           bb =  gfac*(pxx*pyy+pxx*pzz+pyy*pzz-pxy**2-pxz**2-pyz**2) 
           p  = bb-third*aa*aa
           d  = two*sqrt(third*max(-p,zero))-third*aa
C
           deltax(i) = one / sqrt(d)
         END DO
 
       END IF
C
      RETURN
C
 1000 FORMAT(/' ZERO OR NEGATIVE VOLUME : 3D-ELEMENT NB',i10/)
 2000 FORMAT(/' ZERO OR NEGATIVE VOLUME : DELETE 3D-ELEMENT NB',i10/)