Sources Esope | Cast3M

j3load
C J3LOAD    SOURCE    CHAT      05/01/13    00:46:59     5004
      SUBROUTINE J3LOAD(NUM,WWORK,NPTO,TOL,IRET)
C----------------------------------------------------
C     TRANSFERT DES FACE 3D EN FACE 2D
C
C     PP 9/97
C     Pierre Pegon/JRC Ispra
C----------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION NUM(2,NPTO),P1(3)
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
C
      SEGMENT WWORK
        REAL*8 PORIG(3),VNORM(3),VI(3),VJ(3)
        INTEGER FWORK
        INTEGER TWORK(NTROU)
      ENDSEGMENT
C
      SEGMENT WORK
        REAL*8 XYC(2,NPTO)
        INTEGER IST(3,NPTO)
        REAL*8 DENS(NPTO)
        INTEGER JUN
      ENDSEGMENT
C
      IRET=0
      IF(NPTO.LE.2)THEN
        WRITE(IOIMP,*)'J3LOAD: nb de point d"une face insuffisant'
        IRET=IRET+1
        RETURN
      ENDIF
      WORK=FWORK
C
C     LE PREMIER POINT DU CONTOURS DEFINIT L'ORIGINE
C
      NUME=NUM(1,1)
      IREF=(NUME-1)*(IDIM+1)
      DO IE1=1,3
        PORIG(IE1)=XCOOR(IREF+IE1)
      ENDDO
      XYC(1,1)=0.D0
      XYC(2,1)=0.D0
      DENS(1)=XCOOR(IREF+4)
C
C     AVEC LE SECOND POINT, ON DEFINIT LE PREMIER VECTEUR
C
      NUME=NUM(1,2)
      IREF=(NUME-1)*(IDIM+1)
      DO IE1=1,3
        P1(IE1)=XCOOR(IREF+IE1)
        VI(IE1)=P1(IE1)-PORIG(IE1)
      ENDDO
      XNORM=SQRT(VI(1)**2+VI(2)**2+VI(3)**2)
      IF(XNORM.LT.TOL)THEN
        WRITE(IOIMP,*)'J3LOAD: 2 point d"une face sont trop proches'
        IRET=IRET+1
        RETURN
      ENDIF
      DO IE1=1,3
        VI(IE1)=VI(IE1)/XNORM
      ENDDO
      XYC(1,2)=XNORM
      XYC(2,2)=0.D0
      DENS(2)=XCOOR(IREF+4)
C
C     AVEC LE POINT SUIVANT, OU UN AUTRE POINT, ON DEFINIT LE VECTEUR
C     NORMAL PUIS LE SECOND VECTEUR
C     WARNING: LA + GD COMPOSANTE DE VNOR EST TJ POSITIVE!
C
      DO IE1=3,NPTO
        NUME=NUM(1,IE1)
        IREF=(NUME-1)*(IDIM+1)
        DO IE2=1,3
          VJ(IE2)=XCOOR(IREF+IE2)-P1(IE2)
          P1(IE2)=XCOOR(IREF+IE2)
        ENDDO
        XNORM=SQRT(VJ(1)**2+VJ(2)**2+VJ(3)**2)
        IF(XNORM.LT.TOL)THEN
          WRITE(IOIMP,*)'J3LOAD: 2 point d"une face sont trop proches'
          IRET=IRET+1
          RETURN
        ENDIF
        DO IE2=1,3
          VJ(IE2)=VJ(IE2)/XNORM
        ENDDO
        VNORM(1)=VI(2)*VJ(3)-VI(3)*VJ(2)
        VNORM(2)=VI(3)*VJ(1)-VI(1)*VJ(3)
        VNORM(3)=VI(1)*VJ(2)-VI(2)*VJ(1)
        XNORM=SQRT(VNORM(1)**2+VNORM(2)**2+VNORM(3)**2)
        IF(XNORM.GT.100*TOL)THEN
          XCMAX=VNORM(1)
          DO IE2=2,3
            IF(ABS(VNORM(IE2)).GT.ABS(XCMAX))THEN
              XCMAX=VNORM(IE2)
            ENDIF
          ENDDO
          XNORM=SIGN(1.D0,XCMAX)*XNORM
          DO IE2=1,3
            VNORM(IE2)=VNORM(IE2)/XNORM
          ENDDO
          VJ(1)=VNORM(2)*VI(3)-VNORM(3)*VI(2)
          VJ(2)=VNORM(3)*VI(1)-VNORM(1)*VI(3)
          VJ(3)=VNORM(1)*VI(2)-VNORM(2)*VI(1)
          GOTO 1
        ENDIF
      ENDDO
      WRITE(IOIMP,*)'J3LOAD: On ne reussit pas a definir le 2nd vecteur'
      IRET=IRET+1
      RETURN
C
C     ON BOUCLE SUR LES POINTS ET ON LES PROJETTE SUR LE PLAN
C     ON TESTE LA DISTANCE DU POINT AU PLAN
C
 1    CONTINUE
      DO IE1=3,NPTO
        NUME=NUM(1,IE1)
        IREF=(NUME-1)*(IDIM+1)
        DO IE2=1,3
          P1(IE2)=XCOOR(IREF+IE2)
        ENDDO
        XDIST=(P1(1)-PORIG(1))*VNORM(1)+(P1(2)-PORIG(2))*VNORM(2)
     >                                 +(P1(3)-PORIG(3))*VNORM(3)
        IF(ABS(XDIST).GT.TOL)THEN
          WRITE(IOIMP,*)'J3LOAD: 1 point d"une face est hors plan'
          IRET=IRET+1
          RETURN
        ENDIF
        XYC(1,IE1)=(P1(1)-PORIG(1))*VI(1)+(P1(2)-PORIG(2))*VI(2)
     >                                   +(P1(3)-PORIG(3))*VI(3)
        XYC(2,IE1)=(P1(1)-PORIG(1))*VJ(1)+(P1(2)-PORIG(2))*VJ(2)
     >                                   +(P1(3)-PORIG(3))*VJ(3)
        DENS(IE1)=XCOOR(IREF+4)
      ENDDO
C
      JUN=0
C
      RETURN
      END