Sources Esope | Cast3M

coq2ma
C COQ2MA    SOURCE    CHAT      06/03/29    21:17:32     5360
      SUBROUTINE COQ2MA(Q,E,DIM3,RHO,IMATMA,IFOUR,NIFOU,LRE,
     +     AMM,IARR,XDPGE,YDPGE)
C
C***********************************************************************
C          MATRICE DE MASSE D'UN ELEMENT COQUE A 2 NOEUDS
C          CODE X.DE MAZANCOURT   DECEMBRE 85
C***********************************************************************
C ENTREES
C     Q(3,2) : COORDONNEES DES 2 NOEUDS , R PUIS Z
C     E : EPAISSEUR ( CONSTANTE )
C     RHO : MASSE VOLUMIQUE
C     IMATMA=1-MATRICE MASSE- OU 0-VECTEUR MASSE-(CAS DE MATR DIAGONALE)
C            AVANT L'APPEL IL EST MIS A 1
C     IFOUR=IFOUR DE CCOPTIO
C     NIFOU=NIFOUR DE CCOPTIO
C     LRE  = TAILEE DE LA MATRICE MASSE 6 SI IFOU.LE.0 8 SI IFOU.GT.0
C SORTIE
C    AMM(LRE,LRE) = MATRICE DE MASSE
C     IARR:INDICATEUR D'ERREUR,VAUT ZERO NORMALEMENT
C***********************************************************************
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
      DIMENSION Q(3,*),AM(9,9),AMM(LRE,*)
      DIMENSION II(6)
      DATA II/1,2,4,5,6,8/
C
C     ----------   INITIALISATION
      CALL ZERO(AM,9,9)
      CALL ZERO(AMM,LRE,LRE)
C
      IARR=0
      R=.5D0*(Q(1,1)+Q(1,2))
      Y=.5D0*(Q(2,1)+Q(2,2))
      XD=XDPGE-R
      YD=Y-YDPGE
      SP=Q(1,2)-Q(1,1)
      CP=Q(2,2)-Q(2,1)
      D2=CP*CP+SP*SP
      D=SQRT(D2)
      E2=E*E
C
C     -----------------------------TESTS D'ERREUR
C
      IF(D.EQ.0) THEN
         IARR=1
         GOTO 60
      ENDIF
      IF(IFOUR.GE.0) THEN
         IF(ABS(R/D).LE.1.D-03) THEN
            IARR=2
            GOTO 60
         ENDIF
      ELSE
         R=1.D0
      ENDIF
C
C     ----------   NORME
C ********** CASTEM20 SOMME(COS**2)=1
C
      COEF=RHO*E*D*DIM3
      IF(IMATMA.NE.0) GOTO 20
C     ----------   VECTEUR MASSE
      RM=.5D0*COEF*R
      IF(IFOUR.EQ.0.OR.(IFOUR.EQ.1.AND.NIFOU.EQ.0)) THEN
         RM=2*XPI*RM
      ELSEIF(IFOUR.EQ.1.AND.NIFOU.NE.0) THEN
         RM=XPI*RM
      ENDIF
      DO 15 I=1,2
         K=4*(I-1)
         DO 151 J=1,3
            K=K+1
            AM(K,K)=RM
 151     CONTINUE
 15   CONTINUE
      RETURN
C
C     ----------   PRELIMINAIRES
C
 20   CONTINUE
      A=COEF*R/420.D0
      B=COEF*SP/840.D0
      IF(IFOUR.LT.0) B=0.D0
      IF(IFOUR.EQ.0.OR.(IFOUR.EQ.1.AND.NIFOU.EQ.0)) THEN
         A=2*XPI*A
         B=2*XPI*B
      ELSEIF(IFOUR.EQ.1.AND.NIFOU.NE.0) THEN
         A=XPI*A
         B=XPI*B
      ENDIF
      SP=SP/D
      CP=CP/D
      SP2=SP*SP
      CP2=CP*CP
      SPCP=SP*CP
      DSP=D*SP
      DCP=D*CP
      ASP2=A*SP2
      ACP2=A*CP2
      BSP2=B*SP2
      BCP2=B*CP2
      ASPCP=A*SPCP
      BSPCP=B*SPCP
      AUX=22.D0*A-8.D0*B
      AUX1=13.D0*A-B
C     ----------   CALCULS
      AM(1,1)=156.D0*ACP2-84.D0*BCP2+140.D0*ASP2-70.D0*BSP2
      AM(1,2)=14.D0*BSPCP-16.D0*ASPCP
      AM(1,4)=-DCP*AUX
      AM(1,5)=54.D0*ACP2+70.D0*ASP2
      AM(1,6)=16.D0*ASPCP
      AM(1,8)=DCP*AUX1
      AM(2,2)=156.D0*ASP2-84.D0*BSP2+140.D0*ACP2-70.D0*BCP2
      AM(2,4)=DSP*AUX
      AM(2,5)=AM(1,6)
      AM(2,6)=54.D0*ASP2+70.D0*ACP2
      AM(2,8)=-DSP*AUX1
      AM(3,3)=140.D0*A-70.D0*B
      AM(3,7)=70.D0*A
      AM(4,4)=D2*(4.D0*A-B)
      AUX=22.D0*A+8.D0*B
      AUX1=13.D0*A+B
      AM(4,5)=-DCP*AUX1
      AM(4,6)=DSP*AUX1
      AM(4,8)=-3.D0*D2*A
      AM(5,5)=156.D0*ACP2+84.D0*BCP2+140.D0*ASP2+70.D0*BSP2
      AM(5,6)=-16.D0*ASPCP-14.D0*BSPCP
      AM(5,8)=DCP*AUX
      AM(6,6)=156.D0*ASP2+84.D0*BSP2+140.D0*ACP2+70.D0*BCP2
      AM(6,8)=-DSP*AUX
      AM(7,7)=140.D0*A+70.D0*B
      AM(8,8)=D2*(4.D0*A+B)
C     ----------   ON COMPLETE PAR SYMETRIE
      DO 50 I=2,8
         JFIN=I-1
         DO 501 J=1,JFIN
            AM(I,J)=AM(J,I)
 501     CONTINUE
 50   CONTINUE
 60   CONTINUE
      IF(IFOUR.GT.0) THEN
         DO  500 IA= 1,8
            DO  5001 IB= 1,8
               AMM(IA,IB)=AM(IA,IB)
 5001       CONTINUE
 500     CONTINUE
      ELSE IF(IFOUR.LE.0) THEN
         DO  600 IA= 1,6
            DO  6001 IB= 1,6
               AMM(IA,IB)=AM(II(IA),II(IB))
 6001       CONTINUE
 600     CONTINUE
      ENDIF
      IF(IFOUR.EQ.-3)THEN
         AMM(7,7)=COEF
         AMM(7,8)=COEF*XD
         AMM(7,9)=COEF*YD
         AMM(8,7)=AMM(7,8)
         AMM(8,8)=COEF*((D2*CP2/12.D0)+(E2*SP2/12.D0)+(XD*XD))
         AMM(8,9)=COEF*(SPCP*((E2/12.D0)-(D2/12.D0))+(XD*YD))
         AMM(9,7)=AMM(7,9)
         AMM(9,8)=AMM(8,9)
         AMM(9,9)=COEF*((E2*CP2/12.D0)+(D2*SP2/12.D0)+(YD*YD))
      ENDIF
      RETURN
      END