Sources Esope | Cast3M

cq2ksg
C CQ2KSG    SOURCE    CHAT      05/01/12    22:26:43     5004
      SUBROUTINE CQ2KSG(COO,EP,DIM3,IFOU,AN,NBPGAU,XSTRS,RAYON,
     .     RAY,TGAUS,CGAUS,ALF11,ALF12,ALF21,ALF22,A4,S3,AM,A1,A,LRE,REL
     $     )
*-----------------------------------------------------------------------
C     CALCUL DE LA MATRICE KSIGMA POUR UNE COQUE A 2 NOEUDS
C     FEVRIER 86 ; X. DE MAZANCOURT
C     PROGRAMME REPRIS A 90% SUR LE COQFLA DE INCA
*-----------------------------------------------------------------------
C     ENTREES:
C     COO:COORDONNEES DE L'ELEMENT
C     EP:EPAISSEUR DE LA COQUE
C     IFOU = VALEUR DE IFOUR DE CCOPTIO
C     AN:NUMERO DE L'HARMONIQUE DE FOURIER
C     XSTRS:CONTRAINTES
C     NBPGAU:NOMBRE DE POINTS DE GAUSS
C     TGAUS:POSITION DES POINTS DE GAUSS
C     CGAUS:POIDS DES POINTS DE GAUSS
C     RAYON,RAY,ALF:TABLEAUX DE TRAVAIL
C     LRE = TAILLE DE LA MATRICE KSIGMA
C     SORTIES:
C     REL:MATRICE KSIGMA
*-----------------------------------------------------------------------
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
      PARAMETER(UN=1.D0,UNDEMI=.5D0,UNQUAR=.25D0)
      DIMENSION A4(8,*),S3(8,*),AM(8,*),A1(8,*)
      DIMENSION COO(3,*),A(8,*),TGAUS(*),CGAUS(*)
      DIMENSION ALF11(*),ALF12(*),ALF21(*),ALF22(*),RAYON(*),REL(LRE,*)
      DIMENSION RAY(*),XSTRS(*),II(6)
      DATA II/1,2,4,5,6,8/
C     --------------------------------PRELIMINAIRES
      CALL ZERO(REL,LRE,LRE)
      RAYON(1)=COO(1,1)
      RAYON(3)=COO(1,2)
      RAYON(2)=(RAYON(1)+RAYON(3))*UNDEMI
      DR=RAYON(3)-RAYON(1)
      DZ=COO(2,2)-COO(2,1)
      DS=SQRT(DR*DR+DZ*DZ)
      COSP=DR/DS
      SINP=DZ/DS
      XS2=UNDEMI*DS
      XS4=UNQUAR*DS
C     --------------------------------
      CALL RESO1K(A4,DS,COSP,SINP,S3,A)
      CALL ZERO(ALF12,5,1)
      CALL ZERO(ALF22,7,1)
C     --------------------------------BOUCLE SUR LES POINTS DE GAUSS
      DO 9 IN=1,NBPGAU
C         ID1=1+(IN-1)*6
C         ID2=2+(IN-1)*6
         RAY(IN)=RAYON(2)+XS2*COSP*TGAUS(IN)
         X=TGAUS(IN)*XS2
         X2=X*X
         X3=X2*X
         X4=X3*X
C         ANS=EP*XSTRS(ID1)
C         ANT=EP*XSTRS(ID2)
C
C         ON PEUT FAIRE VARIER LES CONTRAINTES AVEC LE POINT DE GAUSS
C         MAIS ON A CHOISI DE PRENDRE CELLES DU 2} POINT DE GAUSS...
C
         IF (IFOU.EQ.1) THEN
            ANS=EP*XSTRS(7)
            ANT=EP*XSTRS(8)
         ELSE
            ANS=EP*XSTRS(5)
            ANT=EP*XSTRS(6)
         ENDIF
         ALF11(1)=UN
         ALF11(2)=X
         ALF11(3)=X2
         ALF11(4)=X3
         ALF11(5)=X4
         DO 12 NC=1,5
            ALF21(NC)=ALF11(NC)
 12      CONTINUE
         ALF21(6)=X4*X
         ALF21(7)=X4*X2
         RRRR=RAY(IN)
         IF(IFOU.EQ.-3) RRRR=UN
*
*   LE CAS DES DEF. PLANES GENERALISEES SERA A TESTER ||||
*
         IF(IFOU.EQ.-1.OR.IFOU.EQ.-2) THEN
            ANS=ANS*XS4*CGAUS(IN)
            ANT=0.D0
         ELSE
            ANS=ANS*RRRR*XS4*CGAUS(IN)
            ANT=(ANT*XS4*CGAUS(IN))/RRRR
         ENDIF
         IF (IFOU.EQ.0.OR.(IFOU.EQ.1.AND.AN.EQ.0)) THEN
            ANS=ANS*2*XPI
            ANT=ANT*2*XPI
         ELSEIF (IFOU.EQ.1.AND.AN.NE.0) THEN
            ANS=ANS*XPI
            ANT=ANT*XPI
         ENDIF
         ANS=ANS*DIM3
         ANT=ANT*DIM3
         CALL AMUL1K(ALF11,5,ANS,ALF12)
         CALL AMUL1K(ALF21,7,ANT,ALF22)
C     ---------------------------------------
 9    CONTINUE
C
C     TRANSFERT DE ALF12 (RESP.ALF22) DANS ALF11 (RESP.ALF21)
C
      CALL SHIF1K(ALF12,ALF11,5)
      CALL SHIF1K(ALF22,ALF21,7)
C
C     CALCUL DE A1,TRANSPOSEE DE A4
C
      DO 15 NC=1,8
         DO 13 NCC=1,8
            A1(NCC,NC)=A4(NC,NCC)
 13      CONTINUE
 15   CONTINUE
      CALL ZERO(A,8,8)
      PROD1=COSP*COSP
      PROD2=SINP*SINP
      PROD3=SINP*COSP
      A(1,1)=PROD1*ALF21(1)
      A(1,2)=PROD1*ALF21(2)
      A(1,5)=PROD3*ALF21(1)
      A(1,6)=PROD3*ALF21(2)
      A(1,7)=PROD3*ALF21(3)
      A(1,8)=PROD3*ALF21(4)
      A(2,1)=A(1,2)
      A(2,2)=ALF11(1)+PROD1*ALF21(3)
      A(2,5)=PROD3*ALF21(2)
      A(2,6)=PROD3*ALF21(3)
      A(2,7)=PROD3*ALF21(4)
      A(2,8)=PROD3*ALF21(5)
      A(3,3)=ALF21(1)
      A(3,4)=ALF21(2)
      A(4,3)=A(3,4)
      A(4,4)=ALF11(1)+ALF21(3)
      A(5,1)=A(1,5)
      A(5,2)=A(2,5)
      A(5,5)=PROD2*ALF21(1)
      A(5,6)=PROD2*ALF21(2)
      A(5,7)=PROD2*ALF21(3)
      A(5,8)=PROD2*ALF21(4)
      A(6,1)=A(1,6)
      A(6,2)=A(2,6)
      A(6,5)=A(5,6)
      A(6,6)=ALF11(1)+A(5,7)
      A(6,7)=ALF11(2)*2.D0+A(5,8)
      A(6,8)=ALF11(3)*3.D0+PROD2*ALF21(5)
      A(7,1)=A(1,7)
      A(7,2)=A(2,7)
      A(7,5)=A(5,7)
      A(7,6)=A(6,7)
      A(7,7)=ALF11(3)*4.D0+PROD2*ALF21(5)
      A(7,8)=ALF11(4)*6.D0+PROD2*ALF21(6)
      A(8,1)=A(1,8)
      A(8,2)=A(2,8)
      A(8,5)=A(5,8)
      A(8,6)=A(6,8)
      A(8,7)=A(7,8)
      A(8,8)=ALF11(5)*9.D0+PROD2*ALF21(7)
 20   CONTINUE
      CALL MULMAT (S3,A,A4,8,8,8)
      CALL MULMAT (A,A1,S3,8,8,8)
      CALL ZERO(AM,8,8)
      AM(1,1)=ALF21(1)
      AM(3,3)=ALF21(1)
      AM(5,5)=ALF21(1)
      AM(1,2)=ALF21(2)
      AM(2,1)=ALF21(2)
      AM(3,4)=ALF21(2)
      AM(4,3)=ALF21(2)
      AM(5,6)=ALF21(2)
      AM(6,5)=ALF21(2)
      AM(2,2)=ALF21(3)
      AM(4,4)=ALF21(3)
      AM(6,6)=ALF21(3)
      AM(5,7)=ALF21(3)
      AM(7,5)=ALF21(3)
      AM(5,8)=ALF21(4)
      AM(8,5)=ALF21(4)
      AM(6,7)=ALF21(4)
      AM(7,6)=ALF21(4)
      AM(7,7)=ALF21(5)
      AM(6,8)=ALF21(5)
      AM(8,6)=ALF21(5)
      AM(7,8)=ALF21(6)
      AM(8,7)=ALF21(6)
      AM(8,8)=ALF21(7)
      CALL MULMAT (S3,AM,A4,8,8,8)
      CALL MULMAT (AM,A1,S3,8,8,8)
      AN2=AN*AN
      DO 30 NC=1,8
         DO 31 NCC=1,8
            A(NCC,NC)=A(NCC,NC)+AN2*AM(NCC,NC)
 31      CONTINUE
 30   CONTINUE
      PROD1=2.D0*COSP
      PROD2=2.D0*SINP
      CALL ZERO(AM,8,8)
      AM(1,3)=PROD1*ALF21(1)
      AM(1,4)=PROD1*ALF21(2)
      AM(2,3)=AM(1,4)
      AM(2,4)=PROD1*ALF21(3)
      AM(3,1)=AM(1,3)
      AM(3,2)=AM(2,3)
      AM(3,5)=PROD2*ALF21(1)
      AM(3,6)=PROD2*ALF21(2)
      AM(3,7)=PROD2*ALF21(3)
      AM(3,8)=PROD2*ALF21(4)
      AM(4,1)=AM(1,4)
      AM(4,2)=AM(2,4)
      AM(4,5)=AM(3,6)
      AM(4,6)=AM(3,7)
      AM(4,7)=AM(3,8)
      AM(7,4)=AM(4,7)
      AM(4,8)=PROD2*ALF21(5)
      AM(5,3)=AM(3,5)
      AM(5,4)=AM(4,5)
      AM(6,3)=AM(3,6)
      AM(6,4)=AM(4,6)
      AM(7,3)=AM(3,7)
      AM(8,3)=AM(3,8)
      AM(8,4)=AM(4,8)
 19   CONTINUE
      CALL MULMAT (S3,AM,A4,8,8,8)
      CALL MULMAT (AM,A1,S3,8,8,8)
C
C     CALCUL FINAL DE LA MATRICE KSIGMA
C
      IF(IFOU.EQ.1) THEN
         DO 32 NC=1,8
            DO 33 NCC=1,8
               REL(NCC,NC) =(A(NCC,NC)+AN*AM(NCC,NC))*2.D0
 33         CONTINUE
 32      CONTINUE
      ELSE IF(IFOU.EQ.0.OR.IFOU.EQ.-1.OR.IFOU.EQ.-2)  THEN
         DO 82 NC=1,6
            IC=II(NC)
            DO 83 NCC=1,6
               ICC=II(NCC)
               REL(NCC,NC) =(A(ICC,IC)+AN*AM(ICC,IC))*2.D0
 83         CONTINUE
 82      CONTINUE
      ENDIF
      RETURN
      END