Sources Esope | Cast3M

conge4
C CONGE4    SOURCE    CHAT      05/01/12    22:17:15     5004
      SUBROUTINE CONGE4 (L1,L2,RCONGE,NCONGE, NL1,LRAC,NL2)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
************************************************************************
*
*                             C O N G E 4
*                             -----------
*
* FONCTION:
* ---------
*
*     CREER UN CONGE DE RACCORDEMENT ENTRE DEUX LIGNES NON COPLANAIRES
*
*     TYPE DU CONGE : DOUBLE COUDE
*
*
* MODULES UTILISES:
* -----------------
*
-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC CCREEL
*
* PARAMETRES:   (E)=ENTREE   (S)=SORTIE   (+ = CONTENU DANS UN COMMUN)
* -----------
*
*     L1      (E)  PREMIERE LIGNE A RACCORDER
*     L2      (E)  DEUXIEME LIGNE A RACCORDER
*     RCONGE  (E)  RAYON DU CONGE DU RACCORDEMENT
*     NCONGE  (E)  FACTEUR DE DECOUPAGE DU CONGE (CONVENTIONS
*                  CLASSIQUES SUR SON SIGNE)
*                  = 0  SI NON FOURNI
*     NL1     (S)  LIGNE APPUYEE SUR L1 ET ABOUTISSANT SUR LE CONGE
*     LRAC    (S)  CONGE DE RACCORDEMENT
*     NL2     (S)  LIGNE APPUYEE SUR L2 ET COMMENCANT SUR LE CONGE
*
      INTEGER NCONGE
      REAL*8 RCONGE
*
* CONSTANTES:
* -----------
*
*
      PARAMETER (DEMI=0.5D0)
*
* VARIABLES:
* ----------
*
      INTEGER N,NBNOEU,NBELEM,NL,ND,
     &        ITYPLM,NBSOUS,NBREF,NBNN
      REAL*8 MU,L
*
*
* AUTEUR, DATE DE CREATION:
* -------------------------
*
*     GILLES DUVERGER       06 OCTOBRE 1988
*
* LANGAGE:
* --------
*
*     ESOPE + FORTRAN77
*
************************************************************************
*
      MELEME=L1
      SEGACT,MELEME
      NL=LISOUS(/1)
      IF (NL.NE.0) THEN
         CALL ERREUR(25)
         RETURN
      END IF
      NBNOEU=NUM(/1)
      NBELEM=NUM(/2)
      IP9=NUM(NBNOEU,NBELEM)
      IP8=NUM(1,NBELEM)
      SEGDES,MELEME
*
      MELEME=L2
      SEGACT,MELEME
      NL=LISOUS(/1)
      IF (NL.NE.0) THEN
         CALL ERREUR(25)
         RETURN
      END IF
      NBNOEU=NUM(/1)
      IP1=NUM(1,1)
      IP2=NUM(NBNOEU,1)
      SEGDES,MELEME
*
      CALL EXCOO1(IP9,X9,Y9,Z9,D9)
      CALL EXCOO1(IP8,X8,Y8,Z8,D8)
      CALL EXCOO1(IP1,X1,Y1,Z1,D1)
      CALL EXCOO1(IP2,X2,Y2,Z2,D2)
*
*     VECTEURS DIRECTEURS DE L1 ET L2
*
      XN9 = SQRT((X8-X9)**2+(Y8-Y9)**2+(Z8-Z9)**2)
      XN1 = SQRT((X1-X2)**2+(Y1-Y2)**2+(Z1-Z2)**2)
      A1=(X8-X9)/XN9
      B1=(Y8-Y9)/XN9
      C1=(Z8-Z9)/XN9
      A2=(X2-X1)/XN1
      B2=(Y2-Y1)/XN1
      C2=(Z2-Z1)/XN1
*
*     RECHERCHE DES POINTS M1 SUR L1 ET M2 SUR L2 TELS QUE
*     D(M1;M2) = D(L1;L2)
*
      B=-(A1*A2+B1*B2+C1*C2)
*
      D=(X1-X9)*A1+(Y1-Y9)*B1+(Z1-Z9)*C1
      E=-((X1-X9)*A2+(Y1-Y9)*B2+(Z1-Z9)*C2)
      ALPHA=(B*E-D)/(B**2-1)
      BETA=(B*D-E)/(B**2-1)
*
      XM1=X9+ALPHA*A1
      YM1=Y9+ALPHA*B1
      ZM1=Z9+ALPHA*C1
      CALL CREPO1(XM1,YM1,ZM1,IPM1)
      XM2=X1+BETA*A2
      YM2=Y1+BETA*B2
      ZM2=Z1+BETA*C2
      CALL CREPO1(XM2,YM2,ZM2,IPM2)
*
*     DISTANCE ENTRE LES DEUX DROITES
*
      D2D=SQRT((XM2-XM1)**2+(YM2-YM1)**2+(ZM2-ZM1)**2)
*
*     TEST DE LA VALEUR DE RCONGE
*
      IF (RCONGE.LT.(D2D*DEMI)) THEN
          CALL ERREUR(401)
          RETURN
      END IF
*
*     CALCUL PREALABLE
*
      BB=D2D*D2D
      AA=2*(1+B)
      MU=-(1+B)
*
      S=(BB*AA+2*AA*(BB-4*RCONGE*RCONGE)-2*MU*MU*(BB+4*RCONGE*RCONGE))
     &/(AA*(AA-MU*MU))
      T=(AA*(BB-4*RCONGE*RCONGE)*(BB-4*RCONGE*RCONGE)+2*AA*BB*(BB-4*
     &RCONGE*RCONGE)-MU*MU*(BB+4*RCONGE*RCONGE)*(BB+4*RCONGE*RCONGE))/
     &(AA*AA*(AA-MU*MU))
      U=BB*(BB-4*RCONGE*RCONGE)/(AA*AA*(AA-MU*MU))
*
      CALL DEGRE3(U,T,S,XR1,XI1,XR2,XI2,XR3,XI3)
*
      QQ=T/3-S*S/9
      RR=(T*S-3*U)/6-S*S*S/27
      D=QQ*QQ*QQ+RR*RR
      IF (D .LT. 0) THEN
         X=MAX(XR1,XR2,XR3)
      ELSE
         X=XR1
      END IF
      L=SQRT(AA*X+BB)*COS(XPI-2*ATAN(SQRT(AA*X+BB)/(2*RCONGE)))/MU
*
*     RECHERCHE DU POINT P1 SUR L1
*
      XP1=XM1+L*A1
      YP1=YM1+L*B1
      ZP1=ZM1+L*C1
      CALL CREPO1(XP1,YP1,ZP1,IPP1)
*
*     RECHERCHE DU POINT P2 SUR L2
*
*   JE PENSE QUE LAMBDA VAUT 1     PV
      LAMBDA=1
      XP2=XM2+LAMBDA*L*A2
      YP2=YM2+LAMBDA*L*B2
      ZP2=ZM2+LAMBDA*L*C2
      CALL CREPO1(XP2,YP2,ZP2,IPP2)
*
*     RECHERCHE DU POINT O ,MILIEU DE P1 P2
*
      XO=(XP1+XP2)*DEMI
      YO=(YP1+YP2)*DEMI
      ZO=(ZP1+ZP2)*DEMI
      CALL CREPO1(XO,YO,ZO,IPO)
*
*     CREATION DE LA DROITE P1-O
*
      ND=1
      CALL ECRENT(ND)
      CALL ECROBJ('POINT',IPO)
      CALL ECROBJ('POINT',IPP1)
      CALL LIGNE(1,1,DEN1,DEN2,N)
      IF (IERR.NE.0) RETURN
      CALL LIROBJ('MAILLAGE',L3,1,IRETOU)
      IF (IERR.NE.0) RETURN
*
*     CREATION DE LA DROITE O-P2
*
      ND=1
      CALL ECRENT(ND)
      CALL ECROBJ('POINT',IPP2)
      CALL ECROBJ('POINT',IPO)
      CALL LIGNE(1,1,DEN1,DEN2,N)
      IF (IERR.NE.0) RETURN
      CALL LIROBJ('MAILLAGE',L4,1,IRETOU)
      IF (IERR.NE.0) RETURN
*
*     CREATION DU CONGE DOUBLE
*
      NCONGE=(NCONGE+1)/2
      CALL CONGE2(L1,L4,RCONGE,NCONGE,NL1,LRAC1,NL4)
      IF (IERR.NE.0) RETURN
      CALL CONGE2(L3,L2,RCONGE,NCONGE,NL3,LRAC2,NL2)
      IF (IERR.NE.0) RETURN
*
*     LES 2 ARCS DE CERCLE ONT UN POINT EN COMMUN:
      MELEME=LRAC1
      SEGACT,MELEME
      NBNOEU=NUM(/1)
      NBELEM=NUM(/2)
      IPO1=NUM(NBNOEU,NBELEM)
      SEGDES,MELEME
      MELEME=LRAC2
      SEGACT,MELEME
      NUM(1,1)=IPO1
      SEGDES,MELEME
*     ON REUNIT LES 2 DEMI-CERCLES:
      CALL ECROBJ('MAILLAGE',LRAC2)
      CALL ECROBJ('MAILLAGE',LRAC1)
      CALL PRFUSE
      IF (IERR.NE.0) RETURN
      CALL LIROBJ('MAILLAGE',LRAC,1,IRETOU)
      IF (IERR.NE.0) RETURN
      MELEME=LRAC1
      SEGSUP,MELEME
      MELEME=LRAC2
      SEGSUP,MELEME
*
      MELEME=L3
      SEGSUP,MELEME
      MELEME=NL3
      SEGSUP,MELEME
      MELEME=L4
      SEGSUP,MELEME
      MELEME=NL4
      SEGSUP,MELEME
*
 
      END