Sources Esope | Cast3M

bmatst
C BMATST    SOURCE    PV        20/09/16    21:15:05     10713          
 
C=======================================================================
C=                          B M A T S T                                =
C=                          -----------                                =
C=  Fonction :                                                         =
C=  ----------                                                         =
C=   Calcul de la matrice B reliant les deformations en un point d'un  =
C=   element fini aux ddls de deplacement aux noeuds de cet element    =
C=   Le jacobien est egalement evalue au point de Gauss pour verifier  =
C=   ulterieurement si l'element fini n'est pas trop distordu.         =
C=                                                                     =
C=  Parametres :  (E)=Entree  (S)=Sortie                               =
C=  ------------                                                       =
C=   iGau     (E)   Numero du point de Gauss considere                 =
C=   NBPGAU   (E)   Nombre de points de Gauss de l'element fini        =
C=   POIGAU,QSIGAU   (E)   |  Poids et coordonnees de reference des    =
C=   ETAGAU,DZEGAU   (E)   |  differents points de Gauss de l'element  =
C=   ITEL     (E)   Type de l'element fini (cf. NOMTP dans bdata.eso)  =
C=   MFR      (E)   Formulation associee a l'element fini              =
C=   NBNO     (E)   Nombre de noeuds de l'element fini                 =
C=   LRE      (E)   Nombre de DDL de l'element fini                    =
C=   IFOU     (E)   Mode de calcul utilise (cf. IFOUR dans CCOPTIO)    =
C=   NST      (E)   Nombre de composantes de deformations              =
C=   NN       (E)   Numero du mode de Fourier (si IFOU=1)              =
C=   DIM3     (E)   Epaisseur dans l'hypothese des contraintes planes  =
C=   XEL      (E)   Coordonnees des noeuds de l'element fini etudie    =
C=   SHPTOT   (E)   Fonctions de forme et leurs derivees               =
C=   SHP      (S)   Fonctions de forme et leurs derivees actuelles     =
C=   BGENE    (S)   Matrice de gradients (B) calcule au point de Gauss =
C=   DJac     (S)   Jacobien au point de Gauss etudie                  =
C=   XDPGE,YDPGE   (E)   Coordonnees du noeud support (mode PLAN GENE) =
C=                                                                     =
C=  Remarque :                                                         =
C=  ----------                                                         =
C=   Lors de l'entree dans le sous-programme, SHPTOT(2 a 4,*) contient =
C=   les DERIVEES des fonctions de forme par rapport aux coordonnees   =
C=   de REFERENCE Qsi,Eta,Dzeta.                                       =
C=   En sortie du sous-programme, SHP(2 a 4,*) contient les DERIVEES   =
C=   des fonctions de FORME par rapport aux coordonnees REELLES x,y,z. =
C=   Pour les elements finis massifs incompressibles, la matrice B est =
C=   modifiee ulterieurement par appel au sous-programme BBAR.         =
C=======================================================================
 
      SUBROUTINE BMATST (iGau,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU,
     &                   ITEL,MFR,NBNO,LRE,IFOU,NST,NN,DIM3,XEL,
     &                   SHPTOT,SHP,BGENE,DJAC,XDPGE,YDPGE)
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
 
-INC CCREEL
 
      DIMENSION XEL(3,*),BGENE(NST,*),SHP(6,*),SHPTOT(6,NBNO,*),
     &          POIGAU(*),QSIGAU(*),ETAGAU(*),DZEGAU(*)
 
      RR=XZero
** semble inutile
**    CALL ZERO(BGENE,NST,LRE)
 
      IFOR=IFOU+4
      GOTO ( 5,10,10,20,30,40,50,50,50,50,50,50,50,50,50,
     &      60,60,60,70),IFOR
      GOTO 100
 
C ===
C  1 - Elements MASSIFS 2D PLAN GENE
C ===
 5    XXX=XZero
      YYY=XZero
      DO i=1,NBNO
        SHP(1,i)=SHPTOT(1,i,iGau)
        SHP(2,i)=SHPTOT(2,i,iGau)
        SHP(3,i)=SHPTOT(3,i,iGau)
        XXX=XXX+SHP(1,i)*XEL(1,i)
        YYY=YYY+SHP(1,i)*XEL(2,i)
      ENDDO
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NN,2,DIM3,RR,DJAC)
      K=1
      DO i=1,NBNO
        BGENE(1,K)=SHP(2,i)
        BGENE(2,K+1)=SHP(3,i)
        BGENE(4,K)=SHP(3,i)
        BGENE(4,K+1)=SHP(2,i)
        K=K+2
      ENDDO
      BGENE(3,K)=1.D0
      BGENE(3,K+1)=XDPGE-XXX
      BGENE(3,K+2)=YYY-YDPGE
      GOTO 100
C ===
C  2 - Elements MASSIFS 2D PLAN DEFO ou PLAN CONT
C ===
 10   DO i=1,NBNO
        SHP(1,i)=SHPTOT(1,i,iGau)
        SHP(2,i)=SHPTOT(2,i,iGau)
        SHP(3,i)=SHPTOT(3,i,iGau)
      ENDDO
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NN,2,DIM3,RR,DJAC)
      K=1
      DO i=1,NBNO
        BGENE(1,K)=SHP(2,i)
        BGENE(2,K+1)=SHP(3,i)
        BGENE(4,K)=SHP(3,i)
        BGENE(4,K+1)=SHP(2,i)
        K=K+2
      ENDDO
      GOTO 100
C ===
C  3 - Elements MASSIFS 2D AXISymetrique
C ===
 20   DO i=1,NBNO
        SHP(1,i)=SHPTOT(1,i,iGau)
        SHP(2,i)=SHPTOT(2,i,iGau)
        SHP(3,i)=SHPTOT(3,i,iGau)
      ENDDO
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NN,2,DIM3,RR,DJAC)
      K=1
      DO i=1,NBNO
        BGENE(1,K)=SHP(2,i)
        BGENE(2,K+1)=SHP(3,i)
        BGENE(3,K)=SHP(1,i)/RR
        BGENE(4,K)=SHP(3,i)
        BGENE(4,K+1)=SHP(2,i)
        K=K+2
      ENDDO
      GOTO 100
C ===
C  4 - Elements MASSIFS 2D FOURIER
C ===
 30   DO i=1,NBNO
        SHP(1,i)=SHPTOT(1,i,iGau)
        SHP(2,i)=SHPTOT(2,i,iGau)
        SHP(3,i)=SHPTOT(3,i,iGau)
      ENDDO
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NN,2,DIM3,RR,DJAC)
      XNSUR=DBLE(NN)/RR
      K=1
      DO i=1,NBNO
        BGENE(1,K)=SHP(2,i)
        BGENE(2,K+1)=SHP(3,i)
        BGENE(3,K)=SHP(1,i)/RR
        BGENE(3,K+2)=SHP(1,i)*XNSUR
        BGENE(4,K)=SHP(3,i)
        BGENE(4,K+1)=SHP(2,i)
        BGENE(5,K)=-SHP(1,i)*XNSUR
        BGENE(5,K+2)=SHP(2,i)-SHP(1,i)/RR
        BGENE(6,K+1)=-SHP(1,i)*XNSUR
        BGENE(6,K+2)=SHP(3,i)
        K=K+3
      ENDDO
      GOTO  100
C ===
C  6 - Elements MASSIFS 3D
C ===
 40   DO i=1,NBNO
        SHP(1,i)=SHPTOT(1,i,iGau)
        SHP(2,i)=SHPTOT(2,i,iGau)
        SHP(3,i)=SHPTOT(3,i,iGau)
        SHP(4,i)=SHPTOT(4,i,iGau)
      ENDDO
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NN,3,DIM3,RR,DJAC)
      K=1
      DO i=1,NBNO
        BGENE(1,K)=SHP(2,i)
        BGENE(2,K+1)=SHP(3,i)
        BGENE(3,K+2)=SHP(4,i)
        BGENE(4,K)=SHP(3,i)
        BGENE(4,K+1)=SHP(2,i)
        BGENE(5,K)=SHP(4,i)
        BGENE(5,K+2)=SHP(2,i)
        BGENE(6,K+1)=SHP(4,i)
        BGENE(6,K+2)=SHP(3,i)
        K=K+3
      ENDDO
      GOTO 100
C ===
C  7 - Elements MASSIFS 1D
C ===
C= 7.1 - Modes 1D PLAN
 50   DO i=1,NBNO
        SHP(1,i)=SHPTOT(1,i,iGau)
        SHP(2,i)=SHPTOT(2,i,iGau)
      ENDDO
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NN,1,DIM3,RR,DJAC)
      DO i=1,NBNO
        BGENE(1,i)=SHP(2,i)
      ENDDO
      IF (IFOU.EQ.7.OR.IFOU.EQ.8) THEN
        BGENE(2,i)=1.
      ELSE IF (IFOU.EQ.9.OR.IFOU.EQ.10) THEN
        BGENE(3,i)=1.
      ELSE IF (IFOU.EQ.11) THEN
        BGENE(2,i)=1.
        BGENE(3,i+1)=1.
      ENDIF
      GOTO 100
C= 7.2 - Modes 1D AXIS
 60   DO i=1,NBNO
        SHP(1,i)=SHPTOT(1,i,iGau)
        SHP(2,i)=SHPTOT(2,i,iGau)
      ENDDO
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NN,1,DIM3,RR,DJAC)
      DO i=1,NBNO
        BGENE(1,i)=SHP(2,i)
        BGENE(3,i)=SHP(1,i)/RR
      ENDDO
      IF (IFOU.EQ.14) BGENE(2,i)=1.
      GOTO 100
C= 7.3 - Mode 1D SPHE
 70   DO i=1,NBNO
        SHP(1,i)=SHPTOT(1,i,iGau)
        SHP(2,i)=SHPTOT(2,i,iGau)
      ENDDO
      CALL DEVOLU(XEL,SHP,MFR,NBNO,IFOU,NN,1,DIM3,RR,DJAC)
      DO i=1,NBNO
        BGENE(1,i)=SHP(2,i)
        ZZ=SHP(1,i)/RR
        BGENE(2,i)=ZZ
        BGENE(3,i)=ZZ
      ENDDO
      GOTO 100
 100  CONTINUE
 
      RETURN
      END