s4deri3.F File Reference

#include "implicit_f.inc"
#include "mvsiz_p.inc"
#include "scr03_c.inc"
#include "vect01_c.inc"
#include "param_c.inc"
#include "com01_c.inc"
#include "scr17_c.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	s4deri3 (vol, veul, geo, igeo, rx, ry, rz, sx, sy, sz, tx, ty, tz, x1, x2, x3, x4, y1, y2, y3, y4, z1, z2, z3, z4, px1, px2, px3, px4, py1, py2, py3, py4, pz1, pz2, pz3, pz4, jac_i, deltax, det, ngl, ngeo, mxt, pm, voldp)

Function/Subroutine Documentation

s4deri3()

subroutine s4deri3	(		vol,
			veul,
			geo,
		integer, dimension(npropgi,*)	igeo,
			rx,
			ry,
			rz,
			sx,
			sy,
			sz,
			tx,
			ty,
			tz,
		double precision, dimension(mvsiz)	x1,
		double precision, dimension(mvsiz)	x2,
		double precision, dimension(mvsiz)	x3,
		double precision, dimension(mvsiz)	x4,
		double precision, dimension(mvsiz)	y1,
		double precision, dimension(mvsiz)	y2,
		double precision, dimension(mvsiz)	y3,
		double precision, dimension(mvsiz)	y4,
		double precision, dimension(mvsiz)	z1,
		double precision, dimension(mvsiz)	z2,
		double precision, dimension(mvsiz)	z3,
		double precision, dimension(mvsiz)	z4,
			px1,
			px2,
			px3,
			px4,
			py1,
			py2,
			py3,
			py4,
			pz1,
			pz2,
			pz3,
			pz4,
			jac_i,
			deltax,
			det,
		integer, dimension(*)	ngl,
		integer, dimension(*)	ngeo,
		integer, dimension(mvsiz)	mxt,
			pm,
		double precision, dimension(*)	voldp )

Definition at line 37 of file s4deri3.F.

      USE message_mod
C-----------------------------------------------
C   I m p l i c i t   T y p e s
C-----------------------------------------------
#include      "implicit_f.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      "scr03_c.inc"
#include      "vect01_c.inc"
#include      "param_c.inc"
#include      "com01_c.inc"
#include      "scr17_c.inc"
C-----------------------------------------------
C   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER IGEO(NPROPGI,*), MXT(MVSIZ)
      double precision
     .   x1(mvsiz), x2(mvsiz), x3(mvsiz), x4(mvsiz), 
     .   y1(mvsiz), y2(mvsiz), y3(mvsiz), y4(mvsiz),
     .   z1(mvsiz), z2(mvsiz), z3(mvsiz), z4(mvsiz),voldp(*) 
C     REAL
      my_real
     .   vol(*), veul(lveul,*),geo(npropg,*), pm(npropm,*),
     .   rx(*), ry(*), rz(*), sx(*), sy(*), sz(*), tx(*), ty(*), tz(*), det(*),
     .   px1(mvsiz), px2(mvsiz), px3(mvsiz), px4(mvsiz),  
     .   py1(mvsiz), py2(mvsiz), py3(mvsiz), py4(mvsiz),  
     .   pz1(mvsiz), pz2(mvsiz), pz3(mvsiz), pz4(mvsiz),  
     .   deltax(mvsiz),jac_i(10,mvsiz)
      INTEGER NGL(*), NGEO(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .   b1(mvsiz), b2(mvsiz), b3(mvsiz), b4(mvsiz), 
     .   c1(mvsiz), c2(mvsiz), c3(mvsiz), c4(mvsiz),
     .   d1(mvsiz), d2(mvsiz), d3(mvsiz), d4(mvsiz) 
      double precision
     .   x41, y41, z41, x42, y42, z42, x43, y43, z43,b1dp,c1dp,d1dp
      my_real
     .   d ,pxx, pyy, pzz, pxy, pyz, pxz, gfac, aa, bb, p, ld
C=======================================================================
      DO i=lft,llt
       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)
       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
       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)
       vol(i) = det(i)
      ENDDO
C
      DO i=lft,llt
        IF (det(i) <= zero) THEN
          IF (igeo(11,ngeo(i)) /= 0) THEN
            CALL ancmsg(msgid=245,
     .                  msgtype=msgerror,
     .                  anmode=aninfo,
     .                  i1=ngl(i))
            vol(i) = em20    ! to prevent crash
            det(i) = em20    
          ELSE
            CALL ancmsg(msgid=635,
     .                  msgtype=msgwarning,
     .                  anmode=aninfo,
     .                  i1=ngl(i))
          ENDIF
        ENDIF
      ENDDO
C
      DO i=lft,llt
        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
      ENDDO
 
 
       IF(idt1sol==0)THEN
 
         DO i=lft,llt
             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=lft,llt
           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=lft,llt
           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=lft,llt
           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
 
      IF (ismstr==10.OR.ismstr==12) THEN
       DO i=lft,llt
        jac_i(1,i) = px1(i)
        jac_i(2,i) = px2(i)
        jac_i(3,i) = px3(i)
        jac_i(4,i) = py1(i)
        jac_i(5,i) = py2(i)
        jac_i(6,i) = py3(i)
        jac_i(7,i) = pz1(i)
        jac_i(8,i) = pz2(i)
        jac_i(9,i) = pz3(i)
        jac_i(10,i) = vol(i)
       ENDDO
      END IF
C
      IF (jeul * (1 - imulti_fvm) /= 0 .AND. nxref==0) THEN
        DO i=lft,llt
          veul(1,i) =-px1(i)
          veul(2,i) = px2(i)
          veul(3,i) =-px3(i)
          veul(4,i) = px4(i)
          veul(5,i) =-py1(i)
          veul(6,i) = py2(i)
          veul(7,i) =-py3(i)
          veul(8,i) = py4(i)
          veul(9,i) =-pz1(i)
          veul(10,i)= pz2(i)
          veul(11,i)=-pz3(i)
          veul(12,i)= pz4(i)
          veul(13,i)= deltax(i)
C------------------------------------------
C     CALCUL DE LA NORMALE A CHAQUE FACE
C------------------------------------------
          veul(14,i) = zero
          veul(15,i) = half*b1(i)
          veul(16,i) = zero
          veul(17,i) = half*b2(i)
          veul(18,i) = half*b3(i)
          veul(19,i) = half*b4(i)
C
          veul(20,i) = zero
          veul(21,i) = half*c1(i)
          veul(22,i) = zero
          veul(23,i) = half*c2(i)
          veul(24,i) = half*c3(i)
          veul(25,i) = half*c4(i)
C
          veul(26,i) = zero
          veul(27,i) = half*d1(i)
          veul(28,i) = zero
          veul(29,i) = half*d2(i)
          veul(30,i) = half*d3(i)
          veul(31,i) = half*d4(i)
        ENDDO
C
        IF (igeo(11,ngeo(lft)) == 15) THEN
          DO i=lft,llt
            vol(i)=vol(i)*geo(1,ngeo(i))    
            voldp(i)=voldp(i)*geo(1,ngeo(i))    
          ENDDO
        ENDIF
      ENDIF
C-----------
      RETURN