C TUFIST    SOURCE    CHAT      05/01/13    03:56:20     5004
      SUBROUTINE TUFIST(VALMAT,XDDL,CAR,DD,BPSS,XDDLOC,XSTRS,KERRE)
C=================================================================
C     CALCULE LES CONTRAINTES DU TUYAU FISSURE
C  ENTREES
C     VALMAT   = TABLEAU DE MATERIAU
C     XDDL(12) = D.D.L. DU TUYAU FISSURE  GLOBAUX
C     CAR(9)   = CARACTERISTIQUES
C                      CAR(1) = RAYON  EXTERIEUR DU TUYAU
C                      CAR(2) = EPAISSEUR DU TUYAU
C                      CAR(3),CAR(4),CAR(5)  AXE DU TUYAU
C                      CAR(6),CAR(7),CAR(8)  DIRECTION DE LA FISSU
C                      CAR(9)  = ANGLE DE LA FISSURE
C     DD(6,6)  = MATRICE DE HOOKE
C  TRAVAIL
C     BPSS(3,3)=     ...
C     XDDLOC(12) = D.D.L. LOCAUX
C  SORTIES
C     XSTRS(7) = CONTRAINTES
C                    XSTRS(1),XSTRS(2),XSTRS(3)  FORCES
C                    XSTRS(4),XSTRS(5),XSTRS(6)  MOMENTS
C                    XSTRS(7)  FACTEUR D INTENSITE DES CONTRAINTES
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 TUYAU TROP EPAI
C     CODE EBERSOLT SEPTEMBRE 85
C
C=================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
      PARAMETER(UN=1.D0,UNS2=.5D0,DEUX=2.D0)
      PARAMETER(EPS=1.D-10)
C
      DIMENSION XDDL(*),V1(3),V2(3),BPSS(3,*),XDDLOC(*)
      DIMENSION XSTRS(*),DD(6,*),BV1(3),BV2(3),CAR(*)
      DIMENSION VALMAT(*)
C
C   D6 = COEFFICIENT INTERVENANT POUR LES TERMES DE PENALISATION.
C
      KERRE=0
      YOU=VALMAT(1)
      RAYO=CAR(1)
      EPAI=CAR(2)
      V1(1)=CAR(3)
      V1(2)=CAR(4)
      V1(3)=CAR(5)
      V2(1)=CAR(6)
      V2(2)=CAR(7)
      V2(3)=CAR(8)
      ANGL=CAR(9)
*
      RAYMO= RAYO - (EPAI/DEUX)
      TETA1= ANGL*UNS2
C     CONVERSION DE TETA1 EN RADIAN
      TETA= (TETA1*XPI)/180.D0
      TESPI = TETA/XPI
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
      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     V1 REPRESENTE L'AXE LOCAL DES Z
C     V2 REPRESENTE L'AXE LOCAL DES Y
C
      DO 20 I=1,3
         BPSS(1,I) =BV1(I)
         BPSS(2,I) =BV2(I)
 20   CONTINUE
      BPSS(3,1) = (BPSS(1,2)*BPSS(2,3))-(BPSS(1,3)*BPSS(2,2))
      BPSS(3,2) = (BPSS(1,3)*BPSS(2,1))-(BPSS(1,1)*BPSS(2,3))
      BPSS(3,3) = (BPSS(1,1)*BPSS(2,2))-(BPSS(1,2)*BPSS(2,1))
C
C     TRANSFORMATION DES DDL GLOBAUX EN DDL LOCAUX X AXE TUYAU
C                                                  Y POSITION FISSURE
C                                                  Z COMPLEMENT
      CALL MATVEC(XDDL,XDDLOC,BPSS,4)
C
C     TRANCHANT AXE Y
C
      XSTRS(2)= DD(2,2)*(XDDLOC(8 )-XDDLOC(2))
C
C     TRANCHANT AVEC L AXE   Z
C
      XSTRS(3)= DD(3,3)*(XDDLOC(9 )-XDDLOC(3))
C
C     TORSION AXE X
C
      XSTRS(4)= DD(4,4)*(XDDLOC(10)-XDDLOC(4))
C
C     FLEXION AXE Y
C
      XSTRS(5)= DD(5,5)*(XDDLOC(11)-XDDLOC(5))
C
      XX1= XDDLOC(7 )-XDDLOC(1)
      XX2= XDDLOC(12)-XDDLOC(6)
C
C     OUVERTURE TUFI ET ROTATION TUFI
C
      XSTRS(1)= XX1*DD(1,1)+XX2*DD(1,6)
      XSTRS(6)= XX1*DD(6,1)+XX2*DD(6,6)
C
C     CALCUL DE A   COEFFICIENT ZAHOR
C
      RSURT=RAYMO / EPAI
      IF(RSURT.LE.10.D0.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
      CALL TUFIFP(TESPI,AXX,FP,FM,FMP,FOP,FOM)
C
C
C     FACTEUR D INTENSITE DES CONTRAINTES
C
      IF(TETA1.LE.(0.5D0))THEN
         XSTRS(7)=XZERO
         GOTO 666
      ENDIF
      SQQ= XPI * RAYMO * TETA
      SQQ= SQRT(SQQ)
      XEX= SQQ * FOP/(DEUX * XPI * RAYMO *EPAI)
      XFL= SQQ * FOM/(XPI * RAYMO * RAYMO *EPAI)
      XSTRS(7)=XEX * XSTRS(1) - XFL * XSTRS(6)
C
C     CALCUL DES AIRES DE BRECHE NOTE TECHNIQUE DRE/STRE/LMA 85/695
C
      SIGM = XSTRS(1)/( DEUX * XPI * RAYMO * EPAI )
      SIGF = XSTRS(6)/( XPI * RAYMO * RAYMO * EPAI )
      XIM =  XPI * RAYMO * RAYMO * DEUX * TETA * TETA * FP /YOU
      XIF =  XIM * ( .75D0 +(.25D0 * COS ( TETA )))
      XSTRS(8)=XIM * SIGM  - XIF * SIGF
C
 666  CONTINUE
      RETURN
      END