C TUFIR1    SOURCE    CHAT      05/01/13    03:56:11     5004
      SUBROUTINE TUFIR1(REG,RAYON,EPAI,V1,V2,ANDEG,D,REL,BPSS,KERRE)
C=======================================================================
C   CALCUL DE LA MATRICE DE RIGIDITE DE L'ELEMENT DE TUYAU FISSURE
C   DANS LE REPERE GLOBAL
C   ENTREES:
C     WORK  = TABLEAU CONTENANT LES CARACTERISTIQUES GEOMETRIQUES DE
C             L'ELEMENT
C     D(2,2)= MATRICE DE HOOKE
C  TRAVAIL
C     REL(12,12)   = MATRICE DE RIGIDITE  LOCALE
C     V1(3),V2(3)  = VECTEURS  ORIENTANT LE TUYAU FISSURE
C     BPSS(3,3)    = MATRICE DE PASSAGE
C   SORTIE:
C     REG(12,12)=  MATRICE DE RIGIDITE  GLOBALE
C     KERRE = CODE D'ERREUR CORRESPONDANT A LA POSITION RELATIVE DES
C             DEUX VECTEURS V1 ET V2 DEFINISSANT L'ELEMENT TUFI:
C             0 SI PAS DE PROBLEME, 1 SINON.
C             1 SI V2 EST NUL OU SI V1 ET V2 SONT COLINEAIRES.
C             2 SI V1 EST NUL.
C             3 SI LE TUYAU FISSURE EST TRP EPAI R/T  INF A 2
C=======================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
      DIMENSION REL(12,*),V1(*),V2(*),D(2,*),BPSS(3,*)
      DIMENSION REG(12,*),BV1(3),BV2(3)
C
      DATA EPS/1.D-10/
      DATA UN,UNS2/1.D0,.5D0/
C
C   D6 = COEFFICIENT INTERVENANT POUR LES TERMES DE PENALISATION.
C
      DATA DEUX,D6/2.D0,1.D06/
C
      KERRE = 0
      TETA1 = ANDEG*UNS2
      YOU   = D(1,1)
      XNU   = D(2,2)
      CISA  = YOU /(DEUX *(UN + XNU))
C=======================================================================
C   LES CARACTERISTIQUES REPRESENTENT:
C      EPAI   = EPAISSEUR DU TUYAU
C      RAYON  = RAYON EXTERIEUR DU TUYAU
C      TETA1  = DEMI-ANGLE D'OUVERTURE DE LA FISSURE EN DEGRE
C      V1     = VECTEUR ORIENTANT LES NOEUDS 1 ET 2 (AXE DU TUYAU)
C      V2     = VECTEUR ORIENTANT L'AXE DE LA FISSURE
C      YOU    = MODULE D'YOUNG
C      CISA   = MODULE DE CISAILLEMENT
C   TABLEAU DE TRAVAIL:
C      BPSS   = MATRICE DE PASSAGE
C   ON OBTIENT:
C      REG(12,12) = MATRICE DE RIGIDITE  AXES GLOBAUX
C
C=======================================================================
C
C     CONVERSION DE TETA1 EN RADIAN
C
      TETA = (TETA1 * XPI)/180.D0
C
C     ON MET DANS 'RAYON' LE RAYON MOYEN DU TUYAU.
C
      RAYMO = RAYON - EPAI*UNS2
C
C     ON CALCULE A
C
      RSURT= RAYMO / EPAI
      IF(RSURT.LE.10.0D0.AND.RSURT.GE.4.9D0) THEN
         AXX= (.125D0*RSURT - .25D0 ) **.25D0
      ELSE IF(RSURT.GT.10.D0.AND.RSURT.LE.35.D0) THEN
         AXX= (.4D0*RSURT - 3.D0 ) **.25D0
      ELSE
         KERRE=3
      ENDIF
      IF(KERRE.EQ.3) GOTO 666
C
C     NORMALISATION DES VECTEURS V1 ET V2.
C
      XNV1 = XZERO
      XNV2 = XZERO
      DO 5 I=1,3
         XNV1 = XNV1 + (V1(I)*V1(I))
         XNV2 = XNV2 + (V2(I)*V2(I))
 5    CONTINUE
C
C     MISE A ZERO DE REG
C
      CALL ZERO(REG,12,12)
      IF (XNV1.LT.EPS) THEN
         KERRE = 2
         GOTO 666
      ENDIF
      IF (XNV2.LT.EPS) THEN
         KERRE = 1
         GOTO 666
      ENDIF
      XNV1 = UN/SQRT(XNV1)
      XNV2 = UN/SQRT(XNV2)
      DO 10 I=1,3
         BV1(I) = V1(I)*XNV1
         BV2(I) = V2(I)*XNV2
 10   CONTINUE
C
C     ORTHOGONALISATION ET RENORMALISATION DE V2.
C
      PS = (BV1(1)*BV2(1))+(BV1(2)*BV2(2))+(BV1(3)*BV2(3))
C
C     TEST SUR LA COLINEARITE DE V1 ET V2.
C
      IF(ABS(PS).GE.(.99D0)) THEN
         KERRE = 1
         GOTO 666
      ELSE
         XNV2 = UN/SQRT(UN-(PS**2))
         DO 15 I=1,3
            BV2(I) = XNV2*(BV2(I)-(PS*BV1(I)))
 15      CONTINUE
      ENDIF
C
C     CALCUL DE LA MATRICE DE PASSAGE
C      LOCAL = BPSS * GLOBAL
C    BV1 REPRESENTE L'AXE LOCAL DES Z
C    BV2 REPRESENTE L'AXE LOCAL DES Y
C
      DO 20 I=1,3
         BPSS(3,I) =BV1(I)
         BPSS(2,I) =BV2(I)
 20   CONTINUE
      BPSS(1,1) = (BPSS(2,2)*BPSS(3,3))-(BPSS(2,3)*BPSS(3,2))
      BPSS(1,2) = (BPSS(2,3)*BPSS(3,1))-(BPSS(2,1)*BPSS(3,3))
      BPSS(1,3) = (BPSS(2,1)*BPSS(3,2))-(BPSS(2,2)*BPSS(3,1))
C
C     CALCUL DES TERMES DE LA MATRICE DE RIGIDITE     REPERE LOCAL
C
      TESPI = TETA/XPI
      CALL ZERO(REL,12,12)
C
C     ON TESTE SUR L'IMPORTANCE DE L'ANGLE D'OUVERTURE DE LA FISSURE. SI
C     CE DERNIER EST INFERIEUR A UN DEGRE, ON ANNULE TOUS LES COEFFI-
C     CIENTS DE COUPLAGE ET ON PENALISE TOUS LES DEGRES DE LIBERTE. DANS
C     LE CAS CONTRAIRE, SEULS SONT PERMIS LA ROTATION AUTOUR DE L'AXE
C     LOCAL DES X ET L'ALLONGEMENT SUIVANT L'AXE LOCAL DES Z ,CES DEUX
C     DEGRES DE LIBERTE ETANT COUPLES.
C
      IF (TETA1.LE.(0.5D0)) THEN
         REL(3,3) = DEUX * YOU * XPI * EPAI * RAYMO      * D6
         REL(4,4) =        YOU * XPI * EPAI * (RAYMO**2) * D6
      ELSE
C
         CALL TUFIFP(TESPI,AXX,FP,FM,FMP,FOP,FOM)
         DENOM = FP*FM - FMP*FMP
         COEF = (YOU * XPI * EPAI)/(DENOM * (TETA**2))
         REL(3,3) =  COEF * DEUX * FM
         REL(3,4) = -COEF * RAYMO * FMP
         REL(4,3) =  REL(3,4)
         REL(4,4) =  COEF * (RAYMO**2) * FP * UNS2
      ENDIF
C
C     INITIALISATION DES TERMES DE PENALISATION
C
      REL(1,1) = DEUX * CISA * XPI * EPAI * RAYMO      * D6
      REL(2,2) = REL(1,1)
      REL(5,5) =         YOU * XPI * EPAI * (RAYMO**2) * D6
      REL(6,6) = DEUX * CISA * XPI * EPAI * (RAYMO**2) * D6
C
C     CALCUL DE LA MATRICE (12,12) A PARTIR DE LA MATRICE (6,6)
C
      DO 30 I=1,6
         DO 301 J=1,6
            REL(I+6,J+6) =  REL(I,J)
            REL(I+6,J  ) = -REL(I,J)
            REL(I  ,J+6) = -REL(I,J)
 301     CONTINUE
 30   CONTINUE
      CALL POUROT(REG,12,BPSS,REL)
 666  CONTINUE
      RETURN
      END