C BBAR SOURCE SP204843 25/07/11 21:15:01 12317 SUBROUTINE BBAR(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, 1 MELE,NBNO,LRE,IFOU,NST,XE,DJAC,A,BB,BGENE) C======================================================================= C= B B A R = C= ----------- = C= Fonction : = C= ---------- = C= CALCUL DE LA MATRICE B-BARRE - ELEMENTS INCOMPRESSIBLES 'IC--' = C= - EN 2D : ELEMENTS ICT3, ICT6, ICQ4, ICQ8 = C= - EN 3D : ELEMENTS ICC8, IC20, ICT4, ICY5 = C= ELEMENTS IC10, ICP6, IC15, IC13 -> A voir = 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= MELE (E) Type de l'element fini (cf. NOMTP dans bdata.eso) = 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= XE (E) Coordonnees des noeuds de l'element fini etudie = C= DJAC (S) Jacobien au point de Gauss etudie = C= BGENE (E/S) Matrice de gradients (B) calcule au point de Gauss = C= A,BB (S) Matrices utiles pour la methode BBar = C======================================================================= C LRE = NOMBRE DE COLONNES DE LA MATRICE B C A = COEFFICIENTS DE MODIFICATION C BB(3,LRE) = MATRICE B-BARRE DILATATION C BGENE(NST,LRE) = MATRICE B C======================================================================= IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) PARAMETER (XZero=0.D0, XUn=1.D0, X1s2=0.5D0, & X1s3=0.333333333333333333333333333333333333333333D0) DIMENSION XE(3,*),BGENE(NST,*), & POIGAU(*),QSIGAU(*),ETAGAU(*),DZEGAU(*), & A(4,*),BB(3,*) * WRITE(6,*) 'Entree dans le sousprogramme BBAR' * WRITE(6,*) 'Element',MELE,NST,IFOU * WRITE(6,*) 'Point Gauss: xsi=',QSIGAU(IGAU),'eta=',ETAGAU(IGAU), * & 'dze=',DZEGAU(IGAU) C----------------------------------------------------------------------- C---- Elements incompressibles (MFR = 31) ------------------------------ C----------------------------------------------------------------------- C NOM : ICT3, ICQ4, ICT6, ICQ8, ICC8, ICT4, ICP6, IC20, IC10, IC15 C MELE : 69 , 70 , 71 , 72 , 73 , 74 , 75 , 76 , 77 , 78 C NOM : ICY5, IC13 C MELE : 273 , 274 , C----------------------------------------------------------------------- C= 0 = Quelques initialisations ======================================== IFR = IFOU+4 C= Indices des composantes diagonales du gradient / deformation C= Attention : NST = 9 pour gradient et NST = 6/4/3 pour deformations I1BG = 1 IF (NST.EQ.9) THEN I2BG = 5 I3BG = 9 ELSE I2BG = 2 I3BG = 3 ENDIF C= 1 = Branchement selon l'Element Fini BBAR ============================ IF (MELE.EQ.273) GOTO 273 IF (MELE.EQ.274) GOTO 274 GOTO ( 69, 70, 71, 72, 73, 74, 75, 76, 77, 78 ), (MELE-68) GOTO 666 C----------------------------------------------------------------------- C----- Element massif bidimensionnel ICT3 ------------------------------ C----------------------------------------------------------------------- 69 CONTINUE *W WRITE(*,*) 'Element ICT3',IFOU GOTO (691,691,691,693),IFR GOTO 666 *----- Contraintes planes ou deformations planes (IFOU=-2,-1) 691 CONTINUE K = 1 DO i = 1,NBNO BGENE(I2BG,K) = X1s2 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I1BG,K) = X1s2 * (BB(1,K)+BGENE(I1BG,K)) K = K+1 BGENE(I1BG,K) = X1s2 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I2BG,K) = X1s2 * (BB(2,K)+BGENE(I2BG,K)) K = K+1 ENDDO GOTO 999 *----- Axisymetrie (IFOU=0) 693 CONTINUE K = 1 DO i = 1,NBNO r_z = X1s3 * (BB(1,K)-BGENE(I1BG,K)+BB(3,K)-BGENE(I3BG,K)) BGENE(I1BG,K) = BGENE(I1BG,K) + r_z BGENE(I2BG,K) = r_z BGENE(I3BG,K) = BGENE(I3BG,K) + r_z K = K + 1 r_z = X1s3 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I1BG,K) = r_z BGENE(I2BG,K) = BGENE(I2BG,K) + r_z BGENE(I3BG,K) = r_z K = K + 1 ENDDO GOTO 999 C----------------------------------------------------------------------- C----- Element massif bidimensionnel ICQ4 ------------------------------ C----------------------------------------------------------------------- 70 CONTINUE *W WRITE(*,*) 'Element ICQ4',IFOU GOTO (700,701,701,703),IFR GOTO 666 *----- Deformations planes generalisees (IFOU=-3) 700 CONTINUE K = 1 DO i = 1,NBNO BGENE(I2BG,K) = X1s2 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I1BG,K) = BGENE(I1BG,K) + BGENE(I2BG,K) BGENE(I3BG,K) = BGENE(I2BG,K) K = K + 1 BGENE(I1BG,K) = X1s2 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I2BG,K) = BGENE(I2BG,K) + BGENE(I1BG,K) BGENE(I3BG,K) = BGENE(I1BG,K) K = K + 1 ENDDO BGENE(I3BG,K+1) = BB(I3BG,K+1) BGENE(I3BG,K+2) = BB(I3BG,K+2) GOTO 999 *----- Contraintes planes ou deformations planes (IFOU=-2,-1) 701 CONTINUE K = 1 DO i = 1,NBNO BGENE(I2BG,K) = X1s2 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I1BG,K) = X1s2 * (BB(1,K)+BGENE(I1BG,K)) K = K + 1 BGENE(I1BG,K) = X1s2 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I2BG,K) = X1s2 * (BB(2,K)+BGENE(I2BG,K)) K = K + 1 ENDDO GOTO 999 *----- Axisymetrie (IFOU=0) 703 CONTINUE K = 1 DO i = 1,NBNO r_z = X1s3 * (BB(1,K)-BGENE(I1BG,K)+BB(3,K)-BGENE(I3BG,K)) BGENE(I1BG,K) = BGENE(I1BG,K) + r_z BGENE(I2BG,K) = r_z BGENE(I3BG,K) = BGENE(I3BG,K) + r_z K = K + 1 r_z = X1s3 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I1BG,K) = r_z BGENE(I2BG,K) = BGENE(I2BG,K) + r_z BGENE(I3BG,K) = r_z K = K + 1 ENDDO GOTO 999 C----------------------------------------------------------------------- C----- Element massif bidimensionnel ICT6 ------------------------------ C----------------------------------------------------------------------- 71 CONTINUE * WRITE(6,*) 'Element ICT6',IFOU GOTO (711,711,711,713),IFR GOTO 666 *----- Contraintes planes ou deformations planes (IFOU=-2,-1) 711 CONTINUE K = 1 DO i = 1,NBNO BGENE(I2BG,K) = X1s2 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I1BG,K) = X1s2 * (BB(1,K)+BGENE(I1BG,K)) K = K + 1 BGENE(I1BG,K) = X1s2 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I2BG,K) = X1s2 * (BB(2,K)+BGENE(I2BG,K)) K = K + 1 ENDDO GOTO 999 *----- Axisymetrie (IFOU=0) 713 CONTINUE K = 1 DO i = 1,NBNO r_z = X1s3 * (BB(1,K)-BGENE(I1BG,K)+BB(3,K)-BGENE(I3BG,K)) BGENE(I1BG,K) = BGENE(I1BG,K) + r_z BGENE(I2BG,K) = r_z BGENE(I3BG,K) = BGENE(I3BG,K) + r_z K = K + 1 r_z = X1s3 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I1BG,K) = r_z BGENE(I2BG,K) = BGENE(I2BG,K) + r_z BGENE(I3BG,K) = r_z K = K + 1 ENDDO GOTO 999 C----------------------------------------------------------------------- C----- Element massif bidimensionnel ICQ8 ------------------------------ C----------------------------------------------------------------------- 72 CONTINUE *W WRITE(*,*) 'Element ICQ8',IFOU QSI_z = QSIGAU(IGAU) ETA_z = ETAGAU(IGAU) GOTO (721,721,721,723),IFR GOTO 666 *----- Contraintes planes ou deformations planes (IFOU=-2,-1) 721 CONTINUE *----- Calculs des composantes de BB-DILATATION BB(1, 1) = A(1, 1) + A(2, 1)*QSI_z + A(3, 1)*ETA_z BB(2, 2) = A(1, 2) + A(2, 2)*QSI_z + A(3, 2)*ETA_z BB(1, 3) = A(1, 3) + A(2, 3)*QSI_z + A(3, 3)*ETA_z BB(2, 4) = A(1, 4) + A(2, 4)*QSI_z + A(3, 4)*ETA_z BB(1, 5) = A(1, 5) + A(2, 5)*QSI_z + A(3, 5)*ETA_z BB(2, 6) = A(1, 6) + A(2, 6)*QSI_z + A(3, 6)*ETA_z BB(1, 7) = A(1, 7) + A(2, 7)*QSI_z + A(3, 7)*ETA_z BB(2, 8) = A(1, 8) + A(2, 8)*QSI_z + A(3, 8)*ETA_z BB(1, 9) = A(1, 9) + A(2, 9)*QSI_z + A(3, 9)*ETA_z BB(2,10) = A(1,10) + A(2,10)*QSI_z + A(3,10)*ETA_z BB(1,11) = A(1,11) + A(2,11)*QSI_z + A(3,11)*ETA_z BB(2,12) = A(1,12) + A(2,12)*QSI_z + A(3,12)*ETA_z BB(1,13) = A(1,13) + A(2,13)*QSI_z + A(3,13)*ETA_z BB(2,14) = A(1,14) + A(2,14)*QSI_z + A(3,14)*ETA_z BB(1,15) = A(1,15) + A(2,15)*QSI_z + A(3,15)*ETA_z BB(2,16) = A(1,16) + A(2,16)*QSI_z + A(3,16)*ETA_z *W WRITE (*,*) 'Ecriture de la matrice BB[ ]' *W WRITE (*,4) (('BB(',I,',',J,')=',BB(I,J),I=1,3),J=1,16) *W 4 FORMAT (3(A,I1,A,I1,A,F8.4,2X)) *----- Calcul des composantes de B-BARRE K = 1 DO i = 1,NBNO BGENE(I2BG,K) = X1s2 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I1BG,K) = X1s2 * (BB(1,K)+BGENE(I1BG,K)) K = K + 1 BGENE(I1BG,K) = X1s2 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I2BG,K) = X1s2 * (BB(2,K)+BGENE(I2BG,K)) K = K + 1 ENDDO GOTO 999 *----- Axisymetrie (IFOU=0) 723 CONTINUE *----- Calculs des composantes de BB-DILATATION BB(1, 1) = A(1, 1) + A(2, 1)*QSI_z + A(3, 1)*ETA_z BB(2, 2) = A(1, 2) + A(2, 2)*QSI_z + A(3, 2)*ETA_z BB(3, 1) = A(1, 3) + A(2, 3)*QSI_z + A(3, 3)*ETA_z BB(1, 3) = A(1, 4) + A(2, 4)*QSI_z + A(3, 4)*ETA_z BB(2, 4) = A(1, 5) + A(2, 5)*QSI_z + A(3, 5)*ETA_z BB(3, 3) = A(1, 6) + A(2, 6)*QSI_z + A(3, 6)*ETA_z BB(1, 5) = A(1, 7) + A(2, 7)*QSI_z + A(3, 7)*ETA_z BB(2, 6) = A(1, 8) + A(2, 8)*QSI_z + A(3, 8)*ETA_z BB(3, 5) = A(1, 9) + A(2, 9)*QSI_z + A(3, 9)*ETA_z BB(1, 7) = A(1,10) + A(2,10)*QSI_z + A(3,10)*ETA_z BB(2, 8) = A(1,11) + A(2,11)*QSI_z + A(3,11)*ETA_z BB(3, 7) = A(1,12) + A(2,12)*QSI_z + A(3,12)*ETA_z BB(1, 9) = A(1,13) + A(2,13)*QSI_z + A(3,13)*ETA_z BB(2,10) = A(1,14) + A(2,14)*QSI_z + A(3,14)*ETA_z BB(3, 9) = A(1,15) + A(2,15)*QSI_z + A(3,15)*ETA_z BB(1,11) = A(1,16) + A(2,16)*QSI_z + A(3,16)*ETA_z BB(2,12) = A(1,17) + A(2,17)*QSI_z + A(3,17)*ETA_z BB(3,11) = A(1,18) + A(2,18)*QSI_z + A(3,18)*ETA_z BB(1,13) = A(1,19) + A(2,19)*QSI_z + A(3,19)*ETA_z BB(2,14) = A(1,20) + A(2,20)*QSI_z + A(3,20)*ETA_z BB(3,13) = A(1,21) + A(2,21)*QSI_z + A(3,21)*ETA_z BB(1,15) = A(1,22) + A(2,22)*QSI_z + A(3,22)*ETA_z BB(2,16) = A(1,23) + A(2,23)*QSI_z + A(3,23)*ETA_z BB(3,15) = A(1,24) + A(2,24)*QSI_z + A(3,24)*ETA_z *W WRITE (*,*) 'Ecriture de la matrice BB[ ]' *W WRITE (*,4) (('BB(',I,',',J,')=',BB(I,J),I=1,3),J=1,16) *----- Calcul des composantes de B-BARRE K = 1 DO i = 1,NBNO r_z = X1s3 * (BB(1,K)-BGENE(I1BG,K)+BB(3,K)-BGENE(I3BG,K)) BGENE(I1BG,K) = BGENE(I1BG,K) + r_z BGENE(I2BG,K) = r_z BGENE(I3BG,K) = BGENE(I3BG,K) + r_z K = K + 1 r_z = X1s3 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I1BG,K) = r_z BGENE(I2BG,K) = BGENE(I2BG,K) + r_z BGENE(I3BG,K) = r_z K = K + 1 ENDDO GOTO 999 C----------------------------------------------------------------------- C----- Element massif tridimensionnel ICC8,ICT4, ICP6 ----------------------------- C----------------------------------------------------------------------- 73 CONTINUE *W WRITE(*,*) 'Element ICC8',IFOU 74 CONTINUE *W WRITE(*,*) 'Element ICT4',IFOU 75 CONTINUE *W WRITE(*,*) 'Element ICP6',IFOU 273 CONTINUE *W WRITE(*,*) 'Element ICY5',IFOU IF (IFOU.NE.2) GOTO 666 * QSI_z = QSIGAU(IGAU) * ETA_z = ETAGAU(IGAU) * DZE_z = DZEGAU(IGAU) K = 1 DO i = 1,NBNO BGENE(I3BG,K) = X1s3 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I2BG,K) = X1s3 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I1BG,K) = X1s3 * (BB(1,K) + 2.D0*BGENE(I1BG,K)) K = K+1 BGENE(I1BG,K) = X1s3 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I3BG,K) = X1s3 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I2BG,K) = X1s3 * (BB(2,K) + 2.D0*BGENE(I2BG,K)) K = K+1 BGENE(I1BG,K) = X1s3 * (BB(3,K)-BGENE(I3BG,K)) BGENE(I2BG,K) = X1s3 * (BB(3,K)-BGENE(I3BG,K)) BGENE(I3BG,K) = X1s3 * (BB(3,K) + 2.D0*BGENE(I3BG,K)) K = K+1 ENDDO GOTO 999 C----------------------------------------------------------------------- C----- Element massif tridimensionnel IC20 IC10 IC15 ------------------- C----------------------------------------------------------------------- 76 CONTINUE *W WRITE(*,*) 'Element IC20',IFOU 77 CONTINUE *W WRITE(*,*) 'Element IC10',IFOU 78 CONTINUE *W WRITE(*,*) 'Element IC15',IFOU 274 CONTINUE *W WRITE(*,*) 'Element IC13',IFOU IF (IFOU.NE.2) GOTO 666 QSI_z = QSIGAU(IGAU) ETA_z = ETAGAU(IGAU) DZE_z = DZEGAU(IGAU) *----- Calculs des composantes de BB-DILATATION K = 0 DO i = 1,NBNO BB(1,K+1) = A(1,K+1) + A(2,K+1)*QSI_z + A(3,K+1)*ETA_z + &A(4,K+1)*DZE_z BB(2,K+2) = A(1,K+2) + A(2,K+2)*QSI_z + A(3,K+2)*ETA_z + &A(4,K+2)*DZE_z BB(3,K+3) = A(1,K+3) + A(2,K+3)*QSI_z + A(3,K+3)*ETA_z + &A(4,K+3)*DZE_z K = K+3 ENDDO *----- Calcul des composantes de B-BARRE K = 1 DO i = 1,NBNO BGENE(I3BG,K) = X1s3 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I2BG,K) = X1s3 * (BB(1,K)-BGENE(I1BG,K)) BGENE(I1BG,K) = X1s3 * (BB(1,K) + 2.D0*BGENE(I1BG,K)) K = K+1 BGENE(I1BG,K) = X1s3 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I3BG,K) = X1s3 * (BB(2,K)-BGENE(I2BG,K)) BGENE(I2BG,K) = X1s3 * (BB(2,K) + 2.D0*BGENE(I2BG,K)) K = K+1 BGENE(I1BG,K) = X1s3 * (BB(3,K)-BGENE(I3BG,K)) BGENE(I2BG,K) = X1s3 * (BB(3,K)-BGENE(I3BG,K)) BGENE(I3BG,K) = X1s3 * (BB(3,K) + 2.D0*BGENE(I3BG,K)) K = K+1 ENDDO GOTO 999 C----------------------------------------------------------------------- C----- ERREUR : Donnees incompatibles ---------------------------------- C----------------------------------------------------------------------- 666 CONTINUE * WRITE(6,*) 'Mode de calcul ',IFOU,' incompatible avec ', * & 'element ',MELE CALL ERREUR(5) C----------------------------------------------------------------------- C----- FIN du sousprogramme -------------------------------------------- C----------------------------------------------------------------------- 999 CONTINUE * WRITE(6,*) 'Sortie du sousprogramme BBAR' RETURN END