Esope source files

   1 : C COQ3M1    SOURCE    CHAT      05/01/12    22:21:24     5004           
   2 :       SUBROUTINE COQ3M1(XE,RHO,XEL,BPSS,BB,REWO,ILUMP)
   3 : c=======================================================================
   4 : c
   5 : c      calcule la matrice de masse de l element coq3
   6 : c
   7 : c   entree
   8 : c      xe(3,3)=coodonnees de l element
   9 : c      rho    =masse volumique * epaisseur
  10 : c      ilump = 1 si l'opérateur LUMP est appelé , 0 sinon
  11 : c   travail
  12 : c      xel(3,3)   =coordonnees locales de l element
  13 : c      bpss(3,3)  =matrice de passage
  14 : c      bb(9)      =stocke les fonctions de formes de flexion
  15 : c   sorties
  16 : c      rewo(18,18)=matrice de masse    repere local puis global
  17 : c
  18 : c      code fevrier 85   ebersolt
  19 : c
  20 : c=======================================================================
  21 :       IMPLICIT INTEGER(I-N)
  22 :       IMPLICIT REAL*8(A-H,O-Z)
  23 :       DIMENSION XE(3,*),XEL(3,*),BPSS(3,*),REWO(18,*),BB(*)
  24 :       DIMENSION REWOLP(18,18)
  25 :       DIMENSION XX(3),YY(3)
  26 :       DATA XX/0.5D0,0.5D0,0.0D0/
  27 :       DATA YY/0.0D0,0.5D0,0.5D0/
  28 :       DATA UNSIX,UN12/.166666666666666666D0,.833333333333333333D-1/
  29 :       DATA UNTIER/.33333333333333333D0/
  30 :       DATA IZERO/0/
  31 : c
  32 : c     matrice de passage
  33 : c
  33 :       CALL VPAST(XE,BPSS)
  34 : c
  35 : c     coordonnees locales
  36 : c
  36 :       CALL VCORLC(XE,XEL,BPSS)
  37 : c
  38 : c     mise a 0 de la matrice de masse
  39 : c
  39 :       CALL ZERO(REWO,18,18)
  40 :       X21=XEL(1,2)-XEL(1,1)
  41 :       Y31=XEL(2,3)-XEL(2,1)
  42 :       SURF=X21*Y31*.5D0*RHO
  43 : c
  44 : c     termes diagonaux  de la membrane
  45 : c
  46 :       DO 100 IA=1,13,6
  47 :           REWO(IA  ,IA  )=UNSIX
  48 :           REWO(IA+1,IA+1)=UNSIX
  49 :           DO 110 IB=1,13,6
  50 :               IF(IB.EQ.IA) GOTO 110
  51 :               REWO(IA  ,IB  )=UN12
  52 :               REWO(IA+1,IB+1)=UN12
  53 :   110     CONTINUE
  54 :   100 CONTINUE
  55 : c
  56 : c     partie flexion   integration numerique
  57 : c
  58 : c     ia numero du point de gauss
  59 : c     ib numero d un noeud
  60 : c     ic numero d un noeud
  61 : c
  62 :       DO 200 IA=1,3
  62 :           CALL MFDKT(XX(IA),YY(IA),XEL,BB)
  63 :           DO 210 IB=1,3
  64 :               IBDERI=(IB-1)*6+2
  65 :               IBDEBB=(IB-1)*3
  66 :               DO 210 IC=1,3
  67 :                   ICDERI=(IC-1)*6+2
  68 :                   ICDEBB=(IC-1)*3
  69 :                   DO 220 ID=1,3
  70 :                    DO 220 IE=1,3
  71 :                     REWO(IBDERI+ID,ICDERI+IE)=REWO(IBDERI+ID,ICDERI+IE)
  72 :      1               +BB(IBDEBB+ID)*BB(ICDEBB+IE)*UNTIER
  73 :   220             CONTINUE
  74 :   210     CONTINUE
  75 :   200 CONTINUE
  76 : c
  77 : c     multiplication par rho*surf
  78 : c
  79 :       DO 300 IA=1,18
  80 :           DO 300 IB=1,18
  81 :               REWO(IA,IB)=REWO(IA,IB)*SURF
  82 :   300 CONTINUE
  83 : c
  84 : c    diagonalisation dans le cas de l'opérateur LUMP
  85 : c
  86 : c      REWO est rangé dans l'ordre i noeud x(UX UY UZ RX RY RZ) ....
  87 : c
  88 :       IF ( ILUMP .EQ. 1 ) THEN
  88 :         CALL LUMP3(REWO)
  89 :       ENDIF  
  90 : c          
  91 : c
  92 : c     changement de repere
  93 : c
  93 :       CALL TRANSK(REWO,BPSS,18,3,IZERO)
  94 : c
  95 :       RETURN
  96 :       END
  97 :  
  98 :