C BBCALC SOURCE SP204843 25/07/11 21:15:02 12317 SUBROUTINE BBCALC(MELE,NBNN,MFR,IDIM,XE, 1 NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, 2 NSTRS,LRE,IFOUR,NHRM,A,BB, 3 SHPTOT,SHP,BGENE,XDPGE,YDPGE,P,NOER) ******************************************************************* * projection de dN/dx dans l espace de dimension inferieure c======================================================================= c calcul des composantes de la matrice b-barre-dlatation c - en 2d: elements ict3,icq4,ict6 c - en 3d: elements .... c calcul des coefficients de modification de la matrice b-barre-dilatation c - en 2d: elements icq8 c - en 3d: elements .... c coefficients {a}={a1,a2,a3} c======================================================================= c input c mele = numero de l'element dans nomtp c nbnn = nombre de noeuds c mfr = formulation de l'element c idim = dimension espace c xe = coordonnees de l'element c nbpgau = nombre de points de gauss pour la rigidite c poigau = poids des points de gauss c qsigau,etagau,dzegau = coordonnees des points de gauss c nstrs = nombre de composantes de contraintes c lre = nombre de colonnes de la matrice b c ifour = ifour dans ccoptio c nhrm = numero du mode de fourier c xdpge,ydpge = coordonnees du point autour duquel se fait c le mouvement de la section c c output c tableau des coefficients de modification: a(4,*) c tableau des composantes bbarre-dilatation: bb(3,*) c======================================================================= IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC CCREEL PARAMETER (XZer=0.D0, XUn=1.D0, X1s2=0.5D0, X1s4 = 0.25D0, & X2s3 = 0.666666666666666666666666666666666666666667D0, & X1s3 = 0.333333333333333333333333333333333333333333D0, & X1s6 = 0.166666666666666666666666666666666666666667D0) DIMENSION XE(3,*), POIGAU(*),QSIGAU(*),ETAGAU(*),DZEGAU(*), & SHP(6,*),SHPTOT(6,NBNN,*), & A(4,*),BB(3,*),BGENE(NSTRS,*),P(4,*) DIMENSION XNUM1(8),XNUM2(8),XNUM3(8) * write(*,*) 'Entree dans le sousprogramme BBCALC' * write(*,*) 'Element :',MELE * write(*,*) 'Nombre de points Gauss : ', nbpgau * Code erreur et pas CALL ERREUR pour rester en fortran 77 NOER = 0 IF (MFR.NE.31) GOTO 666 C======================================================================= C= Elements incompressibles implementes : 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 = IFOUR+4 IF (IFR.LE.4 .OR. IFR.GE.7) THEN KINC = 2 ELSE KINC = 3 ENDIF CALL ZERO(BB,3,60) IF (MELE.EQ.72 .OR. MELE.GE.76) THEN CALL ZERO(A,4,60) CALL ZERO(P,4,4) ENDIF C= 1 = Branchement selon l'Element Fini BBAR ============================ IF (MELE.EQ.273 .OR. MELE.EQ.274) GOTO 100 GOTO ( 10, 20, 30, 40,100,100,100,100,100,100), (MELE-68) GOTO 999 3 FORMAT (4(A,I1,A,I1,A,F10.4,2X)) 4 FORMAT (3(A,I1,A,I1,A,F10.4,2X)) C======================================================================= C========== Elements MASSIFS INCOMPRESSIBLES BIDIMENSIONNELS =========== C======================================================================= *----------------------------------------------------------------------- *----------------------------- Element ICT3 ---------------------------- *----------------------------------------------------------------------- 10 CONTINUE * write(6,*) 'Element ICT3',ifour * write(6,*) 'Ecriture de la matrice xe[ ]' * write(6,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,3) GOTO (101,101,101,103),IFR GOTO 999 *--------contraintes planes ou deformations planes (ifour=-2,-1) 101 CONTINUE *--------donnee des composantes de bb-dil r_z = - XE(1,1)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,2)*XE(2,3) & + XE(1,1)*XE(2,2) + XE(1,3)*XE(2,1) - XE(1,3)*XE(2,2) * AIRE = X1s2 * r_z r_z = XUn / r_z BB(1,1) = (XE(2,2)-XE(2,3)) * r_z BB(2,2) = (XE(1,3)-XE(1,2)) * r_z BB(1,3) = (XE(2,3)-XE(2,1)) * r_z BB(2,4) = (XE(1,1)-XE(1,3)) * r_z BB(1,5) = (XE(2,1)-XE(2,2)) * r_z BB(2,6) = (XE(1,2)-XE(1,1)) * r_z * write(6,*) 'calcul de aire pour bbdil' * write(6,*) 'aire = ',AIRE * write(6,*) 'ecriture de la matrice bb' * write(6,3) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn) GOTO 666 *--------axisymetrique (ifour=0) 103 CONTINUE *--------donnee des composantes de bb-dil r_z = - XE(1,1)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,2)*XE(2,3) & + XE(1,1)*XE(2,2) + XE(1,3)*XE(2,1) - XE(1,3)*XE(2,2) * AIRE = X1s6 * r_z * (XE(1,1)+XE(1,2)+XE(1,3)) r_z = XUn / r_z BB(1,1) = r_z * (XE(2,2)-XE(2,3)) BB(2,2) = r_z * (XE(1,3)-XE(1,2)) BB(1,3) = r_z * (XE(2,3)-XE(2,1)) BB(2,4) = r_z * (XE(1,1)-XE(1,3)) BB(1,5) = r_z * (XE(2,1)-XE(2,2)) BB(2,6) = r_z * (XE(1,2)-XE(1,1)) *-- Integration analytique des termes BB(3,1),BB(3,3),BB(3,5) BB(3,1) = XUn / (XE(1,1)+XE(1,2)+XE(1,3)) BB(3,3) = BB(3,1) BB(3,5) = BB(3,1) * write(6,*) 'calcul de aire pour bbdil' * write(6,*) 'aire = ',AIRE * write(6,*) 'ecriture de la matrice bb' * write(6,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn) GOTO 666 *----------------------------------------------------------------------- *----------------------------- Element ICQ4 ---------------------------- *----------------------------------------------------------------------- 20 CONTINUE * write(6,*) 'Element ICQ4',ifour * write(6,*) 'Ecriture de la matrice xe[ ]' * write(6,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,4) GOTO (201,201,201,203),IFR GOTO 999 *--------contraintes planes ou deformations planes (ifour=-2,-1) 201 CONTINUE *--------donnee des composantes de bb-dil r_z = - XE(1,4)*XE(2,3) - XE(1,2)*XE(2,1) + XE(1,1)*XE(2,2) & - XE(1,1)*XE(2,4) - XE(1,3)*XE(2,2) + XE(1,2)*XE(2,3) & + XE(1,3)*XE(2,4) + XE(1,4)*XE(2,1) * AIRE = X1s2 * r_z r_z = XUn / r_z BB(1,1) = r_z * (XE(2,2)-XE(2,4)) BB(2,2) = r_z * (XE(1,4)-XE(1,2)) BB(1,3) = r_z * (XE(2,3)-XE(2,1)) BB(2,4) = r_z * (XE(1,1)-XE(1,3)) BB(1,5) = r_z * (XE(2,4)-XE(2,2)) BB(2,6) = r_z * (XE(1,2)-XE(1,4)) BB(1,7) = r_z * (XE(2,1)-XE(2,3)) BB(2,8) = r_z * (XE(1,3)-XE(1,1)) IF (IFR.EQ.1) THEN XGXE = X1s4 * (XE(1,1) + XE(1,2) + XE(1,3) + XE(1,4)) YGXE = X1s4 * (XE(2,1) + XE(2,2) + XE(2,3) + XE(2,4)) BB(3,9) = Xun BB(3,10) = XDPGE - XGXE BB(3,11) = YGXE - YDPGE ENDIF * write (*,*) 'calcul de aire pour bbdil' * write (*,*) 'aire = ',aire * write (*,*) 'ecriture de la matrice bb' * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn) GOTO 666 * *--------axisymetriques(ifour=0) 203 CONTINUE r_z = XE(1,1)*( XE(1,1)*XE(2,2)-XE(1,1)*XE(2,4)-XE(1,4)*XE(2,4) & +XE(1,4)*XE(2,1)-XE(1,2)*XE(2,1)+XE(1,2)*XE(2,2)) & + XE(1,2)*(-XE(1,2)*XE(2,1)+XE(1,2)*XE(2,3)+XE(1,3)*XE(2,3) & -XE(1,3)*XE(2,2)) & + XE(1,3)*( XE(1,3)*XE(2,4)+XE(1,4)*XE(2,4)-XE(1,4)*XE(2,3) & -XE(1,3)*XE(2,2)) & + XE(1,4)*XE(1,4)*(XE(2,1)-XE(2,3)) * AIRE = x1s6 * r_z r_z = xUn / r_z BB(1,1) = ( XE(2,2)*(XE(1,2)+XE(1,1))-XE(2,4)*(XE(1,1)+XE(1,4)) & +X1s2*( XE(1,3)*(XE(2,2)-XE(2,4)) & +XE(1,4)*(XE(2,2)+XE(2,3)) & -XE(1,2)*(XE(2,3)+XE(2,4))) ) * r_z BB(2,2) = ( XE(1,4)*(XE(1,1)+XE(1,4)) & -XE(1,2)*(XE(1,1)+XE(1,2)) ) * r_z BB(1,3) = ( XE(2,3)*(XE(1,2)+XE(1,3))-XE(2,1)*(XE(1,2)+XE(1,1)) & +X1s2*( XE(1,4)*(XE(2,3)-XE(2,1)) & +XE(1,1)*(XE(2,3)+XE(2,4)) & -XE(1,3)*(XE(2,4)+XE(2,1))) ) * r_z BB(2,4) = ( XE(1,1)*(XE(1,2)+XE(1,1)) & -XE(1,3)*(XE(1,3)+XE(1,2)) ) * r_z BB(1,5) = ( XE(2,4)*(XE(1,3)+XE(1,4))-XE(2,2)*(XE(1,3)+XE(1,2)) & +X1s2*( XE(1,1)*(XE(2,4)-XE(2,2)) & +XE(1,2)*(XE(2,4)+XE(2,1)) & -XE(1,4)*(XE(2,1)+XE(2,2))) ) * r_z BB(2,6) = ( XE(1,2)*(XE(1,2)+XE(1,3)) & -XE(1,4)*(XE(1,4)+XE(1,3)) ) * r_z BB(1,7) = ( XE(2,1)*(XE(1,4)+XE(1,1))-XE(2,3)*(XE(1,4)+XE(1,3)) & +X1s2*( XE(1,2)*(XE(2,1)-XE(2,3)) & +XE(1,3)*(XE(2,1)+XE(2,2)) & -XE(1,1)*(XE(2,2)+XE(2,3))) ) * r_z BB(2,8) = ( XE(1,3)*(XE(1,3)+XE(1,4)) & -XE(1,1)*(XE(1,1)+XE(1,4)) ) * r_z *- Integration analytique des termes BB(3,1),BB(3,3),BB(3,5),BB(3,7) BB(3,1) = ( XE(1,1)*(XE(2,2)-XE(2,4))+XE(2,1)*(XE(1,4)-XE(1,2)) & +X1s2*( XE(1,2)*(XE(2,3)+XE(2,4)) & +XE(1,3)*(XE(2,4)-XE(2,2)) & -XE(1,4)*(XE(2,3)+XE(2,2))) ) * r_z BB(3,3) = ( XE(1,2)*(XE(2,3)-XE(2,1))+XE(2,2)*(XE(1,1)-XE(1,3)) & +X1s2*( XE(1,3)*(XE(2,1)+XE(2,4)) & +XE(1,4)*(XE(2,1)-XE(2,3)) & -XE(1,1)*(XE(2,3)+XE(2,4))) ) * r_z BB(3,5) = ( XE(1,3)*(XE(2,4)-XE(2,2))+XE(2,3)*(XE(1,2)-XE(1,4)) & +X1s2*( XE(1,4)*(XE(2,1)+XE(2,2)) & +XE(1,1)*(XE(2,2)-XE(2,4)) & -XE(1,2)*(XE(2,4)+XE(2,1))) ) * r_z BB(3,7) = ( XE(1,4)*(XE(2,1)-XE(2,3))+XE(2,4)*(XE(1,3)-XE(1,1)) & +X1s2*( XE(1,1)*(XE(2,2)+XE(2,3)) & +XE(1,2)*(XE(2,3)-XE(2,1)) & -XE(1,3)*(XE(2,1)+XE(2,2))) ) * r_z * write (*,*) 'calcul de aire pour bbdil' * write (*,*) 'aire = ',aire * write (*,*) 'ecriture de la matrice bb' * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn) GOTO 666 *----------------------------------------------------------------------- *----------------------------- Element ICT6 ---------------------------- *----------------------------------------------------------------------- 30 CONTINUE * write(*,*) 'Element ICT6',ifour * write (*,*) 'ecriture de la matrice xe[ ]' * write (*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,6) GOTO (301,301,301,303),IFR GOTO 999 * *--------contraintes planes ou deformations planes(ifour=-2,-1) 301 CONTINUE *--------donnee des composantes de bb-dil r_z = ( XE(1,1)*(XE(2,2)-XE(2,6)) + XE(1,2)*(XE(2,3)-XE(2,1)) & + XE(1,3)*(XE(2,4)-XE(2,2)) + XE(1,4)*(XE(2,5)-XE(2,3)) & + XE(1,5)*(XE(2,6)-XE(2,4)) + XE(1,6)*(XE(2,1)-XE(2,5)) ) & + X1s4 * ( XE(1,1)*(XE(2,5)-XE(2,3)) & + XE(1,3)*(XE(2,1)-XE(2,5)) & + XE(1,5)*(XE(2,3)-XE(2,1)) ) AIRE = X2s3 * r_z r_z = xUn / r_z BB(1, 1) = ( (XE(2,2)-XE(2,6)) + X1s4 * (XE(2,5)-XE(2,3)) ) * r_z BB(2, 2) = ( (XE(1,6)-XE(1,2)) + X1s4 * (XE(1,3)-XE(1,5)) ) * r_z BB(1, 3) = (XE(2,3)-XE(2,1)) * r_z BB(2, 4) = (XE(1,1)-XE(1,3)) * r_z BB(1, 5) = ( (XE(2,4)-XE(2,2)) + X1s4 * (XE(2,1)-XE(2,5)) ) * r_z BB(2, 6) = ( (XE(1,2)-XE(1,4)) + X1s4 * (XE(1,5)-XE(1,1)) ) * r_z BB(1, 7) = (XE(2,5)-XE(2,3)) * r_z BB(2, 8) = (XE(1,3)-XE(1,5)) * r_z BB(1, 9) = ( (XE(2,6)-XE(2,4)) + X1s4 * (XE(2,3)-XE(2,1)) ) * r_z BB(2,10) = ( (XE(1,4)-XE(1,6)) + X1s4 * (XE(1,1)-XE(1,3)) ) * r_z BB(1,11) = (XE(2,1)-XE(2,5)) * r_z BB(2,12) = (XE(1,5)-XE(1,1)) * r_z * write (*,*) 'calcul de aire pour bbdil' * write (*,*) 'aire = ',aire * write (*,*) 'ecriture de la matrice bb' * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn) GOTO 666 * *--------axisymetriques (ifour=0) 303 CONTINUE *--------donnee des composantes de bb-dil *-----integration numerique des composantes bb-dil------ AIRE = XZer DO i = 1,nbnn XNUM1(i) = XZer XNUM2(i) = XZer XNUM3(i) = XZer ENDDO DO IGAU = 1,NBPGAU DO i = 1,nbnn SHP(1,i) = SHPTOT(1,i,IGAU) SHP(2,i) = SHPTOT(2,i,IGAU) SHP(3,i) = SHPTOT(3,i,IGAU) ENDDO CALL DISTRR(XE,SHP,nbnn,RR) CALL DISTR2(XE,SHP,nbnn,RR2) CALL DISTR3(XE,SHP,nbnn,RR3) CALL DISTZ2(XE,SHP,nbnn,ZZ2) CALL DISTZ3(XE,SHP,nbnn,ZZ3) CALL JACOBI(XE,SHP,2,nbnn,DJAC) r_z = POIGAU(IGAU) * RR r_zz2 = r_z * ZZ2 r_zz3 = r_z * ZZ3 r_rr2 = r_z * RR2 r_rr3 = r_z * RR3 r_z = POIGAU(IGAU) * DJAC AIRE = AIRE + (r_z * RR) DO i = 1,nbnn XNUM1(i) = XNUM1(i) + ( SHPTOT(2,i,IGAU)*r_zz3 & - SHPTOT(3,i,IGAU)*r_zz2 ) XNUM2(i) = XNUM2(i) + ( SHPTOT(3,i,IGAU)*r_rr2 & - SHPTOT(2,i,IGAU)*r_rr3 ) XNUM3(i) = XNUM3(i) + ( SHPTOT(1,i,IGAU)*r_z ) ENDDO ENDDO r_z = XUn / AIRE DO i = 1,nbnn k = 2*i BB(1,k-1) = XNUM1(i) * r_z BB(2,k) = XNUM2(i) * r_z BB(3,k-1) = XNUM3(i) * r_z ENDDO * write (*,*) 'calcul de aire pour bbdil' * write (*,*) 'aire = ',aire * write (*,*) 'ecriture de la matrice bb' * write (*,4) (('bb(',i,',',j,')=',bb(i,j),i=1,3),j=1,2*nbnn) GOTO 666 *----------------------------------------------------------------------- *----------------------------- Element ICQ8 ---------------------------- *----------------------------------------------------------------------- 40 CONTINUE * write(*,*) 'Element ICQ8',ifour * write(*,*) 'Ecriture de la matrice xe[ ]' * write(*,3) (('xe(',i,',',j,')=',xe(i,j),i=1,2),j=1,6) GOTO (401,401,401,403),IFR GOTO 999 *--------contraintes planes ou deformations planes (ifour=-2,-1) 401 CONTINUE *--------donnee de la matrice symetrique p(3,3) (1 fois par element) P(1,1) = X2s3 * ( XE(1,1)*(XE(2,2)-XE(2,8)) 1 + XE(1,2)*(XE(2,3)-XE(2,1)) 1 +XE(1,3)*(XE(2,4)-XE(2,2))+XE(1,4)*(XE(2,5)-XE(2,3)) 2 +XE(1,5)*(XE(2,6)-XE(2,4))+XE(1,6)*(XE(2,7)-XE(2,5)) 3 +XE(1,7)*(XE(2,8)-XE(2,6))+XE(1,8)*(XE(2,1)-XE(2,7)) ) 4 + X1s6 * ( (XE(1,1)-XE(1,5))*(XE(2,7)-XE(2,3)) 5 +(XE(1,3)-XE(1,7))*(XE(2,1)-XE(2,5)) ) P(1,2) = 5.D0/9.D0*(-XE(2,2)*(XE(1,1)+XE(1,3)) 1 +XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,6)*(-XE(2,5)-XE(2,7)) 1 +XE(2,6)*(XE(1,5)+XE(1,7))) 2 +4.D0/9.D0*(XE(1,8)*(XE(2,7)+XE(2,2)-XE(2,6)-XE(2,1)) 2 +XE(2,8)*(-XE(1,7)-XE(1,2)+XE(1,6)+XE(1,1)) 2 +XE(1,4)*(XE(2,2)-XE(2,6)+XE(2,5)-XE(2,3)) 2 +XE(2,4)*(-XE(1,2)+XE(1,6)-XE(1,5)+XE(1,3))) 3 +7.D0/45.D0*(XE(1,2)*(XE(2,5)+XE(2,7)) 3 -XE(2,2)*(XE(1,5)+XE(1,7))-XE(1,6)*(XE(2,1)+XE(2,3)) 3 +XE(2,6)*(XE(1,1)+XE(1,3))) 4 +8.D0/15.D0*(XE(1,6)*XE(2,2)-XE(1,2)*XE(2,6)) 5 +7.D0/90.D0*(XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5) 5 +XE(1,5)*XE(2,3)-XE(1,1)*XE(2,7)) 6 +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7) 6 -XE(1,5)*XE(2,1)+XE(1,1)*XE(2,5)) P(1,3) = 5.D0/9.D0*(-XE(1,8)*(XE(2,7)+XE(2,1)) 1 +XE(2,8)*(XE(1,7)+XE(1,1)) 1 +XE(1,4)*(XE(2,3)+XE(2,5)) 1 -XE(2,4)*(XE(1,3)+XE(1,5))) 2 +4.D0/9.D0*((XE(1,8)-XE(1,4))*(XE(2,2)+XE(2,6)) 2 +(-XE(2,8)+XE(2,4))*(XE(1,2)+XE(1,6)) 2 +XE(1,2)*(XE(2,1)-XE(2,3))+XE(2,2)*(-XE(1,1)+XE(1,3)) 2 +XE(1,6)*(XE(2,7)-XE(2,5))+XE(2,6)*(-XE(1,7)+XE(1,5))) 3 +7.D0/45.D0*(-XE(1,8)*(XE(2,3)+XE(2,5)) 3 +XE(2,8)*(XE(1,3)+XE(1,5)) 3 +XE(1,4)*(XE(2,1)+XE(2,7)) 3 -XE(2,4)*(XE(1,1)+XE(1,7))) 4 +1.D0/30.D0*(-XE(1,7)*XE(2,3)+XE(1,3)*XE(2,7) 4 +XE(1,5)*XE(2,1)-XE(1,1)*XE(2,5)) 5 +7.D0/90.D0*(XE(1,7)*XE(2,5)-XE(1,3)*XE(2,1) 5 -XE(1,5)*XE(2,7)+XE(1,1)*XE(2,3)) 6 +8.D0/15.D0*(XE(1,8)*XE(2,4)-XE(1,4)*XE(2,8)) P(1,4) = XZer P(2,1) = P(1,2) P(2,2) = 16.D0/45.D0*(XE(2,6)*(XE(1,5)-XE(1,7)) 1 +XE(1,6)*(XE(2,7)-XE(2,5))+XE(1,2)*(XE(2,3)-XE(2,1)) 1 +XE(2,2)*(XE(1,1)-XE(1,3))) 2 +14.D0/45.D0*(XE(2,4)*(XE(1,3)-XE(1,5)) 2 +XE(1,4)*(XE(2,5)-XE(2,3))+XE(1,8)*(XE(2,1)-XE(2,7)) 2 +XE(2,8)*(XE(1,7)-XE(1,1))) 3 +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6)) 3 +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4)) 3 +XE(1,6)*(XE(2,4)-XE(2,8))) 4 +2.D0/45.D0*(XE(2,2)*(XE(1,5)-XE(1,7)) 4 +XE(1,2)*(XE(2,7)-XE(2,5)) +XE(1,6)*(XE(2,3)-XE(2,1)) 4 +XE(2,6)*(XE(1,1)-XE(1,3))) 5 +4.D0/45.D0*(XE(2,4)*(XE(1,1)-XE(1,7)) 5 +XE(1,4)*(XE(2,7)-XE(2,1)) +XE(1,8)*(XE(2,3)-XE(2,5)) 5 +XE(2,8)*(XE(1,5)-XE(1,3))) 6 +17.D0/90.D0*(XE(1,7)*XE(2,5)+XE(1,3)*XE(2,1) 6 -XE(1,5)*XE(2,7)-XE(1,1)*XE(2,3)) 7 +1.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5) 7 +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7)) P(2,3) = 1.D0/3.D0*(XE(1,5)*(XE(2,6)-XE(2,4)) 1 +XE(1,8)*(XE(2,7)+XE(2,1))+XE(1,7)*(XE(2,6)-XE(2,8)) 1 -XE(1,2)*(XE(2,1)+XE(2,3))+XE(1,1)*(XE(2,2)-XE(2,8)) 1 +XE(1,4)*(XE(2,5)+XE(2,3))-XE(1,6)*(XE(2,7)+XE(2,5)) 1 +XE(1,3)*(XE(2,2)-XE(2,4))) 2 +4.D0/9.D0*(-XE(1,4)*(XE(2,2)+XE(2,6)) 2 -XE(1,8)*(XE(2,6)+XE(2,2))+XE(1,2)*(XE(2,8)+XE(2,4)) 2 +XE(1,6)*(XE(2,4)+XE(2,8))) 3 +1.D0/9.D0*(XE(1,7)*(XE(2,2)-XE(2,4)) 3 +XE(1,8)*(XE(2,3)+XE(2,5)) 3 +XE(1,5)*(XE(2,2)-XE(2,8))+XE(1,3)*(XE(2,6)-XE(2,8)) 3 -XE(1,2)*(XE(2,5)+XE(2,7))+XE(1,1)*(XE(2,6)-XE(2,4)) 3 +XE(1,4)*(XE(2,1)+XE(2,7))-XE(1,6)*(XE(2,1)+XE(2,3))) P(2,4) = XZer P(3,1) = P(1,3) P(3,2) = P(2,3) P(3,3) = 14.D0/45.D0*(XE(1,6)*(XE(2,7)-XE(2,5)) 1 +XE(2,6)*(XE(1,5)-XE(1,7))+XE(1,2)*(XE(2,3)-XE(2,1)) 1 +XE(2,2)*(XE(1,1)-XE(1,3))) 2 +16.D0/45.D0*(XE(1,8)*(XE(2,1)-XE(2,7)) 2 +XE(2,8)*(XE(1,7)-XE(1,1))+XE(1,4)*(XE(2,5)-XE(2,3)) 2 +XE(2,4)*(XE(1,3)-XE(1,5))) 3 +8.D0/45.D0*(XE(1,4)*(XE(2,2)-XE(2,6)) 3 +XE(1,8)*(XE(2,6)-XE(2,2))+XE(1,2)*(XE(2,8)-XE(2,4)) 3 +XE(1,6)*(XE(2,4)-XE(2,8))) 4 +2.D0/45.D0*(XE(2,4)*(XE(1,7)-XE(1,1)) 4 +XE(1,8)*(XE(2,5)-XE(2,3))+XE(2,8)*(-XE(1,5)+XE(1,3)) 4 +XE(1,4)*(XE(2,1)-XE(2,7))) 5 +4.D0/45.D0*(XE(2,2)*(XE(1,7)-XE(1,5)) 5 +XE(1,2)*(XE(2,5)-XE(2,7))+XE(1,6)*(XE(2,1)-XE(2,3)) 5 +XE(2,6)*(XE(1,3)-XE(1,1))) 6 +1.D0/90.D0*(-XE(1,5)*XE(2,7)+XE(1,3)*XE(2,1) 6 +XE(1,7)*XE(2,5)-XE(1,1)*XE(2,3)) 7 +17.D0/90.D0*(-XE(1,7)*XE(2,1)-XE(1,3)*XE(2,5) 7 +XE(1,5)*XE(2,3)+XE(1,1)*XE(2,7)) P(3,4) = XZer P(4,1) = XZer P(4,2) = XZer P(4,3) = XZer P(4,4) = XZer * write (*,*) 'ecriture de la matrice p[ ]' * write (*,4) (('p(',i,',',j,')=',p(i,j),i=1,3),j=1,3) *-----------------donnee des vecteurs {t}={t1,t2,t3}-------------------- A(1,1) = 2.D0/3.D0*(XE(2,2)-XE(2,8)) 1 +1.D0/6.D0*(XE(2,7)-XE(2,3)) A(2,1) = 1.D0/9.D0*(4.D0*XE(2,8)-5.D0*XE(2,2)) 1 +7.D0/90.D0*(-XE(2,7)+2.D0*XE(2,6)) 1 +1.D0/30.D0*XE(2,5) A(3,1) = 1.D0/9.D0*(-4.D0*XE(2,2)+5.D0*XE(2,8)) 1 +7.D0/90.D0*(XE(2,3)-2.D0*XE(2,4)) 1 -1.D0/30.D0*XE(2,5) A(1,2) = 2.D0/3.D0*(XE(1,8)-XE(1,2)) 1 +1.D0/6.D0*(XE(1,3)-XE(1,7)) A(2,2) = 1.D0/9.D0*(-4.D0*XE(1,8)+5.D0*XE(1,2)) 1 +7.D0/90.D0*(XE(1,7)-2.D0*XE(1,6)) 1 -1.D0/30.D0*XE(1,5) A(3,2) = 1.D0/9.D0*(4.D0*XE(1,2)-5.D0*XE(1,8)) 1 +7.D0/90.D0*(-XE(1,3)+2.D0*XE(1,4)) 1 +1.D0/30.D0*XE(1,5) A(1,3) = 2.D0/3.D0*(-XE(2,1)+XE(2,3)) A(2,3) = 5.D0/9.D0*(XE(2,1)+XE(2,3)) 1 +7.D0/45.D0*(XE(2,5)+XE(2,7)) 1 -4.D0/9.D0*(XE(2,4)+XE(2,8)) 1 -8.D0/15.D0*XE(2,6) A(3,3) = 4.D0/9.D0*(XE(2,1)-XE(2,3)+XE(2,4)-XE(2,8)) A(1,4) = 2.D0/3.D0*(-XE(1,3)+XE(1,1)) A(2,4) = -5.D0/9.D0*(XE(1,1)+XE(1,3)) 1 -7.D0/45.D0*(XE(1,5)+XE(1,7)) 1 +4.D0/9.D0*(XE(1,4)+XE(1,8)) 1 +8.D0/15.D0*XE(1,6) A(3,4) = 4.D0/9.D0*(-XE(1,1)+XE(1,3)-XE(1,4)+XE(1,8)) A(1,5) = 2.D0/3.D0*(XE(2,4)-XE(2,2)) 1 +1.D0/6.D0*(-XE(2,5)+XE(2,1)) A(2,5) = 1.D0/9.D0*(4.D0*XE(2,4)-5.D0*XE(2,2)) 1 +7.D0/90.D0*(-XE(2,5)+2.D0*XE(2,6)) 1 +1.D0/30.D0*XE(2,7) A(3,5) = 1.D0/9.D0*(4.D0*XE(2,2)-5.D0*XE(2,4)) 1 +7.D0/90.D0*(-XE(2,1)+2.D0*XE(2,8)) 1 +1.D0/30.D0*XE(2,7) A(1,6) = 2.D0/3.D0*(-XE(1,4)+XE(1,2)) 1 +1.D0/6.D0*(XE(1,5)-XE(1,1)) A(2,6) = 1.D0/9.D0*(-4.D0*XE(1,4)+5.D0*XE(1,2)) 1 +7.D0/90.D0*(XE(1,5)-2.D0*XE(1,6)) 1 -1.D0/30.D0*XE(1,7) A(3,6) = 1.D0/9.D0*(-4.D0*XE(1,2)+5.D0*XE(1,4)) 1 +7.D0/90.D0*(XE(1,1)-2.D0*XE(1,8)) 1 -1.D0/30.D0*XE(1,7) A(1,7) = 2.D0/3.D0*(-XE(2,3)+XE(2,5)) A(2,7) = 4.D0/9.D0*(XE(2,5)-XE(2,3)+XE(2,2)-XE(2,6)) A(3,7) = 5.D0/9.D0*(XE(2,3)+XE(2,5)) 1 +7.D0/45.D0*(XE(2,1)+XE(2,7)) 1 -4.D0/9.D0*(XE(2,2)+XE(2,6)) 1 -8.D0/15.D0*XE(2,8) A(1,8) = 2.D0/3.D0*(XE(1,3)-XE(1,5)) A(2,8) = 4.D0/9.D0*(-XE(1,5)+XE(1,3)-XE(1,2)+XE(1,6)) A(3,8) = -5.D0/9.D0*(XE(1,3)+XE(1,5)) 1 -7.D0/45.D0*(XE(1,1)+XE(1,7)) 1 +4.D0/9.D0*(XE(1,2)+XE(1,6)) 1 +8.D0/15.D0*XE(1,8) A(1,9) = 2.D0/3.D0*(-XE(2,4)+XE(2,6)) 1 +1.D0/6.D0*(XE(2,3)-XE(2,7)) A(2,9) = 1.D0/9.D0*(-4.D0*XE(2,4)+5.D0*XE(2,6)) 1 +7.D0/90.D0*(XE(2,3)-2.D0*XE(2,2)) 1 -1.D0/30.D0*XE(2,1) A(3,9) = 1.D0/9.D0*(4.D0*XE(2,6)-5.D0*XE(2,4)) 1 +7.D0/90.D0*(-XE(2,7)+2.D0*XE(2,8)) 1 +1.D0/30.D0*XE(2,1) A(1,10) = 2.D0/3.D0*(XE(1,4)-XE(1,6)) 1 +1.D0/6.D0*(-XE(1,3)+XE(1,7)) A(2,10) = 1.D0/9.D0*(4.D0*XE(1,4)-5.D0*XE(1,6)) 1 +7.D0/90.D0*(-XE(1,3)+2.D0*XE(1,2)) 1 +1.D0/30.D0*XE(1,1) A(3,10) = 1.D0/9.D0*(-4.D0*XE(1,6)+5.D0*XE(1,4)) 1 +7.D0/90.D0*(XE(1,7)-2.D0*XE(1,8)) 1 -1.D0/30.D0*XE(1,1) A(1,11) = 2.D0/3.D0*(-XE(2,5)+XE(2,7)) A(2,11) = -5.D0/9.D0*(XE(2,7)+XE(2,5)) 1 -7.D0/45.D0*(XE(2,1)+XE(2,3)) 1 +4.D0/9.D0*(XE(2,4)+XE(2,8)) 1 +8.D0/15.D0*XE(2,2) A(3,11) = 4.D0/9.D0*(XE(2,4)-XE(2,5)-XE(2,8)+XE(2,7)) A(1,12) = 2.D0/3.D0*(XE(1,5)-XE(1,7)) A(2,12) = 5.D0/9.D0*(XE(1,7)+XE(1,5)) 1 +7.D0/45.D0*(XE(1,1)+XE(1,3)) 1 -4.D0/9.D0*(XE(1,4)+XE(1,8)) 1 -8.D0/15.D0*XE(1,2) A(3,12) = 4.D0/9.D0*(-XE(1,4)+XE(1,5)+XE(1,8)-XE(1,7)) A(1,13) = 2.D0/3.D0*(-XE(2,6)+XE(2,8)) 1 +1.D0/6.D0*(XE(2,5)-XE(2,1)) A(2,13) = 1.D0/9.D0*(-4.D0*XE(2,8)+5.D0*XE(2,6)) 1 +7.D0/90.D0*(XE(2,1)-2.D0*XE(2,2)) 1 -1.D0/30.D0*XE(2,3) A(3,13) = 1.D0/9.D0*(-4.D0*XE(2,6)+5.D0*XE(2,8)) 1 +7.D0/90.D0*(XE(2,5)-2.D0*XE(2,4)) 1 -1.D0/30.D0*XE(2,3) A(1,14) = 2.D0/3.D0*(XE(1,6)-XE(1,8)) 1 +1.D0/6.D0*(-XE(1,5)+XE(1,1)) A(2,14) = 1.D0/9.D0*(4.D0*XE(1,8)-5.D0*XE(1,6)) 1 +7.D0/90.D0*(-XE(1,1)+2.D0*XE(1,2)) 1 +1.D0/30.D0*XE(1,3) A(3,14) = 1.D0/9.D0*(4.D0*XE(1,6)-5.D0*XE(1,8)) 1 +7.D0/90.D0*(-XE(1,5)+2.D0*XE(1,4)) 1 +1.D0/30.D0*XE(1,3) A(1,15) = 2.D0/3.D0*(-XE(2,7)+XE(2,1)) A(2,15) = 4.D0/9.D0*(-XE(2,6)+XE(2,7)+XE(2,2)-XE(2,1)) A(3,15) = -5.D0/9.D0*(XE(2,1)+XE(2,7)) 1 -7.D0/45.D0*(XE(2,3)+XE(2,5)) 1 +4.D0/9.D0*(XE(2,2)+XE(2,6)) 1 +8.D0/15.D0*XE(2,4) A(1,16) = 2.D0/3.D0*(XE(1,7)-XE(1,1)) A(2,16) = 4.D0/9.D0*(XE(1,6)-XE(1,7)-XE(1,2)+XE(1,1)) A(3,16) = 5.D0/9.D0*(XE(1,1)+XE(1,7)) 1 +7.D0/45.D0*(XE(1,3)+XE(1,5)) 1 -4.D0/9.D0*(XE(1,2)+XE(1,6)) 1 -8.D0/15.D0*XE(1,4) * write (*,*) 'ecriture de la matrice t[ ]' * write (*,4) (('t(',i,',',j,')=',a(i,j),i=1,3),j=1,16) *--------------resolution matricielle de [p]{a}={t}--------------------- * write (*,*) 'Appel a GAUSSJ' N = 3 NP = 4 M = 16 MP = 60 CALL GAUSSJ(P,N,NP,A,M,MP) GOTO 666 *--------axisymetriques (ifour=0) 403 CONTINUE *--------donnee de la matrice symetrique p(3,3) (1 fois par element) CALL ZERO(P,4,4) DO IGAU = 1,NBPGAU DO I = 1,nbnn SHP(1,I) = SHPTOT(1,I,IGAU) SHP(2,I) = SHPTOT(2,I,IGAU) SHP(3,I) = SHPTOT(3,I,IGAU) ENDDO CALL DISTRR(XE,SHP,nbnn,RR) CALL JACOBI(XE,SHP,2,nbnn,DJAC) r_z = POIGAU(IGAU)*RR*DJAC P(1,1) = P(1,1) + r_z P(2,1) = P(2,1) + (r_z * QSIGAU(IGAU)) P(3,1) = P(3,1) + (r_z * ETAGAU(IGAU)) P(2,2) = P(2,2) + (r_z * QSIGAU(IGAU) * QSIGAU(IGAU)) P(3,2) = P(3,2) + (r_z * QSIGAU(IGAU) * ETAGAU(IGAU)) P(3,3) = P(3,3) + (r_z * ETAGAU(IGAU) * ETAGAU(IGAU)) ENDDO P(1,2) = P(2,1) P(1,3) = P(3,1) P(2,3) = P(3,2) * write (*,*) 'ecriture de la matrice p[ ]' * write (*,4) (('p(',i,',',j,')=',p(i,j),i=1,3),j=1,3) *-----------calcul des vecteurs {t}={t1,t2,t3}--integration numerique--- CALL ZERO(A,4,24) DO 4031 IGAU = 1,NBPGAU DO I = 1,nbnn SHP(1,I) = SHPTOT(1,I,IGAU) SHP(2,I) = SHPTOT(2,I,IGAU) SHP(3,I) = SHPTOT(3,I,IGAU) ENDDO CALL DISTRR(XE,SHP,nbnn,RR) CALL DISTR2(XE,SHP,nbnn,RR2) CALL DISTR3(XE,SHP,nbnn,RR3) CALL DISTZ2(XE,SHP,nbnn,ZZ2) CALL DISTZ3(XE,SHP,nbnn,ZZ3) CALL JACOBI(XE,SHP,2,nbnn,DJAC) K = 0 DO I = 1,nbnn K = K+1 r_z = POIGAU(IGAU)*RR & *(SHPTOT(2,I,IGAU)*ZZ3-SHPTOT(3,I,IGAU)*ZZ2) A(1,K) = A(1,K) + r_z A(2,K) = A(2,K) + r_z * QSIGAU(IGAU) A(3,K) = A(3,K) + r_z * ETAGAU(IGAU) K = K+1 r_z = POIGAU(IGAU)*RR & *(-SHPTOT(2,I,IGAU)*RR3+SHPTOT(3,I,IGAU)*RR2) A(1,K) = A(1,K) + r_z A(2,K) = A(2,K) + r_z * QSIGAU(IGAU) A(3,K) = A(3,K) + r_z * ETAGAU(IGAU) K = K+1 r_z = POIGAU(IGAU)*SHPTOT(1,I,IGAU)*DJAC A(1,K) = A(1,K) + r_z A(2,K) = A(2,K) + r_z * QSIGAU(IGAU) A(3,K) = A(3,K) + r_z * ETAGAU(IGAU) ENDDO 4031 CONTINUE * write(*,*) 'Ecriture de la matrice t[ ]' * write(*,4) (('t(',i,',',j,')=',a(i,j),i=1,3),j=1,24) *------------resolution matricielle de [p]{a}={t}----------------- N = 3 NP = 4 M = 24 MP = 60 * write (*,*) 'appel a GAUSSJ' CALL GAUSSJ(P,N,NP,A,M,MP) * write(*,*) 'Ecriture de la matrice a[ ]' * write(*,4) (('a(',i,',',j,')=',a(i,j),i=1,3),j=1,24) GOTO 666 C======================================================================= C========== Elements MASSIFS INCOMPRESSIBLES TRIDIMENSIONNELS ========== C======================================================================= 100 CONTINUE IF (ifour.NE.2) GOTO 999 DIM3 = 1.D0 FACAR = 1.D0 FACSCA = 1.D0 ESTEL = XZero C= Boucle sur les points d'integration DO iGau=1,NBPGAU * BGENE au point IGAU CALL BMATST(IGAU,NBPGAU,POIGAU,QSIGAU,ETAGAU,DZEGAU, & MELE,MFR,NBNN,LRE,IFOUR,NSTRS,NHRM,DIM3, & XE,SHPTOT,SHP,BGENE,DJAC,XDPGE,YDPGE) C= Calcul de la composante integree en ce point de Gauss DJAC=ABS(DJAC)*POIGAU(iGau)*FACAR*DIM3 ESTEL=ESTEL+FACSCA*DJAC *----- Element ICQ8 ICCU20(76) ICTE10(77) ICPR15(78) ICPY13(274) IF (MELE.EQ.72 .or. (MELE.ge.76.and.MELE.le.78) .or. & MELE.EQ.274) then * second membre T K = 1 jj = 0 DO i = 1,NBNN jj = jj + 1 A(1,jj) = A(1,jj) + BGENE(1,K)*DJAC A(2,jj) = A(2,jj) + BGENE(1,K)*qsigau(igau)*DJAC A(3,jj) = A(3,jj) + BGENE(1,K)*etagau(igau)*DJAC if (mele.ge.76) & A(4,jj) = A(4,jj) + BGENE(1,K)*dzegau(igau)*DJAC jj = jj + 1 A(1,jj) = A(1,jj) + BGENE(2,K+1)*DJAC A(2,jj) = A(2,jj) + BGENE(2,K+1)*qsigau(igau)*DJAC A(3,jj) = A(3,jj) + BGENE(2,K+1)*etagau(igau)*DJAC if (mele.ge.76) then A(4,jj) = A(4,jj) + BGENE(2,K+1)*dzegau(igau)*DJAC jj = jj + 1 A(1,jj) = A(1,jj) + BGENE(3,K+2)*DJAC A(2,jj) = A(2,jj) + BGENE(3,K+2)*qsigau(igau)*DJAC A(3,jj) = A(3,jj) + BGENE(3,K+2)*etagau(igau)*DJAC A(4,jj) = A(4,jj) + BGENE(3,K+2)*dzegau(igau)*DJAC endif K = K + KINC ENDDO * matrice M P(1,1) = P(1,1) + DJAC P(1,2) = P(1,2) + qsigau(igau)*DJAC P(1,3) = P(1,3) + etagau(igau)*DJAC P(2,2) = P(2,2) + qsigau(igau)*qsigau(igau)*DJAC P(2,3) = P(2,3) + qsigau(igau)*etagau(igau)*DJAC P(3,3) = P(3,3) + etagau(igau)*etagau(igau)*DJAC if (mele.ge.76) then P(1,4) = P(1,4) + dzegau(igau)*DJAC P(2,4) = P(2,4) + qsigau(igau)*dzegau(igau)*DJAC P(3,4) = P(3,4) + etagau(igau)*dzegau(igau)*DJAC P(4,4) = P(4,4) + dzegau(igau)*dzegau(igau)*DJAC endif else if(mele.le.71) then K = 1 DO i = 1,NBNN bb(1,K) = bb(1,K) + BGENE(1,K)*DJAC bb(2,K+1) = bb(2,K+1) + BGENE(2,K+1)*DJAC bb(3,K) = bb(3,K) + BGENE(3,K)*DJAC K = K + KINC ENDDO else if ( (mele.ge.73.and.mele.le.75) &.or.mele.eq.273) then K = 1 DO i = 1,NBNN bb(1,K) = bb(1,K) + BGENE(1,K)*DJAC bb(2,K+1) = bb(2,K+1) + BGENE(2,K+1)*DJAC bb(3,K+2) = bb(3,K+2) + BGENE(3,K+2)*DJAC K = K + KINC ENDDO endif ENDDO if (MELE.EQ.72) then P(1,4) = XZer P(2,1) = P(1,2) P(2,4) = XZer P(3,1) = P(1,3) P(3,2) = P(2,3) P(3,4) = XZer P(4,1) = XZer P(4,2) = XZer P(4,3) = XZer P(4,4) = XZer ELSE IF (mele.GE.76) THEN P(2,1) = P(1,2) P(3,1) = P(1,3) P(3,2) = P(2,3) P(4,1) = P(1,4) P(4,2) = P(2,4) P(4,3) = P(3,4) ENDIF * ** synthese * if (MELE.EQ.72) then N = 3 NP = 4 M = 16 MP = 60 * write (*,*) 'appel a GAUSSJ' CALL GAUSSJ(P,N,NP,A,M,MP) * write(*,*) 'Ecriture de la matrice a[ ]' * write(*,4) (('a(',i,',',j,')=',a(i,j),i=1,3),j=1,24) elseif (mele.ge.76.and.mele.ne.273) then N = 4 NP = 4 M = NBNN*3 MP = 60 * write (*,*) 'appel a GAUSSJ' CALL GAUSSJ(P,N,NP,A,M,MP) elseif(mele.le.71) then K = 1 DO i = 1,NBNN bb(1,K) = bb(1,K)/ESTEL bb(2,K+1) = bb(2,K+1)/ESTEL bb(3,K) = bb(3,K)/ESTEL K = K + KINC ENDDO elseif((mele.ge.73.and.mele.le.75).or. &mele.eq.273) then K = 1 DO i = 1,NBNN bb(1,K) = bb(1,K)/ESTEL bb(2,K+1) = bb(2,K+1)/ESTEL bb(3,K+2) = bb(3,K+2)/ESTEL K = K + KINC ENDDO endif GOTO 666 C======================================================================= C======================= ERREUR : Cas non prevus ======================= C======================================================================= 999 CONTINUE C write(6,*) 'Erreur d aiguillage - cas non prevu',MELE,IFOUR C CALL ERREUR(922) NOER = 922 666 CONTINUE C write(6,*) 'bbcalc:MELE,IFOUR',MELE,IFOUR RETURN END