i11dst3.F File Reference

#include "implicit_f.inc"
#include "mvsiz_p.inc"

Go to the source code of this file.

Functions/Subroutines
subroutine	i11dst3 (jlt, gap, cand_s, cand_m, irects, irectm, nx, ny, nz, n1, n2, m1, m2, jlt_new, x, igap, gap_s, gap_m, gapv2, gap_s_l, gap_m_l, drad, dgapload)

Function/Subroutine Documentation

i11dst3()

subroutine i11dst3	(	integer	jlt,
			gap,
		integer, dimension(*)	cand_s,
		integer, dimension(*)	cand_m,
		integer, dimension(2,*)	irects,
		integer, dimension(2,*)	irectm,
			nx,
			ny,
			nz,
		integer, dimension(*)	n1,
		integer, dimension(*)	n2,
		integer, dimension(*)	m1,
		integer, dimension(*)	m2,
		integer	jlt_new,
			x,
		integer	igap,
			gap_s,
			gap_m,
			gapv2,
			gap_s_l,
			gap_m_l,
			drad,
		intent(in)	dgapload )

Definition at line 28 of file i11dst3.F.

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   D u m m y   A r g u m e n t s
C-----------------------------------------------
      INTEGER JLT,JLT_NEW,IGAP
      INTEGER IRECTS(2,*), IRECTM(2,*),CAND_S(*),CAND_M(*),
     . N1(*),N2(*),M1(*),M2(*)
      my_real , INTENT(IN) :: dgapload
      my_real
     .     gap,drad,
     . nx(*),ny(*),nz(*),x(3,*),gap_s(*),gap_m(*), gapv2(*),
     . gap_s_l(*),gap_m_l(*)
C-----------------------------------------------
C   L o c a l   V a r i a b l e s
C-----------------------------------------------
      INTEGER I
      my_real
     .     xs12,ys12,zs12,xm12,ym12,zm12,xa,xb,
     .     xs2,xm2,xsm,xs2m2,ys2,ym2,ysm,ys2m2,zs2,zm2,zsm,zs2m2,
     .     xx,yy,zz,als,alm,det,h1s,h2s,h1m,h2m,
     .     gap2,drad2,
     .     pene2(mvsiz)
C-----------------------------------------------
       gap2=(gap+dgapload)*(gap+dgapload)
       drad2=drad*drad
C
      IF(igap==0)THEN
        DO i=1,jlt
          gapv2(i)=gap2
        ENDDO
      ELSE
        DO i=1,jlt
          gapv2(i)=gap_s(cand_s(i))+gap_m(cand_m(i))
         IF(igap == 3)  
     .     gapv2(i)=min(gap_s_l(cand_s(i))+gap_m_l(cand_m(i)),gapv2(i))+dgapload
          gapv2(i)=max(gapv2(i)*gapv2(i),gap2)
        ENDDO
      ENDIF
C--------------------------------------------------------
C  
C--------------------------------------------------------
C       F = [A*X1+(1-A)*X2-B*X3-(1-B)*X4]^2 + [..Y..]^2 + [..Z..]^2
C       DF/DA = 0 = (X1-X2)(A(X1-X2)+X2-X4 +B(X4-X3))+...
C       DF/DA = 0 = A(X1-X2)^2 +X2-X4 + B(X1-X2)(X4-X3))+...
C       DF/DA = 0 = A[(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2] 
C                 + B[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
C                 +   (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4) 
C       DF/DB = 0 = (X4-X3)(A(X1-X2)+X2-X4 +B(X4-X3))+...
C       DF/DB = 0 = B[(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2] 
C                 + A[(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
C                 +   (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4) 
C       XS2 = [(X1-X2)^2 + (Y1-Y2)^2 + (Z1-Z2)^2]
C       XM2 = [(X4-X3)^2 + (Y4-Y3)^2 + (Z4-Z3)^2]
C       XSM = [(X1-X2)(X4-X3) + (Y1-Y2)(Y4-Y3) + (Z1-Z2)(Z4-Z3)]
C       XA = (X1-X2)(X2-X4) + (Y1-Y2)(Y2-Y4) + (Z1-Z2)(Z2-Z4)
C       XB = (X4-X3)(X2-X4) + (Y4-Y3)(Y2-Y4) + (Z4-Z3)(Z2-Z4)
C       A XS2 + B XSM +   XA = 0
C       A XSM + B XM2 +   XB = 0
C
C       A = -(XA + B XSM)/XS2
C       -(XA + B XSM)*XSM + B XM2*XS2 +   XB*XS2 = 0
C       -B XSM*XSM + B XM2*XS2 +   XB*XS2-XA*XSM  = 0
C       B*(XM2*XS2 - XSM*XSM) = -XB*XS2+XA*XSM  
C       B = (XA*XSM-XB*XS2) / (XM2*XS2 - XSM*XSM)
C       A = (XB*XSM-XA*XM2) / (XM2*XS2 - XSM*XSM)
      DO i=1,jlt
       n1(i)=irects(1,cand_s(i))
       n2(i)=irects(2,cand_s(i))
       m1(i)=irectm(1,cand_m(i))
       m2(i)=irectm(2,cand_m(i))
       xs12 = x(1,n2(i))-x(1,n1(i))
       ys12 = x(2,n2(i))-x(2,n1(i))
       zs12 = x(3,n2(i))-x(3,n1(i))
       xs2 = xs12*xs12 + ys12*ys12 + zs12*zs12
       xm12 = x(1,m2(i))-x(1,m1(i))
       ym12 = x(2,m2(i))-x(2,m1(i))
       zm12 = x(3,m2(i))-x(3,m1(i))
       xm2 = xm12*xm12 + ym12*ym12 + zm12*zm12
       xsm = - (xs12*xm12 + ys12*ym12 + zs12*zm12)
       xs2m2 = x(1,m2(i))-x(1,n2(i))
       ys2m2 = x(2,m2(i))-x(2,n2(i))
       zs2m2 = x(3,m2(i))-x(3,n2(i))
       xa =  xs12*xs2m2 + ys12*ys2m2 + zs12*zs2m2
       xb = -xm12*xs2m2 - ym12*ys2m2 - zm12*zs2m2 
       det = xm2*xs2 - xsm*xsm
       det = max(em20,det)
C
       h1m = (xa*xsm-xb*xs2) / det
C
       xs2 = max(xs2,em20)
       xm2 = max(xm2,em20)
       h1m=min(one,max(zero,h1m))
       h1s = -(xa + h1m*xsm) / xs2
       h1s=min(one,max(zero,h1s))
       h1m = -(xb + h1s*xsm) / xm2
       h1m=min(one,max(zero,h1m))
C
       h2s = one - h1s
       h2m = one - h1m
C !!!!!!!!!!!!!!!!!!!!!!!
C  PENE = GAP^2 - DIST^2 UTILISE POUR TESTER SI NON NUL
C!!!!!!!!!!!!!!!!!!!!!!!!
       nx(i) =  h1s*x(1,n1(i)) + h2s*x(1,n2(i))
     .      - h1m*x(1,m1(i)) - h2m*x(1,m2(i))
       ny(i) =  h1s*x(2,n1(i)) + h2s*x(2,n2(i))
     .      - h1m*x(2,m1(i)) - h2m*x(2,m2(i))
       nz(i) =  h1s*x(3,n1(i)) + h2s*x(3,n2(i))
     .      - h1m*x(3,m1(i)) - h2m*x(3,m2(i))
       pene2(i) = max(gapv2(i),drad2) - nx(i)*nx(i) - ny(i)*ny(i) - nz(i)*nz(i)
       pene2(i) = max(zero,pene2(i)) 
C
      ENDDO
      DO i=1,jlt
        IF(pene2(i)/=zero)THEN
          jlt_new = jlt_new + 1
          cand_s(jlt_new) = cand_s(i)
          cand_m(jlt_new) = cand_m(i)
          n1(jlt_new)  = n1(i)
          n2(jlt_new)  = n2(i)
          m1(jlt_new)  = m1(i)
          m2(jlt_new)  = m2(i)
          nx(jlt_new)  = nx(i)
          ny(jlt_new)  = ny(i)
          nz(jlt_new)  = nz(i)
          gapv2(jlt_new)  = gapv2(i)
        ENDIF
      ENDDO
C
      RETURN