Sources Esope | Cast3M

bcoq2
C BCOQ2     SOURCE    CHAT      06/03/29    21:15:33     5360
      SUBROUTINE BCOQ2(B,NSTRS,DJAC,IGAU,IFOU,XEL,NN,T,P,EXCEN,DIM3,
     .     IARR,XDPGE,YDPGE)
C==================================================================
C           CALCUL DE LA MATRICE B DES COQUES @ 2 NOEUDS
C                ROUTINE FORTRAN PUR
C           GELE   JANVIER 87
C==================================================================
C ENTREES
C     IGAU=NUMERO DU POINT DE GAUSS
C     IFOU=IFOUR DE CCOPTIO
C     XEL=COORDONNEES LOCALE DE L'ELEMENT
C     NN=NUMERO DU MODE DE FOURIER
C     T(IGAU)=POSITION DU POINT DE GAUSS
C     P(IGAU)=POIDS DU POINT DE GAUSS
C     EXCEN  = EXCENTREMENT
C     DIM3   = EPAISSEUR DANS L'AUTRE DIMENSION
C     XDPGE,YDPGE : COORDONNEE DU POINT AUTOUR DUQUEL
C          FAIT LE MOUVEMENT EN DEFO PLAN GENE
C SORTIE
C     B(6,8)=MATRICE B AU POINT DE GAUSSPOUR IFOU GT 0
C     B(4,6)=MATRICE B AU POINT DE GAUSSPOUR IFOU LE 0
C     B(4,9)=MATRICE B AU POINT DE GAUSSPOUR IFOU EQ -3
C     DJAC=JACOBIEN AU POINT DE GAUSS=POIGAU*LONG/2 (*R(IGAU), SI
C          IFOU EST SUPERIEUR OU EGAL A ZERO)
C     IARR=0 SI OK 1 SI LONGUEUR ELEMENT NULLE
C                  2 SI R / D INFERIEUR A 10-3
C==================================================================
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
C    Include contenant quelques constantes dont XPI, XZERO :
-INC CCREEL
      PARAMETER(UNDE=.5D0,UN=1.D0,DEUX=2.D0,TRS=3.D0)
      PARAMETER(QUTR=4.D0,SIX=6.D0,DOUZ=12.D0)
      DIMENSION B(NSTRS,*),T(*),XEL(3,*),P(*)
C
C     ---------------------------------INITIALISATION
      IARR=0
      DJAC=XZERO
      IF(IFOU.GT.0)  THEN
         LRE=8
         CALL ZERO(B,6,8)
      ELSE IF(IFOU.LE.0) THEN
         LRE=6
         CALL ZERO(B,4,6)
         IF(IFOU.EQ.-3) THEN
            LRE=9
            CALL ZERO(B,4,9)
         ENDIF
      ENDIF
C
      D=SQRT((XEL(1,2)-XEL(1,1))**2+(XEL(2,2)-XEL(2,1))**2)
      IF(D.EQ.0) THEN
         IARR=1
         GOTO 4
      ENDIF
      DD=UN/D
      RO=(XEL(1,1)+XEL(1,2))*UNDE
      R1=(XEL(2,1)+XEL(2,2))*UNDE
      SP=(XEL(1,2)-XEL(1,1))/D
      CP=(XEL(2,2)-XEL(2,1))/D
      SP2=SP*SP
      CP2=CP*CP
      SPCP=SP*CP
C           X FONCTION FORME NOEUD 2     RRRR RAYON
C           Y FONCTION FORME NOEUD 1     D LONGUEUR DD INVERSE LONGUEUR
      X=UNDE+UNDE*T(IGAU)
      Y=UNDE-UNDE*T(IGAU)
      RRRR=RO+UNDE*D*SP*T(IGAU)
C     ---------------------------------
C
C     TEST D'ERREUR
C
      IF(IFOU.GE.0) THEN
         IF(ABS(RRRR/D).LE.1.D-03) THEN
            IARR=2
            GOTO 4
         ENDIF
      ENDIF
C
C     ---------------------------------CALCULS
C
      IF(IFOU.LT.0) RRRR =UN
      U=X/RRRR
      V=Y/RRRR
C
C          AXISYMETRIQUE DEFORMATIONS PLANES CONTRAINTES PLANES
C
      IF(IFOU.LE.0) THEN
C
C          EPSILON S S
C
         B(1,1)=-SP*DD
         B(1,2)=-CP*DD
         B(1,4)=-B(1,1)
         B(1,5)=-B(1,2)
C
C            KSI S S
C
         AUX = SIX*T(IGAU)*DD*DD
         B(3,1)= CP*AUX
         B(3,2)=-SP*AUX
         B(3,3)=(QUTR-SIX*X)*DD
         B(3,4)=-B(3,1)
         B(3,5)=-B(3,2)
         B(3,6)=(DEUX-SIX*X)*DD
C
C          MODIFICATION DEFORMATION PLANE GENERALISEE
C
         IF (IFOU.EQ.-3) THEN
            RRRX=RO+UNDE*D*SP*T(IGAU)
            RRRY=R1+UNDE*D*CP*T(IGAU)
            B(2,7 )=UN
            B(2,8)=XDPGE-RRRX
            B(2,9)=RRRY-YDPGE
            B(4,8)=CP
            B(4,9)=SP
         ENDIF
C
C           AXISYMETRIQUE
C
         IF(IFOU.EQ.0) THEN
C
C          EPSILON  THETA  THETA
C
            B(2,1)= V*(SP2+CP2*Y*(UN+DEUX*X))
            B(2,2)= SPCP*U*Y*T(IGAU)
            B(2,3)=-D*CP*X*Y*V
            B(2,4)= U*(SP2+CP2*X*(TRS-DEUX*X))
            B(2,5)=-B(2,2)
            B(2,6)= D*CP*X*Y*U
C
C          KSI     THETA  THETA
C
            AUX    =-SIX*U*Y*SP*DD
            B(4,1)= CP*AUX
            B(4,2)=-SP*AUX
            B(4,3)= SP*V*(TRS*X-UN)
            B(4,4)=-CP*AUX
            B(4,5)= SP*AUX
            B(4,6)=-SP*U*(TRS*X-DEUX)
         ENDIF
C
C          FOURIER
C
      ELSE IF(IFOU.GT.0) THEN
         AN=DBLE(NN)
C
C               EPSILON S S
C
         B(1,1)=-SP*DD
         B(1,2)=-CP*DD
         B(1,5)=-B(1,1)
         B(1,6)=-B(1,2)
C
C               EPSILON THETA THETA
C
         B(2,1)= V*(SP2+CP2*Y*(UN+DEUX*X))
         B(2,2)= SPCP*U*Y*T(IGAU)
         B(2,3)=+AN*V
         B(2,4)=-D*CP*X*Y*V
         B(2,5)= U*(SP2+CP2*X*(TRS-DEUX*X))
         B(2,6)=-B(2,2)
         B(2,7)=+AN*U
         B(2,8)= D*CP*X*Y*U
C
C              EPSILON S  THETA
C
         B(3,1)=-SP*B(2,3)
         B(3,2)=-CP*B(2,3)
         B(3,3)=-DD-V*SP
         B(3,5)=-SP*B(2,7)
         B(3,6)=-CP*B(2,7)
         B(3,7)=+DD-U*SP
C
C           KSI  S S
C
         AUX= SIX*T(IGAU)*DD*DD
         B(4,1)= CP*AUX
         B(4,2)=-SP*AUX
         B(4,4)=(QUTR-SIX*X)*DD
         B(4,5)=-B(4,1)
         B(4,6)=-B(4,2)
         B(4,8)=(DEUX-SIX*X)*DD
C
C          KSI THETA THETA
C
         AUX1=SIX*U*Y*SP*DD
         AUX=(UN+DEUX*X)*(AN*V)**2+AUX1
         B(5,1)=-CP*AUX
         B(5,2)= SP*AUX
         B(5,3)=-B(2,3)*CP/RRRR
         B(5,4)= D*X*(AN*V)**2-SP*V*(UN-TRS*X)
         AUX=(TRS-DEUX*X)*(AN*U)**2-AUX1
         B(5,5)=-CP*AUX
         B(5,6)= SP*AUX
         B(5,7)=-B(2,7)*CP/RRRR
         B(5,8)=-D*Y*(AN*U)**2-SP*U*(TRS*X-DEUX)
C
C         KSI S  THETA
C
         AUX1=DOUZ*U*Y*AN*DD
         AUX=-AN*SP*(DEUX+QUTR*X)*V*V-AUX1
         B(6,1)=-CP*AUX
         B(6,2)= SP*AUX
         B(6,3)= DEUX*(CP/D+SPCP*V)/RRRR
         B(6,4)= AN*V*((DEUX-SIX*X)-DEUX*D*X*V*SP)
         AUX=AN*SP*(QUTR*X-SIX)*U*U+AUX1
         B(6,5)=-CP*AUX
         B(6,6)= SP*AUX
         B(6,7)= DEUX*(-CP/D+SPCP*U)/RRRR
         B(6,8)= AN*U*((SIX*X-QUTR)+DEUX*D*X*V*SP)
      ENDIF
      IF(IFOU.EQ.0.OR.(IFOU.EQ.1.AND.NN.EQ.0)) THEN
         DJAC=D*UNDE*P(IGAU)*RRRR*2*XPI
      ELSEIF(IFOU.EQ.1.AND.NN.NE.0) THEN
         DJAC=D*UNDE*P(IGAU)*RRRR*XPI
      ELSE
         DJAC=D*UNDE*P(IGAU)*RRRR*DIM3
      ENDIF
*
*  ON MODIFIE LA MATRICE B EN CAS D'EXCENTREMENT
*
      NSTRS2 = NSTRS / 2
      DO 88 IJL=1,NSTRS2
         DO 881 IJC=1,LRE
            B(IJL,IJC)=B(IJL,IJC)+EXCEN*B(IJL+NSTRS2,IJC)
 881     CONTINUE
 88   CONTINUE
*
 4    CONTINUE
      RETURN
      END