Sources Esope | Cast3M

bcoq4
C BCOQ4     SOURCE    CHAT      06/03/29    21:15:38     5360
      SUBROUTINE BCOQ4
     &          (IGAU,XEL,SHPTOT,SHP,BGENE,DJAC,EXCEN,NOPLAN,IARR,IFS)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
************************************************************************
*
*                             B C O Q 4
*                             ---------
*
* FONCTION:
* ---------
*
*    Calcul de la matrice "B" pour un COQ4 excentr{ ou non.
*
* PARAMETRES:   (E)=ENTREE   (S)=SORTIE
* -----------
*
*     IGAU    (E)  Num{ro du point d'int{gration.
*     XEL     (E)  Coordonn{es locales des noeuds de l'{l{ment.
*     SHPTOT  (E)  Fonctions de forme et d{riv{es dans l'espace de
*                  r{f{rence.
*                  SHPTOT(1,... = fonction,
*                  SHPTOT(2,... = d{riv{e par rapport @ Qsi,
*                  SHPTOT(3,... = d{riv{e par rapport @ Eta.
*     EXCEN   (E)  Excentrement.
*     NOPLAN  (E)  = 1 si {l{ment non plan,
*                  = 0 sinon.
*     IFS     (E)  = 1 si demande de Matrice "B" compl}te
*                  (pour SIGMA et EPSILON),
*                  = 0 sinon.
*     SHP     (S)  Fonctions de forme et d{riv{es dans l'espace
*                  g{om{trique.
*                  SHP(1,... = fonction,
*                  SHP(2,... = d{riv{e par rapport @ "X" local,
*                  SHP(3,... = d{riv{e par rapport @ "Y" local.
*                  - au point d'int{gration n.5 si IFS=1,
*                  - au point d'int{gration IGAU courant sinon.
*     DJAC    (S)  Jacobien au point d'int{gration IGAU.
*     BGENE   (S)  Matrice "B" au point d'int{gration IGAU.
*     IARR    (S)  Diff{rent de 0  si erreur.
*
      PARAMETER (LRE=24,NST=8,NBNO=4)
      REAL*8 XEL(3,*),BGENE(NST,*),SHP(6,*),SHPTOT(6,NBNO,*)
*
* CONSTANTES:
* -----------
*
*     NBNO   = Nombre de noeuds.
*     LRE    = Nombre de colonnes de la Matrice "B".
*     NST    = Nombre de composantes de contraintes.
*
* VARIABLES:
* ----------
*
*     FULL   = .TRUE. si on veut une Matrice "B" compl}te.
*     GAUSS  = .TRUE. si utilisation d'un point d'int{gration n.1 @ 4.
*     SHP7   = Fonctions d'interpolation mixte pour le cisaillement X-Z
*     SHP8   = Fonctions d'interpolation mixte pour le cisaillement Y-Z
*
      REAL*8 SHP7(4:5,NBNO),SHP8(4:5,NBNO)
      LOGICAL FULL,GAUSS
*
* MODE DE FONCTIONNEMENT:
* -----------------------
*
* POUR LA RIGIDITE:
*
*     Pour la membrane,
*     - int{gration r{duite (1 point) pour le cisaillement,
*     - int{gration r{duite de tout si {l{ment non plan.
*
*     Pour le cisaillement transverse, voir:
*     DONEA / LAMAIN (1987) - ISPRA
*     ({l{ment de type "A.N.S.")
*
*     Pour le couplage membrane-flexion (excentrement), int{gration
*     r{duite analogue @ membrane, non justifi{e au niveau de la
*     Rigidit{, mais indispensable pour le calcul pratique des forces
*     nodales @ partir des contraintes, par "BSIGMA".
*
* POUR LES CONTRAINTES:
*
*     Pour {viter des surcontraintes dues @ un "hourglass", il ne faut
*     pas calculer aux 4 points de Gauss des contraintes qui ont {t{
*     calcul{es au centre lors de la cr{ation de la Matrice de Rigidit{.
*     Par cons{quent, il faut prendre
*     - toujours      : le cisaillement de membrane au centre,
*     - quand NOPLAN=1: toutes les contraintes de membrane au centre.
*     (centre = point d'int{gration n.5)
*     Calcul au centre dans les memes conditions pour le couplage
*     membrane-flexion.
*
*     Cependant, pour que le champ de contraintes produit soit facile
*     @ lire par l'utilisateur et @ traiter par les op{rateurs tels
*     que:
*     - CHANGER CHPOINT,
*     - TRACER, ...
*     on inscrit aux points de Gauss 1 @ 4 les valeurs de contrainte
*     au point d'int{gration n.5 pour les parties concern{es.
*     De plus, on remplit compl}tement le champ de valeurs au point
*     d'int{gration n.5 .
*
* AUTEUR, DATE DE CREATION:
* -------------------------
*
*     Pascal MANIGOT     17 Juin 1991
*
************************************************************************
*
      IARR=0
      FULL = IFS .EQ. 1
      GAUSS = IGAU .LE. 4
*
*
*     Matrice jacobienne:
*     -------------------
*
      DO 110 I=1,NBNO
         SHP(1,I)=SHPTOT(1,I,IGAU)
         SHP(2,I)=SHPTOT(2,I,IGAU)
         SHP(3,I)=SHPTOT(3,I,IGAU)
  110    CONTINUE
*     END DO
*
      DXDQSI = 0.D0
      DXDETA = 0.D0
      DYDQSI = 0.D0
      DYDETA = 0.D0
      DO 120 I=1,NBNO
         DXDQSI=DXDQSI+SHP(2,I)*XEL(1,I)
         DXDETA=DXDETA+SHP(3,I)*XEL(1,I)
         DYDQSI=DYDQSI+SHP(2,I)*XEL(2,I)
         DYDETA=DYDETA+SHP(3,I)*XEL(2,I)
 120     CONTINUE
*     END DO
      DJAC=DXDQSI*DYDETA-DXDETA*DYDQSI
      IF(DJAC.LE. 0.D0) THEN
         IARR=1
         RETURN
      END IF
*
*     Fonctions de forme dans l'espace g{om{trique:
      DQSIDX= DYDETA/DJAC
      DQSIDY=-DXDETA/DJAC
      DETADX=-DYDQSI/DJAC
      DETADY= DXDQSI/DJAC
      DO 130 I=1,NBNO
         TT=SHP(2,I)*DQSIDX+SHP(3,I)*DETADX
         SHP(3,I)=SHP(2,I)*DQSIDY+SHP(3,I)*DETADY
         SHP(2,I)=TT
  130    CONTINUE
*     END DO
*
*
*     Matrice "B":
*     ------------
*
      CALL ZERO(BGENE,NST,LRE)
*
      IF (GAUSS .OR. FULL) THEN
*        Pr{paration pour cisaillement transverse:
*
         N1 = 1
         N2 = 2
         DX1 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2)
         DY1 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2)
         N1 = 3
         N2 = 4
         DX2 = SHPTOT(2,N1,IGAU)*XEL(1,N1)+SHPTOT(2,N2,IGAU)*XEL(1,N2)
         DY2 = SHPTOT(2,N1,IGAU)*XEL(2,N1)+SHPTOT(2,N2,IGAU)*XEL(2,N2)
         N1 = 1
         N2 = 4
         DX3 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2)
         DY3 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2)
         N1 = 2
         N2 = 3
         DX4 = SHPTOT(3,N1,IGAU)*XEL(1,N1)+SHPTOT(3,N2,IGAU)*XEL(1,N2)
         DY4 = SHPTOT(3,N1,IGAU)*XEL(2,N1)+SHPTOT(3,N2,IGAU)*XEL(2,N2)
*
         D2JAC = 2.D0 * DJAC
*
         SHP7(4,1) = (- DYDETA*DY1 + DYDQSI*DY3) / D2JAC
         SHP7(4,2) = (- DYDETA*DY1 + DYDQSI*DY4) / D2JAC
         SHP7(4,3) = (- DYDETA*DY2 + DYDQSI*DY4) / D2JAC
         SHP7(4,4) = (- DYDETA*DY2 + DYDQSI*DY3) / D2JAC
         SHP7(5,1) = (  DYDETA*DX1 - DYDQSI*DX3) / D2JAC
         SHP7(5,2) = (  DYDETA*DX1 - DYDQSI*DX4) / D2JAC
         SHP7(5,3) = (  DYDETA*DX2 - DYDQSI*DX4) / D2JAC
         SHP7(5,4) = (  DYDETA*DX2 - DYDQSI*DX3) / D2JAC
*
         SHP8(4,1) = (  DXDETA*DY1 - DXDQSI*DY3) / D2JAC
         SHP8(4,2) = (  DXDETA*DY1 - DXDQSI*DY4) / D2JAC
         SHP8(4,3) = (  DXDETA*DY2 - DXDQSI*DY4) / D2JAC
         SHP8(4,4) = (  DXDETA*DY2 - DXDQSI*DY3) / D2JAC
         SHP8(5,1) = (- DXDETA*DX1 + DXDQSI*DX3) / D2JAC
         SHP8(5,2) = (- DXDETA*DX1 + DXDQSI*DX4) / D2JAC
         SHP8(5,3) = (- DXDETA*DX2 + DXDQSI*DX4) / D2JAC
         SHP8(5,4) = (- DXDETA*DX2 + DXDQSI*DX3) / D2JAC
*
      END IF
*
      DO 300 I=1,NBNO
         K = 6*(I-1)
*
         IF (GAUSS .OR. FULL) THEN
            IF (NOPLAN .EQ. 0) THEN
*              Membrane:
               BGENE(1,K+1)=SHP(2,I)
               BGENE(2,K+2)=SHP(3,I)
*              Couplage membrane-flexion:
               BGENE(1,K+5) = SHP(2,I)  * EXCEN
               BGENE(2,K+4) = -SHP(3,I) * EXCEN
            ENDIF
*           Flexion:
            BGENE(4,K+5)=SHP(2,I)
            BGENE(5,K+4)=-SHP(3,I)
            BGENE(6,K+4)=-SHP(2,I)
            BGENE(6,K+5)=SHP(3,I)
*           Cisaillement transverse:
            BGENE(7,K+3)=SHP(2,I)
            BGENE(7,K+4) = SHP7(4,I)
            BGENE(7,K+5) = SHP7(5,I)
            BGENE(8,K+3)=SHP(3,I)
            BGENE(8,K+4) = SHP8(4,I)
            BGENE(8,K+5) = SHP8(5,I)
         ENDIF
*
         IF (.NOT. GAUSS) THEN
            IF (NOPLAN.EQ.1 .OR. FULL) THEN
*              Membrane:
               BGENE(1,K+1)=SHP(2,I)
               BGENE(2,K+2)=SHP(3,I)
*              Couplage membrane-flexion:
               BGENE(1,K+5) = SHP(2,I)  * EXCEN
               BGENE(2,K+4) = -SHP(3,I) * EXCEN
            ENDIF
*           Cisaillement de membrane et de couplage membrane-flexion:
            BGENE(3,K+1)=SHP(3,I)
            BGENE(3,K+2)=SHP(2,I)
            BGENE(3,K+4) = -SHP(2,I) * EXCEN
            BGENE(3,K+5) = SHP(3,I)  * EXCEN
         ENDIF
*
  300    CONTINUE
*     END DO
*
*
      IF (GAUSS .AND. FULL) THEN
*        On finit de compl{ter la Matrice "B", pour les contraintes ou
*        les d{formations, par des valeurs calcul{es au centre:
*
         DO 500 I=1,NBNO
            SHP(1,I)=SHPTOT(1,I, 5  )
            SHP(2,I)=SHPTOT(2,I, 5  )
            SHP(3,I)=SHPTOT(3,I, 5  )
  500       CONTINUE
*        END DO
*        FONCTIONS DE FORME DANS LA G{OM{TRIE R{ELLE:
         CALL JACOBI(XEL,SHP,2,NBNO,DJAC5)
         IF(DJAC5 .LE. 0.D0) THEN
            IARR=1
            RETURN
         END IF
*
         DO 520 I=1,NBNO
            K = 6*(I-1)
            IF (NOPLAN .EQ. 1) THEN
*              Membrane:
               BGENE(1,K+1)=SHP(2,I)
               BGENE(2,K+2)=SHP(3,I)
*              Couplage membrane-flexion:
               BGENE(1,K+5) = SHP(2,I)  * EXCEN
               BGENE(2,K+4) = -SHP(3,I) * EXCEN
            ENDIF
*           Cisaillement de membrane et de couplage membrane-flexion:
            BGENE(3,K+1)=SHP(3,I)
            BGENE(3,K+2)=SHP(2,I)
            BGENE(3,K+4) = -SHP(2,I) * EXCEN
            BGENE(3,K+5) = SHP(3,I)  * EXCEN
  520       CONTINUE
*        END DO
*
      ENDIF
*
      END