C SHB8 SOURCE PV 18/06/18 21:15:29 9860 SUBROUTINE SHB8(ICLE,XE,D,PROPEL,WORK,RE,OUT) IMPLICIT REAL *8(A-H,O-Z) * SUBROUTINE SHB8(ICLE,XE,D,PROPEL,WORK,RE,OUT) C C ELEMENT SHB8 C INTEGER IOP,IBID,IDPLA(5),LAG,IRDC DIMENSION EYG(5),ENU(5) DIMENSION CMATTMP(6,6) DIMENSION XE(24),D(*),PROPEL(*),WORK(30),RE(24,24),OUT(*) DIMENSION XXG5(5),PXG5(5),XCOQ(3,4),BKSIP(3,8,5),B(3,8) DIMENSION XELOC(24),XCENT(3),PPP(3,3),PPPT(3,3) DIMENSION XL(3,4),XXX(3),YYY(3),BGL_LOC(6,24) DIMENSION BGL_LOCT(24,6),TMPTAB(6,24),TMPKE(24,24),CMATLOC(6,6) DIMENSION XELOCTMP(24) DIMENSION XXVB(3),HIJ(6),XKSTAB(24,24),TMPKE2(24,24) DIMENSION GB(8,4),GS(8,4),VB(8,3),XXGB(3,4) DIMENSION XK11(8,8),XK22(8,8),XK33(8,8),RR2(3,3),XK12(8,8) DIMENSION XK21(8,8),XK13(8,8),XK23(8,8),XK31(8,8),XK32(8,8) DIMENSION XR(3,3),XRT(3,3) C DIMENSION DEPS(6),DUSDX(9),UE(3,8),SIG(6),XE_LOCP(24) DIMENSION UDEF(24),XE_LOC(24),XE_LOC12(24) DIMENSION DEPS_LOC(6),SIG_LOC(6) C DIMENSION F(3,8),BSIG(24),SIGMAG(6),ULOC(3,8),PQIALPH(3,4) DIMENSION QIALPHA(3,4),FHG(3,8),RR12(3,3),RR1(3,3),FHG24(24) C DIMENSION DIREC(3),V(3),X56(3),X67(3),X78(3),X85(3),X12(3) DIMENSION X23(3),X34(3),X41(3),V56_78(3),V85_67(3),VSUP(3) DIMENSION V12_34(3),V41_23(3),VINF(3) C DIMENSION XKSIGTMP1(8,8),XKSIGTMP2(8,8) C DIMENSION XCOQ14(3),XCOQ23(3),XCOQ34(3),XCOQ12(3),XNORMALE(3) DIMENSION XCOQ12_34(3),XCOQ14_23(3),XKP(24,24),BKP(3,4) C DIMENSION XMASSE(8,8) C DIMENSION FQ(24) C C AVRIL 2004: PROGRAMMATION DU SHB8 AVEC B_2_CHAPEAU_BAR C DIMENSION CB2BAR(6,6),SIGMOY(6) C CCCCCCCCCCCCCC ENTREES CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ICLE=2 ON CALCULE LA MATRICE DE RAIDEUR C ICLE=3 ON CALCULE LA MATRICE DE MASSE C ICLE=4 ON CALCULE LA MATRICE TANGENTE C ICLE=5 ON CALCULE LES FORCES DE PRESSION C ICLE=7 ON CALCULE LES CONTRAINTES C ICLE=8 ON CALCULE BSIGMA C ICLE=9 ON CALCULE KSIGMA C ICLE=10 ON CALCULE KP C ICLE=11 On calcule epsi a partir de deplacement C icle=12 on calcule epsi à partir des contraintes C ou sigma à partir de epsi suivan propel(3) C C D MATRICE D'ELASTICITE C PROPEL PROPRIETES ELASTIQUES (1)=YUNG (2)=NU C (3)=RO (4)=ALPHA C WORK TABLEAU DE TRAVAIL CCCCCCCCCCCCCC SORTIE CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C RE MATRICE C OUT VECTEUR FORCE OU VECTEUR CONTRAINTE AUX NOEUDS C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C INITIALISATIONS C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C INFOS: C XE EST RANGE COMME CA: (XNOEUD1 YNOEUD1 ZNOEUD1, XNOEUD2 YNOEUD2 ZNOEUD2,...) C DANS SHB8_TEST_NUM: ATTENTION A LA NUMEROTATION DES NOEUDS CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DATA GB/1.D0, 1.D0,-1.D0,-1.D0,-1.D0,-1.D0, 1.D0, 1.D0, & 1.D0,-1.D0,-1.D0, 1.D0,-1.D0, 1.D0, 1.D0,-1.D0, & 1.D0,-1.D0, 1.D0,-1.D0, 1.D0,-1.D0, 1.D0,-1.D0, & -1.D0, 1.D0,-1.D0, 1.D0, 1.D0,-1.D0, 1.D0,-1.D0/ C C VB: COORD DES NOEUDS DANS REPERE DE REFERENCE C DATA VB/-1.D0, 1.D0, 1.D0,-1.D0,-1.D0, 1.D0, 1.D0,-1.D0, & -1.D0,-1.D0, 1.D0, 1.D0,-1.D0,-1.D0, 1.D0, 1.D0, & -1.D0,-1.D0,-1.D0,-1.D0, 1.D0, 1.D0, 1.D0, 1.D0/ C C ON DEFINI LES POINTS GAUSS ET LES POIDS C XXG5(1) = -0.906179845938664D0 XXG5(2) = -0.538469310105683D0 XXG5(3) = 0.D0 XXG5(4) = 0.538469310105683D0 XXG5(5) = 0.906179845938664D0 C PXG5(1) = 0.236926885056189D0 PXG5(2) = 0.478628670499366D0 PXG5(3) = 0.568888888888889D0 PXG5(4) = 0.478628670499366D0 PXG5(5) = 0.236926885056189D0 C UNS8 = 1.D0/8.D0 UNS3 = 1.D0/3.D0 UNS5 = 1.D0/5.D0 UNS6 = 1.D0/6.D0 ZERO = 0.D0 C C TYPE DE LOI DE COMPORTEMENT: C IRDC = 1 : SHB8 TYPE PLEXUS C IRDC = 2 : C.P. C IRDC = 3 : 3D COMPLETE C C IRDC = 1 IRDC=nint(OUT(1)) IF(ICLE.EQ.2)THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ON CALCULE LA RAIDEUR : SORTIE DANS RE C C C C SI IETAN = 1 , ALORS ON CALCULE AUSSI C C LA MATRICE TANGENTE PLASTIQUE C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IETAN = nint(PROPEL(13)) C C INTIALISATION LONGUEUR DES COTES C CALCUL DES COEFF D ELANCEMENT A METTRE DANS LA MATRICE DE CPT C XXL1 = 0.D0 XXL2 = 0.D0 TT1 = 0.D0 TT2 = 0.D0 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C STABILISATION ADAPTATIVE EN FONCTION DE LA DISTORTION DE L'ELEMENT C Cbp POUR L'INSTANT, ON NE MET PAS EN SERVICE => on commente Cbp DO I=1,3 Cbp C DISTANCE ENTRE 1 ET 5 (EPAISSEUR) Cbp TT1 = TT1+(XE(I+12)-XE(I))**2 Cbp C DISTANCE ENTRE 3 ET 7 (EPAISSEUR) Cbp TT2 = TT2+(XE(I+18)-XE(I+6))**2 Cbp C DISTANCE ENTRE 1 ET 2 Cbp XXL1 = XXL1+(XE(I+3)-XE(I))**2 Cbp C DISTANCE ENTRE 2 ET 3 Cbp XXL2 = XXL2+(XE(I+6)-XE(I+3))**2 Cbp ENDDO Cbp XXL1 = SQRT(XXL1) Cbp XXL2 = SQRT(XXL2) Cbp TT1 = 0.5*(SQRT(TT1)+SQRT(TT2)) Cbp COELA1 = 5./6. Cbp COELA2 = 5./6. Cbp C ELANCEMENT DANS DIRECTION 2 Cbp ELT = 6.*TT1/XXL1 Cbp IF(COELA1.GT.ELT) COELA1=ELT Cbp C ELANCEMENT DANS DIRECTION 1 Cbp ELT = 6.*TT1/XXL2 Cbp IF(COELA2.GT.ELT) COELA2=ELT C POUR L'INSTANT, ON NE MET PAS EN SERVICE: COELA1 = 1.D0 COELA2 = 1.D0 C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IF (IETAN.EQ.1) THEN EYG(1) = PROPEL(3) EYG(2) = PROPEL(4) EYG(3) = PROPEL(5) EYG(4) = PROPEL(6) EYG(5) = PROPEL(7) ENU(1) = PROPEL(8) ENU(2) = PROPEL(9) ENU(3) = PROPEL(10) ENU(4) = PROPEL(11) ENU(5) = PROPEL(12) DO I=1,5 C SI IDPLA(I)=1, POINT GAUSS (I) PLASTIQUE IDPLA(I)=nint(RE(I,1)) ENDDO ENDIF C CALL ZDANUL(RE,576) CALL ZDANUL(CMATLOC,36) C CALL ZDANUL(CMATTMP,36) IELEM=nint(WORK(1)) C C PROCEDURE DE TEST DE NUMEROTATION DES ELEMENTS: MARCHE BIEN, MAIS DEJA FAIT EN C AMONT C C CALL SHB8_TEST_NUM(XE,IELEM) C C ON DEFINI CMATLOC: MATRICE DE COMPORTEMENT C IF (IRDC.LE.2) THEN C COMPORTEMENT SHB8 PLEXUS ou C.P. XNU = PROPEL(2) XNU_B = XNU / (1-XNU) XLAMBDA = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2)) XMU = 0.5D0 * PROPEL(1) / (1+PROPEL(2)) CMATLOC(1,1) = XLAMBDA+2*XMU CMATLOC(1,2) = XLAMBDA CMATLOC(2,1) = XLAMBDA CMATLOC(2,2) = XLAMBDA+2*XMU IF (IRDC.EQ.1) THEN C COMPORTEMENT SHB8 PLEXUS CMATLOC(3,3) = PROPEL(1) ENDIF IF (IRDC.EQ.2) THEN C COMPORTEMENT C.P. CMATLOC(3,3) = 0.D0 ENDIF CMATLOC(4,4) = XMU CMATLOC(5,5) = XMU CMATLOC(6,6) = XMU ENDIF C IF (IRDC.EQ.3) THEN C COMPORTEMENT LOI TRIDIM MMC 3D XNU = PROPEL(2) XCOOEF = PROPEL(1)/((1+XNU)*(1-2*XNU)) CMATLOC(1,1) = (1-XNU)*XCOOEF CMATLOC(2,2) = (1-XNU)*XCOOEF CMATLOC(3,3) = (1-XNU)*XCOOEF CMATLOC(1,2) = XNU*XCOOEF CMATLOC(2,1) = XNU*XCOOEF CMATLOC(1,3) = XNU*XCOOEF CMATLOC(3,1) = XNU*XCOOEF CMATLOC(2,3) = XNU*XCOOEF CMATLOC(3,2) = XNU*XCOOEF CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF ENDIF Cbp cas non prevu IF (IRDC.GT.3) CALL ERREUR(5) C C CALL SHIFTD(CMATLOC,CMATTMP,36) C C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE C BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP C BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP C BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP C CALL CUBKSI(5,XXG5,BKSIP) C C DEBUT DE LA BOUCLE SUR LES 5 PTS GAUSS C DO IP=1,5 C C SI POINT PLASTIQUE, ON ITERE AVEC LE MODULE TANGENT C C IF (IETAN.EQ.1) THEN C IF (IDPLA(IP).EQ.1) THEN C CMATLOC(1,1) = CMATTMP(1,1)*EYG(IP)/PROPEL(1) C CMATLOC(2,2) = CMATTMP(2,2)*EYG(IP)/PROPEL(1) C CMATLOC(3,3) = CMATTMP(3,3)*EYG(IP)/PROPEL(1) C CMATLOC(1,2) = CMATTMP(1,2)*EYG(IP)/PROPEL(1) C CMATLOC(2,1) = CMATTMP(2,1)*EYG(IP)/PROPEL(1) C CMATLOC(4,4) = CMATTMP(4,4)*EYG(IP)/PROPEL(1) C CMATLOC(5,5) = CMATTMP(5,5)*EYG(IP)/PROPEL(1) C CMATLOC(6,6) = CMATTMP(6,6)*EYG(IP)/PROPEL(1) C ELSE C CMATLOC(1,1) = CMATTMP(1,1) C CMATLOC(2,2) = CMATTMP(2,2) C CMATLOC(3,3) = CMATTMP(3,3) C CMATLOC(1,2) = CMATTMP(1,2) C CMATLOC(2,1) = CMATTMP(2,1) C CMATLOC(4,4) = CMATTMP(4,4) C CMATLOC(5,5) = CMATTMP(5,5) C CMATLOC(6,6) = CMATTMP(6,6) C ENDIF C ENDIF C C DEFINITION DES 4 POINTS COQUES C ZETA = XXG5(IP) ZLAMB = 0.5*(1.-ZETA) DO I=1,4 DO J=1,3 XCOQ(J,I) = ZLAMB*XE((I-1)*3+J) & + (1.-ZLAMB)*XE((I-1+4)*3+J) END DO END DO C C CALCUL DE PPP 3*3 PASSAGE DE GLOBAL A LOCAL,COQUE C XCENT : COORD GLOBAL DU CENTRE DE L'ELEMENT C CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID) C C CALCUL DE B EN GLOBAL C C ATTENTION A L'ORDRE DE EPSILON: C FARID DANS SON PAPIER: 11 22 33 12 13 23 C HARID DANS PLEXUS: 11 22 33 12 23 13 C C C ON RAJOUTE LES TERMES H1,X . G1 , H2,X . G2 C ET H1,Y . G1 , H2,Y . G2 C AVEC H1 = Y.Z H2 = X.Z C DONC H1,X =0 H2,X = Z C ET H1,Y =Z H2,Y = 0 C C DONC IL NE RESTE PLUS QU'A CALCULER G1 ET G2, ET A AJOUTER A BKSIP C CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC) C C IL FAUT: EPS_LOCAL=BGLOB .U_NODAL_GLOBAL C ON CALCULE BGL_LOC LA MATRICE B(6,24) UGLOBAL ---> EPS LOCAL C CALL ASBGLB(BGL_LOC,B,PPP) C C IL FAUT TRANSPOSER BGL_LOC C CALL AEQBT(BGL_LOCT,BGL_LOC,6,24) C C CALCUL DE LA MATRICE TANGENTE PLASTIQUE: VOIR FICHIER SHB8_MAT_TGE_PLAS.FFF C C C IL NE RESTE PLUS QU'A FAIRE: BGL_LOCT * C * BGL_LOC C CALL ZDANUL(TMPTAB,144) CALL ZDANUL(TMPKE,576) CALL ZDANUL(TMPKE2,576) * CALL MULMAT(6,6,24,CMATLOC,BGL_LOC,TMPTAB) CALL MULMAT( TMPTAB,CMATLOC,BGL_LOC,6,24,6) * CALL MULMAT(24,6,24,BGL_LOCT,TMPTAB,TMPKE2) CALL MULMAT(TMPKE2,BGL_LOCT,TMPTAB,24,24,6) C C ASSEMBLAGE: KE=KE + POIDS*JACOBIAN*TMPKE C DO J=1,8 DO I=1,24 TMPKE(I,(J-1)*3+1)=TMPKE2(I,J) TMPKE(I,(J-1)*3+2)=TMPKE2(I,J+8) TMPKE(I,(J-1)*3+3)=TMPKE2(I,J+16) END DO END DO CALL ZDANUL(TMPKE2,576) DO I=1,8 DO J=1,24 TMPKE2((I-1)*3+1,J)=TMPKE(I,J) TMPKE2((I-1)*3+2,J)=TMPKE(I+8,J) TMPKE2((I-1)*3+3,J)=TMPKE(I+16,J) END DO END DO DO J=1,24 DO I=1,24 RE(I,J)=RE(I,J) + 4.*AJAC*PXG5(IP)*TMPKE2(I,J) END DO END DO END DO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C MATRICE DE STABILISATION : PAS DE BOUCLE SUR LES POINTS DE GAUSS C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ON DEFINIE LES CONSTANTS DE B_2_CHAPEAU_BAR C CALL ZDANUL(CB2BAR,36) CB2BAR(1,1) = 1.D0 CB2BAR(1,2) = 1.D0 CB2BAR(1,3) = ZERO CB2BAR(1,4) = ZERO CB2BAR(1,6) = ZERO C CB2BAR(2,1) = ZERO CB2BAR(2,2) = ZERO CB2BAR(2,3) = 1.D0 CB2BAR(2,4) = 1.D0 CB2BAR(2,6) = ZERO C CB2BAR(3,1) = ZERO CB2BAR(3,2) = ZERO CB2BAR(3,3) = ZERO CB2BAR(3,4) = ZERO CB2BAR(3,6) = ZERO C CB2BAR(4,1) = ZERO CB2BAR(4,2) = ZERO CB2BAR(4,3) = ZERO CB2BAR(4,4) = ZERO C CB2BAR(5,4) = ZERO CB2BAR(5,5) = ZERO CB2BAR(5,6) = 1.D0 C CB2BAR(6,2) = ZERO CB2BAR(6,5) = ZERO CB2BAR(6,6) = 1.D0 C C ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4) C DO I=1,24 XELOC(I) = XE(I) END DO C C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES! C CALL CUBROT(XELOC,RR2) CALL VECTRO(RR2,XELOC,1) CALL CUBEBB(XELOC,B,VOL) C C CALCUL DES FONCTIONS DE FORME GS C XXGB = X * GB C DO J=1,3 DO IA=1,4 XXGB(J,IA) = HOUXGB(XELOC(J),IA) END DO END DO C C GS = (BBB) * XXGB C DO I=1,8 DO J=1,4 GS(I,J)=0.D0 END DO END DO DO J=1,3 DO IA=1,4 DO I=1,8 GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA) END DO END DO END DO C C GS = GB - GS C DO I=1,4 DO J=1,8 GS(J,I)= (GB(J,I) - GS(J,I))*UNS8 END DO END DO C C CALCUL DE XXVB = X * VB C XXVB(1)= - XELOC(1) + XELOC(4) + XELOC(7) - XELOC(10) & - XELOC(13) + XELOC(16) + XELOC(19) - XELOC(22) XXVB(2)= - XELOC(2) - XELOC(5) + XELOC(8) + XELOC(11) & - XELOC(14) - XELOC(17) + XELOC(20) + XELOC(23) XXVB(3)= - XELOC(3) - XELOC(6) - XELOC(9) - XELOC(12) & + XELOC(15) + XELOC(18) + XELOC(21) + XELOC(24) C C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES C HIJ(1) = UNS3*XXVB(2)*XXVB(3)/XXVB(1) HIJ(2) = UNS3*XXVB(1)*XXVB(3)/XXVB(2) HIJ(3) = UNS3*XXVB(2)*XXVB(1)/XXVB(3) HIJ(4) = UNS3*XXVB(3) HIJ(5) = UNS3*XXVB(1) HIJ(6) = UNS3*XXVB(2) C C CALCUL DES COEFS A METTRE DANS KIJ POUR COMPOSER KSTAB C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C ANCIENNE PROG A FARID (PLEXUS): C XK11_1=(XLAMBDA+2*XMU)*HIJ(1)+XMU*HIJ(2) C XK11_2=UNS3*((XLAMBDA+2*XMU)*HIJ(1)+XMU*(HIJ(2)+HIJ(3))) C XK22_1=(XLAMBDA+2*XMU)*HIJ(2)+XMU*HIJ(1) C XK22_2=UNS3*((XLAMBDA+2*XMU)*HIJ(2)+XMU*(HIJ(1)+HIJ(3))) C XK33_1=XMU*(HIJ(1)+HIJ(2)) C XK33_2=UNS3*(PROPEL(1)*HIJ(3)+XMU*(HIJ(1)+HIJ(2))) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ICI IL FAUT PRENDRE LA MOYENNE DES MODULES D'YOUNG TANGENTS C C POUR LE SHB8 PLASTIQUE ET NON_LINEAIRE, METTRE IETAN A 1 DANS RIGID.FF C SINON METTRE IETAN A 0 C IF (IETAN.EQ.1) THEN XYOUNGT = (EYG(1)+EYG(2)+EYG(3)+EYG(4)+EYG(5))/5. ELSE XYOUNGT = PROPEL(1) ENDIF C C POUR RESOUDRE LE CISAILLEMENT TRANSVERSE: C XYOUNGT = XYOUNGT*0.5*(COELA1+COELA2) XLAMBDA = XYOUNGT*PROPEL(2)/(1-PROPEL(2)*PROPEL(2)) XMU = 0.5 * XYOUNGT / (1+PROPEL(2)) XLAMBDA2= XLAMBDA + 2.D0*XMU C ON MET LES LIGNES SUIVANTES EN COMMENTAIRE C C XK11_1 = (XLAMBDA+2.D0*XMU)*HIJ(1) C XK11_2 = UNS3*((XLAMBDA+2.D0*XMU)*HIJ(1)) C XK22_1 = (XLAMBDA+2.D0*XMU)*HIJ(2) C XK22_2 = UNS3*((XLAMBDA+2.D0*XMU)*HIJ(2)) C XK33_1 = 0. C XK11_1 = (CB2BAR(1,1)*(XLAMBDA2*CB2BAR(1,1) + XLAMBDA*CB2BAR(2,1)) &+ CB2BAR(2,1)*(XLAMBDA*CB2BAR(1,1) + XLAMBDA2*CB2BAR(2,1)) &+ XYOUNGT*CB2BAR(3,1)**2)*HIJ(1) + XMU*HIJ(2)*CB2BAR(4,1)**2 XK11_2 = (CB2BAR(1,2)*(XLAMBDA2*CB2BAR(1,2) + XLAMBDA*CB2BAR(2,2)) &+ CB2BAR(2,2)*(XLAMBDA*CB2BAR(1,2) + XLAMBDA2*CB2BAR(2,2)) C &+ XYOUNGT*CB2BAR(4,2)**2)*HIJ(1)*UNS3 &+ XYOUNGT*CB2BAR(3,2)**2)*HIJ(1)*UNS3 &+ XMU*UNS3*(HIJ(2)*CB2BAR(4,2)**2 + HIJ(3)*CB2BAR(6,2)**2) C XK22_1 = (CB2BAR(1,3)*(XLAMBDA2*CB2BAR(1,3) + XLAMBDA*CB2BAR(2,3)) &+ CB2BAR(2,3)*(XLAMBDA*CB2BAR(1,3) + XLAMBDA2*CB2BAR(2,3)) &+ XYOUNGT*CB2BAR(3,3)**2)*HIJ(2) + XMU*HIJ(1)*CB2BAR(4,3)**2 XK22_2 = (CB2BAR(1,4)*(XLAMBDA2*CB2BAR(1,4) + XLAMBDA*CB2BAR(2,4)) &+ CB2BAR(2,4)*(XLAMBDA*CB2BAR(1,4) + XLAMBDA2*CB2BAR(2,4)) &+ XYOUNGT*CB2BAR(3,4)**2)*HIJ(2)*UNS3 &+ XMU*UNS3*(HIJ(1)*CB2BAR(4,4)**2 + HIJ(3)*CB2BAR(5,4)**2) C C XK33_1 = XMU*(HIJ(2)*CB2BAR(5,5)**2 + HIJ(1)*CB2BAR(6,5)) XK33_1 = XMU*(HIJ(2)*CB2BAR(5,5)**2 + HIJ(1)*CB2BAR(6,5)**2) XK33_2 = (CB2BAR(1,6)*(XLAMBDA2*CB2BAR(1,6) + XLAMBDA*CB2BAR(2,6)) &+ CB2BAR(2,6)*(XLAMBDA*CB2BAR(1,6) + XLAMBDA2*CB2BAR(2,6)) &+ XYOUNGT*CB2BAR(3,6)**2)*UNS3*HIJ(3) &+ XMU*UNS3*(HIJ(2)*CB2BAR(5,6)**2 + HIJ(1)*CB2BAR(6,6)**2) C C ON MET LES SUIVANTES EN COMMENTAIRE POUR ANNULER LES TERMES CROISES C C XK13_1 = XYOUNGT*CB2BAR(3,1)*HIJ(6) C XK13_2 = XMU*CB2BAR(6,5)*HIJ(6) C XK23_1 = XYOUNGT*CB2BAR(3,3)*HIJ(5) C XK23_2 = XMU*CB2BAR(5,5)*HIJ(5) C XK31_1 = XMU*CB2BAR(6,5)*HIJ(6) C XK31_2 = XYOUNGT*CB2BAR(3,1)*HIJ(6) C XK32_1 = XMU*CB2BAR(5,5)*HIJ(5) C XK32_2 = XYOUNGT*CB2BAR(3,3)*HIJ(5) C XK13_1 = ZERO XK13_2 = ZERO XK23_1 = ZERO XK23_2 = ZERO XK31_1 = ZERO XK31_2 = ZERO XK32_1 = ZERO XK32_2 = ZERO C C ON MET LES LIGNES SUIVANTES EN COMMENTAIRE C C IF (IRDC.EQ.1) THEN C COMPORTEMENT SHB8 PLEXUS C XK33_2 = XYOUNGT*HIJ(3)*UNS3 C ENDIF C IF (IRDC.EQ.2) THEN C COMPORTEMENT C.P. XK33_2 = 0. ENDIF C IF (IRDC.EQ.3) THEN C COMPORTEMENT LOI TRIDIM MMC 3D XK11_1 = XCOOEF*(1-XNU)*HIJ(1) XK11_2 = UNS3*(XCOOEF*(1-XNU)*HIJ(1)) XK22_1 = XCOOEF*(1-XNU)*HIJ(2) XK22_2 = UNS3*(XCOOEF*(1-XNU)*HIJ(2)) XK33_1 = 0. XK33_2 = XCOOEF*(1-XNU)*HIJ(3)*UNS3 ENDIF C CALL ZDANUL(XK11,64) CALL ZDANUL(XK22,64) CALL ZDANUL(XK12,64) CALL ZDANUL(XK21,64) CALL ZDANUL(XK13,64) CALL ZDANUL(XK23,64) CALL ZDANUL(XK31,64) CALL ZDANUL(XK32,64) CALL ZDANUL(XK33,64) C DO J=1,8 DO I=1,8 C C IL FAUT CALCULER K11 K22 K33 K13 K23 K31 K32 MATRICES 8*8 C XK11(I,J)=XK11_1*GS(I,3)*GS(J,3)+XK11_2*GS(I,4)*GS(J,4) XK22(I,J)=XK22_1*GS(I,3)*GS(J,3)+XK22_2*GS(I,4)*GS(J,4) XK33(I,J)=XK33_1*GS(I,3)*GS(J,3)+XK33_2*GS(I,4)*GS(J,4) C XK13(I,J)=XK13_1*GS(I,3)*GS(J,1)+XK13_2*GS(I,1)*GS(J,3) XK23(I,J)=XK23_1*GS(I,3)*GS(J,2)+XK23_2*GS(I,2)*GS(J,3) XK31(I,J)=XK31_1*GS(I,3)*GS(J,1)+XK31_2*GS(I,1)*GS(J,3) XK32(I,J)=XK32_1*GS(I,3)*GS(J,2)+XK32_2*GS(I,2)*GS(J,3) END DO END DO C C ASSEMBLAGE DE KSTAB C CALL ASKSTB(XKSTAB,XK11,XK22,XK33,XK12,XK21,XK13,XK23, & XK31,XK32) C C C REMISE EN ORDRE DE KSTAB C DO J=1,8 DO I=1,24 TMPKE(I,(J-1)*3+1)=XKSTAB(I,J) TMPKE(I,(J-1)*3+2)=XKSTAB(I,J+8) TMPKE(I,(J-1)*3+3)=XKSTAB(I,J+16) END DO END DO CALL ZDANUL(XKSTAB,576) DO I=1,8 DO J=1,24 XKSTAB((I-1)*3+1,J)=TMPKE(I,J) XKSTAB((I-1)*3+2,J)=TMPKE(I+8,J) XKSTAB((I-1)*3+3,J)=TMPKE(I+16,J) END DO END DO C C IL FAUT REPASSER DANS LE REPERE GLOBAL AVEC RR2^T . KSTAB . RR2 C EN FAIT C'EST RR2. KSTAB .RR2^T CAR RR2 RANGEE PAR LIGNES C CALL ASKSGL(XKSTAB,RR2) C C AJOUT DE KSTAB A KE C DO J=1,24 DO I=1,24 RE(I,J)=RE(I,J)+XKSTAB(I,J) END DO END DO ENDIF C IF(ICLE.EQ.7.or.icle.eq.11)THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C icle=11 calcul de epsi C C ON CALCULE LES CONTRAINTES : SORTIE DANS OUT(30) C C CONTRAINTES LOCALES DANS CHAQUE COUCHE C C SUR LA CONFIGURATION 1 C C LE DEPLACEMENT NODAL A L'AIR D'ETRE DANS WORK(1 A 24) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C UE: INCREMENT DE DEPLACEMENT NODAL, REPERE GLOBAL C C XE: DEBUT DU PAS irdc=3 DO J=1,8 DO I=1,3 UE(I,J)=WORK((J-1)*3+I) END DO END DO IGRDDEP = nint(D(1)) LAG = nint(PROPEL(3)) C C ON DEFINIT CMATLOC LOI MODIFIEE SHB8 C if(icle.eq.7) then CALL ZDANUL(CMATLOC,36) XLAMBDA = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2)) XMU = 0.5 * PROPEL(1) / (1+PROPEL(2)) CMATLOC(1,1) = XLAMBDA+2*XMU CMATLOC(2,2) = XLAMBDA+2*XMU IF (IRDC.EQ.1) THEN C COMPORTEMENT SHB8 PLEXUS CMATLOC(3,3) = PROPEL(1) ENDIF C IF (IRDC.EQ.2) THEN C COMPORTEMENT C.P. CMATLOC(3,3) = 0 ENDIF C CMATLOC(1,2) = XLAMBDA CMATLOC(2,1) = XLAMBDA CMATLOC(4,4) = XMU CMATLOC(5,5) = XMU CMATLOC(6,6) = XMU C IF (IRDC.EQ.3) THEN C COMPORTEMENT LOI TRIDIM MMC 3D XNU = PROPEL(2) XCOOEF = PROPEL(1)/((1+XNU)*(1-2*XNU)) CMATLOC(1,1) = (1-XNU)*XCOOEF CMATLOC(2,2) = (1-XNU)*XCOOEF CMATLOC(3,3) = (1-XNU)*XCOOEF CMATLOC(1,2) = XNU*XCOOEF CMATLOC(2,1) = XNU*XCOOEF CMATLOC(1,3) = XNU*XCOOEF CMATLOC(3,1) = XNU*XCOOEF CMATLOC(2,3) = XNU*XCOOEF CMATLOC(3,2) = XNU*XCOOEF CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF ENDIF endif C C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE C BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP C BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP C BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP C CALL CUBKSI(5,XXG5,BKSIP) C DO IP=1,5 C C DEFINITION DES 4 POINTS COQUES C ZETA = XXG5(IP) ZLAMB = 0.5*(1.-ZETA) DO J=1,3 DO I=1,4 XCOQ(J,I) = ZLAMB*XE((I-1)*3+J) & + (1.-ZLAMB)*XE((I-1+4)*3+J) END DO END DO C C CALCUL DE PPP 3*3 PASSAGE DE GLOBAL A LOCAL,COQUE C XCENT : COORD GLOBAL DU CENTRE DE L'ELEMENT C CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID) C C CALCUL DE B : U_GLOBAL ---> EPS_GLOBAL C CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC) C C CALCUL DE EPS DANS LE REPERE GLOBAL: 1 POUR DEFORMATIONS LINEAIRES C 2 POUR TERMES CARRES EN PLUS DO I=1,6 DEPS(I)=0. ENDDO IF (LAG.EQ.1) THEN C ON AJOUTE LA PARTIE NON-LINEAIRE DE EPS CALL DSDX3D(2,B,UE,DEPS,DUSDX,8) ELSE CALL DSDX3D(1,B,UE,DEPS,DUSDX,8) ENDIF C C SORTIE DE DUSDX DANS PROPEL(1 A 9 * 5 PT DE GAUSS) C POUR UTILISATION ULTERIEURE DANS Q8PKCN_SHB8 C CALL AEQBT(PPPT,PPP,3,3) RR12(1,1) = DUSDX(1) RR12(1,2) = DUSDX(2) RR12(1,3) = DUSDX(3) RR12(2,1) = DUSDX(4) RR12(2,2) = DUSDX(5) RR12(2,3) = DUSDX(6) RR12(3,1) = DUSDX(7) RR12(3,2) = DUSDX(8) RR12(3,3) = DUSDX(9) * CALL MULMAT(3,3,3,PPPT,RR12,RR2) CALL MULMAT(RR2,PPPT,RR12,3,3,3) * CALL MULMAT(3,3,3,RR2,PPP,RR12) CALL MULMAT(RR12,RR2,PPP,3,3,3) DUSDX(1) = RR12(1,1) DUSDX(2) = RR12(1,2) DUSDX(3) = RR12(1,3) DUSDX(4) = RR12(2,1) DUSDX(5) = RR12(2,2) DUSDX(6) = RR12(2,3) DUSDX(7) = RR12(3,1) DUSDX(8) = RR12(3,2) DUSDX(9) = RR12(3,3) DO I=1,9 PROPEL(I+(IP-1)*9)=DUSDX(I) ENDDO DO I=1,6 DEPS_LOC(I) = 0. SIG_LOC(I) = 0. ENDDO CALL CHRP3D(PPP,DEPS,DEPS_LOC,2) C C CALCUL DE SIGMA DANS LE REPERE LOCAL C * CALL MULMAT(6,6,1,CMATLOC,DEPS_LOC,SIG_LOC) if(icle.eq.7) then CALL MULMAT(SIG_LOC,CMATLOC,DEPS_LOC,6,1,6) endif C C CONTRAINTES ECRITES SOUS LA FORME: C [SIG] = [S_11, S_22, S_33, S_12, S_23, S_13] if(icle.eq.7) then DO I=1,6 C ON LAISSE LES CONTRAINTES DANS LE REPERE LOCAL POUR LA PLASTICITE OUT((IP-1)*6+I)=SIG_LOC(I) END DO elseif(icle.eq.11)then DO I=1,6 C ON LAISSE LES CONTRAINTES DANS LE REPERE LOCAL POUR LA PLASTICITE OUT((IP-1)*6+I)=DEPS_LOC(I) END DO endif END DO ENDIF IF(icle.eq.12) then C propel(3)=1 on veut epsi à partir de sigma C propel(3)=2 on veut sigma à partir de epsi C en entrée work contient soit epsilon soit sigma C work(1,2 ... contrainte 1 du point 1, contrainte 2 sdu point 1 .. C propel(1)=young C propel(2)=nu C propel(3) =1 si on veut des contraintes C =2 si on veut des epsilons C C out en sortie out(1,2 .. epsi d1 du point 1 , epsi 2 du point 1 .. irdc=3 icas=nint(propel(3)) C C ON DEFINIT CMATLOC LOI MODIFIEE SHB8 C CALL ZDANUL(CMATLOC,36) XLAMBDA = PROPEL(1)*PROPEL(2)/(1-PROPEL(2)*PROPEL(2)) XMU = 0.5 * PROPEL(1) / (1+PROPEL(2)) CMATLOC(1,1) = XLAMBDA+2*XMU CMATLOC(2,2) = XLAMBDA+2*XMU IF (IRDC.EQ.1) THEN C COMPORTEMENT SHB8 PLEXUS CMATLOC(3,3) = PROPEL(1) ENDIF C IF (IRDC.EQ.2) THEN C COMPORTEMENT C.P. CMATLOC(3,3) = 0 ENDIF C CMATLOC(1,2) = XLAMBDA CMATLOC(2,1) = XLAMBDA CMATLOC(4,4) = XMU CMATLOC(5,5) = XMU CMATLOC(6,6) = XMU C IF (IRDC.EQ.3) THEN C COMPORTEMENT LOI TRIDIM MMC 3D XNU = PROPEL(2) XCOOEF = PROPEL(1)/((1+XNU)*(1-2*XNU)) CMATLOC(1,1) = (1-XNU)*XCOOEF CMATLOC(2,2) = (1-XNU)*XCOOEF CMATLOC(3,3) = (1-XNU)*XCOOEF CMATLOC(1,2) = XNU*XCOOEF CMATLOC(2,1) = XNU*XCOOEF CMATLOC(1,3) = XNU*XCOOEF CMATLOC(3,1) = XNU*XCOOEF CMATLOC(2,3) = XNU*XCOOEF CMATLOC(3,2) = XNU*XCOOEF CMATLOC(4,4) = (1-2*XNU)*0.5*XCOOEF CMATLOC(5,5) = (1-2*XNU)*0.5*XCOOEF CMATLOC(6,6) = (1-2*XNU)*0.5*XCOOEF ENDIF if(icas.eq.2) then call invalm(cmatloc,6,3 , ier,1.d-12) CMATLOC(4,4) =1.d0/CMATLOC(4,4) CMATLOC(5,5) =1.d0/CMATLOC(5,5) CMATLOC(6,6) =1.d0/CMATLOC(6,6) endif write(6,*) ' cmatloc ' write(6,*) (cmatloc(1,iop),iop=1,6) write(6,*) (cmatloc(2,iop),iop=1,6) write(6,*) (cmatloc(3,iop),iop=1,6) write(6,*) (cmatloc(4,iop),iop=1,6) write(6,*) (cmatloc(5,iop),iop=1,6) write(6,*) (cmatloc(6,iop),iop=1,6) write(6,*) 'epsi ' C DO IP=1,5 C C DEFINITION DES 4 POINTS COQUES C write(6,*) (work((ip-1)*6 +iop),iop=1,6) C C CALCUL DE SIGMA DANS LE REPERE LOCAL C CALL MULMAT(out((ip-1)*6+1),CMATLOC,work((ip-1)*6+1),6,1,6) END DO ENDIF IF(ICLE.EQ.8)THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ON CALCULE BSIGMA: SORTIE DANS OUT(24) C C ENTREE DE SIGMA DANS WORK(DIM=30) C C ENTREE DU MATERIAU DANS D(1) ET D(2) C C ENTREE DE QIALPHA PAS PRECEDENT C C DANS RE(1 A 12) C C QIALPHA(J,I)=RE(J+(I-1)*3) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC LAG = nint(D(3)) IRDC=3 CALL ZDANUL(XKSIGTMP2,64) DO J=1,8 DO I=1,3 F(I,J) = 0.D0 END DO END DO C C CALCUL DE BKSIP(3,8,IP) DANS REPERE DE REFERENCE C BKSIP(1,*,IP) = VECTEUR BX AU POINT GAUSS IP C BKSIP(2,*,IP) = VECTEUR BY AU POINT GAUSS IP C BKSIP(3,*,IP) = VECTEUR BZ AU POINT GAUSS IP C CALL CUBKSI(5,XXG5,BKSIP) DO IP=1,5 C C RECHERCHE DE SIGMA DU POINT DE GAUSS GLOBAL C DO I=1,6 C C C'EST LES CONTRAINTES LOCALES POUR POUVOIR TRAITER LA PLASTICITE AVANT C LES CONTRAINTES SONT PASSEES DANS LA CONFI. A LA FIN DU PAS DANS Q8PKCN C SIG_LOC(I)=WORK((IP-1)*6+I) END DO ZETA = XXG5(IP) ZLAMB = 0.5*(1.-ZETA) DO J=1,3 DO I=1,4 XCOQ(J,I) = ZLAMB*XE((I-1)*3+J) & + (1.-ZLAMB)*XE((I-1+4)*3+J) END DO END DO CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,xBID) C C PASSAGE DES CONTRAINTES AU REPERE GLOBAL C CALL CHRP3D(PPP,SIG_LOC,SIGMAG,1) C CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC) C C CALCUL DE BQ.SIGMA SI LAGRANGIEN TOTAL C IF (LAG.EQ.1) THEN CALL ASKSIG(SIGMAG,B,XKSIGTMP1) DO J=1,8 DO I=1,8 XKSIGTMP2(I,J) = XKSIGTMP2(I,J) & + 4D0*AJAC*PXG5(IP)*XKSIGTMP1(I,J) ENDDO ENDDO ENDIF C C CALCUL DE B.SIGMA EN GLOBAL C POIDS = 4D0 * AJAC * PXG5(IP) DO K=1,8 F(1,K) = F(1,K) + POIDS * & (B(1,K)*SIGMAG(1)+B(2,K)*SIGMAG(4)+B(3,K)*SIGMAG(6)) F(2,K) = F(2,K) + POIDS * & (B(1,K)*SIGMAG(4)+B(2,K)*SIGMAG(2)+B(3,K)*SIGMAG(5)) F(3,K) = F(3,K) + POIDS * & (B(1,K)*SIGMAG(6)+B(2,K)*SIGMAG(5)+B(3,K)*SIGMAG(3)) END DO C C SI LAGRANGIEN TOTAL: AJOUT DE FQ A F C END DO IF (LAG.EQ.1) THEN CALL ZDANUL(TMPKE,576) DO KK=1,3 DO I=1,8 DO J=1,8 TMPKE(I+(KK-1)*8,J+(KK-1)*8) = XKSIGTMP2(I,J) ENDDO ENDDO ENDDO CALL ZDANUL(TMPKE2,576) DO J=1,8 DO I=1,24 TMPKE2(I,(J-1)*3+1)=TMPKE(I,J) TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8) TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16) END DO END DO CALL ZDANUL(TMPKE,576) DO I=1,8 DO J=1,24 TMPKE((I-1)*3+1,J)=TMPKE2(I,J) TMPKE((I-1)*3+2,J)=TMPKE2(I+8,J) TMPKE((I-1)*3+3,J)=TMPKE2(I+16,J) END DO END DO * CALL MULMAT(24,24,1,TMPKE,PROPEL,FQ) CALL MULMAT(FQ,TMPKE,PROPEL,24,1,24) DO K=1,8 F(1,K) = F(1,K) + FQ((K-1)*3+1) F(2,K) = F(2,K) + FQ((K-1)*3+2) F(3,K) = F(3,K) + FQ((K-1)*3+3) END DO ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C ON CALCULE LA STABILISATION POUR LE CALCUL DE B.SIGMA C UE EST DANS PROPEL! C ON A BESOIN DU MATERIAU AUSSI! CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C XYOUNG = D(1) XNU = D(2) XLAMBDA = XYOUNG*XNU/(1-XNU*XNU) XMU = 0.5 * XYOUNG / (1+XNU) C C DEPLACEMENT NODAL, REPERE GLOBAL : DANS PROPEL(1 A 24) C DO I=1,24 XE_LOC(I) = XE(I) XE_LOCP(I) = XE(I)-PROPEL(I) XE_LOC12(I) = XE(I)-0.5*PROPEL(I) ENDDO C C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES! C CALL CUBROT(XE_LOC,RR2) CALL CUBROT(XE_LOCP,RR1) CALL CUBROT(XE_LOC12,RR12) CALL VECTRO(RR2,XE_LOC,1) CALL VECTRO(RR1,XE_LOCP,1) CALL VECTRO(RR12,XE_LOC12,1) C C ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4) C ORIGINE DE XE_LOC12 PAS BONNE, MAIS PAS GRAVE : DIFFERENTIELLES DANS B C CALL CUBEBB(XE_LOC12,B,VOL) C C C CALCUL DEPLACEMENT DEFORMANT DANS REPERE 1/2 PAS: C DO I=1,24 UDEF(I)=XE_LOC(I)-XE_LOCP(I) ENDDO C C CALCUL DES FONCTIONS DE FORME GS C XXGB = X * GB C DO IA=1,4 DO J=1,3 XXGB(J,IA) = HOUXGB(XE_LOC12(J),IA) END DO END DO C C GS = (BBB) * XXGB C DO J=1,4 DO I=1,8 GS(I,J)=0.D0 END DO END DO DO J=1,3 DO IA=1,4 DO I=1,8 GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA) END DO END DO END DO C C GS = GB - GS C DO I=1,4 DO J=1,8 GS(J,I)= (GB(J,I) - GS(J,I))*UNS8 END DO END DO C C CALCUL DE XXVB = X * VB C XXVB(1)= - XE_LOC12(1) + XE_LOC12(4) + XE_LOC12(7) - XE_LOC12(10) & - XE_LOC12(13) + XE_LOC12(16) + XE_LOC12(19) - XE_LOC12(22) XXVB(2)= - XE_LOC12(2) - XE_LOC12(5) + XE_LOC12(8) + XE_LOC12(11) & - XE_LOC12(14) - XE_LOC12(17) + XE_LOC12(20) + XE_LOC12(23) XXVB(3)= - XE_LOC12(3) - XE_LOC12(6) - XE_LOC12(9) - XE_LOC12(12) & + XE_LOC12(15) + XE_LOC12(18) + XE_LOC12(21) + XE_LOC12(24) C C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES C HIJ(1) = UNS3*XXVB(2)*XXVB(3) / XXVB(1) HIJ(2) = UNS3*XXVB(1)*XXVB(3) / XXVB(2) HIJ(3) = UNS3*XXVB(2)*XXVB(1) / XXVB(3) HIJ(4) = UNS3*XXVB(3) HIJ(5) = UNS3*XXVB(1) HIJ(6) = UNS3*XXVB(2) C C CALCUL DES DEFORMATIONS GENERALISEES C DO IA=1,4 DO J=1,3 AUX=0.D0 DO I=1,8 AUX = AUX + GS(I,IA)*UDEF((I-1)*3+J) END DO PQIALPH(J,IA) = AUX END DO END DO C C CALCUL DES CONTRAINTES GENERALISEES C DO I=1,4 DO J=1,3 QIALPHA(J,I)=RE(J+(I-1)*3,1) END DO END DO CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C PROGRAMMATION FARID (PLEXUS) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C QIALPHA(1,1) = QIALPHA(1,1) C QIALPHA(2,2) = QIALPHA(2,2) C QIALPHA(3,3) = QIALPHA(3,3) + XMU*((HIJ(2)+HIJ(1))*PQIALPH(3,3)) C QIALPHA(1,2) = QIALPHA(1,2) C QIALPHA(2,3) = QIALPHA(2,3) + ((XLAMBDA+2.D0*XMU)*HIJ(2) C & + XMU*HIJ(1))*PQIALPH(2,3) C QIALPHA(1,3) = QIALPHA(1,3) + ((XLAMBDA+2.D0*XMU)*HIJ(1) C & + XMU*HIJ(2))*PQIALPH(1,3) C QIALPHA(2,1) = QIALPHA(2,1) C QIALPHA(3,2) = QIALPHA(3,2) C QIALPHA(3,1) = QIALPHA(3,1) C QIALPHA(1,4) = QIALPHA(1,4) + UNS3*((XLAMBDA+2.D0*XMU)*HIJ(1) C & + XMU*(HIJ(2)+HIJ(3)))*PQIALPH(1,4) C QIALPHA(2,4) = QIALPHA(2,4) + UNS3*((XLAMBDA+2.D0*XMU)*HIJ(2) C & + XMU*(HIJ(1)+HIJ(3)))*PQIALPH(2,4) C QIALPHA(3,4) = QIALPHA(3,4) + UNS3*(XYOUNG*HIJ(3) C & + XMU*(HIJ(1)+HIJ(2)))*PQIALPH(3,4) CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF (IRDC.EQ.1) THEN QIALPHA(1,1) = QIALPHA(1,1) QIALPHA(2,2) = QIALPHA(2,2) QIALPHA(3,3) = QIALPHA(3,3) QIALPHA(1,2) = QIALPHA(1,2) QIALPHA(2,3) = QIALPHA(2,3) & + ((XLAMBDA+2.D0*XMU)*HIJ(2))*PQIALPH(2,3) QIALPHA(1,3) = QIALPHA(1,3) & + ((XLAMBDA+2.D0*XMU)*HIJ(1))*PQIALPH(1,3) QIALPHA(2,1) = QIALPHA(2,1) QIALPHA(3,2) = QIALPHA(3,2) QIALPHA(3,1) = QIALPHA(3,1) QIALPHA(1,4) = QIALPHA(1,4) & + UNS3*((XLAMBDA+2*XMU)*HIJ(1))*PQIALPH(1,4) QIALPHA(2,4) = QIALPHA(2,4) & + UNS3*((XLAMBDA+2*XMU)*HIJ(2))*PQIALPH(2,4) C PROJECTION LEGAY POUR LE TERME QIALPHA(3,4) C COMPORTEMENT SHB8 PLEXUS CC QIALPHA(3,4) = QIALPHA(3,4) + XYOUNG*HIJ(3)*UNS3*PQIALPH(3,4) C PROJECTION 28 LE TERME QIALPHA(3,4) EST MODIFIE QIALPHA(3,4)=QIALPHA(3,4)+UNS3*(HIJ(1)+HIJ(2))*XMU*PQIALPH(3,4) ENDIF C A MODIFIER POUR LA PROJECTION 28 (CAS DES C.P.) IF (IRDC.EQ.2) THEN QIALPHA(1,1) = QIALPHA(1,1) QIALPHA(2,2) = QIALPHA(2,2) QIALPHA(3,3) = QIALPHA(3,3) QIALPHA(1,2) = QIALPHA(1,2) QIALPHA(2,3) = QIALPHA(2,3) & + ((XLAMBDA+2.D0*XMU)*HIJ(2))*PQIALPH(2,3) QIALPHA(1,3) = QIALPHA(1,3) & + ((XLAMBDA+2.D0*XMU)*HIJ(1))*PQIALPH(1,3) QIALPHA(2,1) = QIALPHA(2,1) QIALPHA(3,2) = QIALPHA(3,2) QIALPHA(3,1) = QIALPHA(3,1) QIALPHA(1,4) = QIALPHA(1,4) & + UNS3*((XLAMBDA+2*XMU)*HIJ(1))*PQIALPH(1,4) QIALPHA(2,4) = QIALPHA(2,4) & + UNS3*((XLAMBDA+2*XMU)*HIJ(2))*PQIALPH(2,4) C COMPORTEMENT C.P. QIALPHA(3,4) = QIALPHA(3,4) ENDIF IF (IRDC.EQ.3) THEN C COMPORTEMENT LOI TRIDIM MMC 3D C A MODIFIER POUR LA PROJECTION 28 (CAS DE LA LOI 3D) XCOOEF = XYOUNG/((1+XNU)*(1-2*XNU)) QIALPHA(1,1) = QIALPHA(1,1) QIALPHA(2,2) = QIALPHA(2,2) QIALPHA(3,3) = QIALPHA(3,3) QIALPHA(1,2) = QIALPHA(1,2) QIALPHA(2,3) = QIALPHA(2,3) & + ((XCOOEF*(1-XNU))*HIJ(2))*PQIALPH(2,3) QIALPHA(1,3) = QIALPHA(1,3) & + ((XCOOEF*(1-XNU))*HIJ(1))*PQIALPH(1,3) QIALPHA(2,1) = QIALPHA(2,1) QIALPHA(3,2) = QIALPHA(3,2) QIALPHA(3,1) = QIALPHA(3,1) QIALPHA(1,4) = QIALPHA(1,4) & + UNS3*((XCOOEF*(1-XNU))*HIJ(1))*PQIALPH(1,4) QIALPHA(2,4) = QIALPHA(2,4) & + UNS3*((XCOOEF*(1-XNU))*HIJ(2))*PQIALPH(2,4) QIALPHA(3,4) = QIALPHA(3,4) & + XCOOEF*(1-XNU)*HIJ(3)*UNS3*PQIALPH(3,4) ENDIF CCCCCCCCCCCCCCCCCCC C C SAUVEGARDE DES FORCES DE STABILISATION C DO I=1,4 DO J=1,3 RE(J+(I-1)*3,1)=QIALPHA(J,I) END DO END DO C C CALCUL DES FORCES DE HOURGLASS FHG C DO I=1,8 DO J=1,3 FHG(J,I)=0.D0 ENDDO END DO DO J=1,3 DO I=1,8 DO IA=1,4 FHG(J,I) = FHG(J,I) + QIALPHA(J,IA)*GS(I,IA) END DO END DO END DO C ON REPASSE AU REPERE GLOBAL DO I=1,3 DO J=1,8 FHG24((J-1)*3+I)=FHG(I,J) ENDDO ENDDO CALL VECTRO(RR12,FHG24,2) C C AJOUT DE LA STABILISATION AU B SIGMA DEJA CALCULE C DO J=1,3 DO I=1,8 F(J,I) = F(J,I) + FHG24((I-1)*3+J) ENDDO ENDDO C C ATTENTION A L'ORDRE DE OUT C DO I=1,3 DO J=1,8 OUT((J-1)*3+I) = F(I,J) END DO END DO ENDIF IF(ICLE.EQ.5)THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C ON CALCULE LES FORCES DE PRESSION C C C C ENTREES: PROPEL(1) : VALEUR DE LA PRESSION C C XE : GEOMETRIE ELEMENT C C WORK(1 A 3) : DIRECTION DE PRESSION C C PROPEL(2) : 0 RECHERCHE DE LA FACE C C 1 OU 2 : FACE INF. OU SUP. (POUR C C PRESSION SUIVEUSE) C C SORTIE: OUT FORCE DE PRESSION C C WORK(4) : NUMERO DE FACE POUR LA PRESSION C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC DO I=1,8 DO J=1,3 F(J,I)=0. END DO END DO C C LECTURE: PRESSION A METTRE AUX 4 NOEUDS, FACE ET DIRECTION C PRESS=PROPEL(1) IFACE=nint(PROPEL(2)) IF(IFACE.EQ.0)THEN DO I=1,3 DIREC(I)=WORK(I) END DO ENDIF C C COORD. DES POINTS MILIEUX C DO J=1,3 X56(J) =( XE((5-1)*3+J)+XE((6-1)*3+J) )/2. X67(J) =( XE((6-1)*3+J)+XE((7-1)*3+J) )/2. X78(J) =( XE((7-1)*3+J)+XE((8-1)*3+J) )/2. X85(J) =( XE((8-1)*3+J)+XE((5-1)*3+J) )/2. X12(J) =( XE((1-1)*3+J)+XE((2-1)*3+J) )/2. X23(J) =( XE((2-1)*3+J)+XE((3-1)*3+J) )/2. X34(J) =( XE((3-1)*3+J)+XE((4-1)*3+J) )/2. X41(J) =( XE((4-1)*3+J)+XE((1-1)*3+J) )/2. END DO C C NORMALE RENTRANTE FACE SUP: V56_78 ^ V85_67 C IF(IFACE.EQ.0.OR.IFACE.EQ.2)THEN DO J=1,3 V56_78(J)=X78(J)-X56(J) V85_67(J)=X67(J)-X85(J) END DO VSUP(1) = V56_78(2) * V85_67(3) - V56_78(3) * V85_67(2) VSUP(2) = V56_78(3) * V85_67(1) - V56_78(1) * V85_67(3) VSUP(3) = V56_78(1) * V85_67(2) - V56_78(2) * V85_67(1) ENDIF C C NORMALE RENTRANTE FACE INF: V41_23 ^ V12_34 C IF(IFACE.EQ.0.OR.IFACE.EQ.1)THEN DO J=1,3 V12_34(J)=X34(J)-X12(J) V41_23(J)=X23(J)-X41(J) END DO VINF(1) = V41_23(2) * V12_34(3) - V41_23(3) * V12_34(2) VINF(2) = V41_23(3) * V12_34(1) - V41_23(1) * V12_34(3) VINF(3) = V41_23(1) * V12_34(2) - V41_23(2) * V12_34(1) ENDIF IF(IFACE.EQ.0) THEN C C RECHERCHE DE LA BONNE FACE: PREMIER PASSAGE C SCALSUP = 0. SCALINF = 0. DO I=1,3 SCALSUP = SCALSUP+VSUP(I)*DIREC(I) SCALINF = SCALINF+VINF(I)*DIREC(I) END DO C C TEST ERREUR C IF((SCALSUP*SCALINF).GE.0)THEN IF ((SCALSUP*SCALINF).EQ.0)THEN WRITE(6,*)'** ERREUR SUR LES SURFACES **' WRITE(6,*)'** DE SHB8 POUR PRESSION **' WRITE(6,*)'** PRODUIT SCALAIRE ÉGAL A ZERO **' WRITE(6,*)'** PROBLÈME DE DIRECTION DE PRESSION **' ENDIF IF ((SCALSUP*SCALINF).GT.0)THEN WRITE(6,*)'** ERREUR SUR LES SURFACES **' WRITE(6,*)'** DE SHB8 POUR PRESSION **' WRITE(6,*)'** LES 2 FACES SONT POSSIBLES **' ENDIF STOP ENDIF C C SI SCAL <0 : PRESSION SUR CETTE FACE C IF(SCALSUP.LT.0)THEN WORK(4) = 1 IFACE = 1 ENDIF IF(SCALINF.LT.0)THEN WORK(4) = 2 IFACE = 2 ENDIF ENDIF IF(IFACE.EQ.2)THEN C C BOUCLE SUR LES 4 NOEUDS DE LA FACE SUPERIEURE C DO I=5,8 DO J=1,3 F(J,I)= PRESS * VSUP(J) / 4. END DO END DO ENDIF IF(IFACE.EQ.1)THEN C C BOUCLE SUR LES 4 NOEUDS DE LA FACE INFERIEURE C DO I=1,4 DO J=1,3 F(J,I)= PRESS * VINF(J) / 4. END DO END DO ENDIF C C ATTENTION A L'ORDRE DE OUT C DO J=1,8 DO I=1,3 OUT((J-1)*3+I)=F(I,J) END DO END DO ENDIF IF(ICLE.EQ.9)THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C CALCUL DE K.SIGMA C C EN ENTREE DANS WORK : SIGMA NORMALEMENT LONGUEUR 30 C C PROPEL(1): 1 POUR RE=RE+KSIGMA C C PROPEL(1): 0 POUR RE=KSIGMA C C SORTIE : RE(24*24)=RE(24*24)+KSIGMA C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C CALCUL DE B (1 2 3) AUX 5 POINTS DE GAUSS C CALL CUBKSI(5,XXG5,BKSIP) DO J=1,8 DO I=1,8 XKSIGTMP2(I,J)=0. ENDDO ENDDO C C C CALCUL DE LA CONTRAINTE MOYENNE SUR LES 5 POINTS DE GAUSS POUR LE C KSIGMA DE STABILISATION CALCULE PLUS LOIN C DO J=1,6 SIGMOY(J)=0.D0 ENDDO C DO J=1,6 DO IP=1,5 SIGMOY(J)=SIGMOY(J)+WORK((IP-1)*6+J) ENDDO SIGMOY(J)=UNS5*SIGMOY(J) ENDDO C C DEBUT DE LA BOUCLE SUR LES 5 PTS GAUSS C DO IP=1,5 C C CALCUL DE B C CALL CUCALB(BKSIP(1,1,IP),XE,B,AJAC) C C CALCUL DE MATRICE DE PASSAGE POUR POUVOIR CALCULER LES CONTRAINTES DANS C LE REPERE GLOBAL C DO I=1,6 C C'EST LES CONTRAINTES LOCALES POUR POUVOIR TRAITER LA PLASTICITE AVANT SIG_LOC(I)=WORK((IP-1)*6+I) END DO ZETA = XXG5(IP) ZLAMB = 0.5*(1.-ZETA) DO J=1,3 DO I=1,4 XCOQ(J,I) = ZLAMB*XE((I-1)*3+J) & + (1.-ZLAMB)*XE((I-1+4)*3+J) END DO END DO CALL RLOSHB(XCOQ,XCENT,PPP,XL,XXX,YYY,XBID) C C PASSAGE DES CONTRAINTES AU REPERE GLOBAL C CALL CHRP3D(PPP,SIG_LOC,SIGMAG,1) CALL ASKSIG(SIGMAG,B,XKSIGTMP1) DO J=1,8 DO I=1,8 XKSIGTMP2(I,J) = XKSIGTMP2(I,J) & + 4D0*AJAC*PXG5(IP)*XKSIGTMP1(I,J) ENDDO ENDDO ENDDO CALL ZDANUL(TMPKE,576) DO KK=1,3 DO I=1,8 DO J=1,8 TMPKE(I+(KK-1)*8,J+(KK-1)*8) = XKSIGTMP2(I,J) ENDDO ENDDO ENDDO C C ON MET DE L'ORDRE: C CALL ZDANUL(TMPKE2,576) DO J=1,8 DO I=1,24 TMPKE2(I,(J-1)*3+1)=TMPKE(I,J) TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8) TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16) END DO END DO CALL ZDANUL(TMPKE,576) DO I=1,8 DO J=1,24 TMPKE((I-1)*3+1,J)=TMPKE2(I,J) TMPKE((I-1)*3+2,J)=TMPKE2(I+8,J) TMPKE((I-1)*3+3,J)=TMPKE2(I+16,J) END DO END DO IPROPEL=nint(PROPEL(1)) IF(IPROPEL.EQ.0) CALL SHIFTD(TMPKE,RE,576) IF(IPROPEL.EQ.1) CALL AEQAPB(RE,TMPKE,576) C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C MATRICE DE STABILISATION POUR K.SIGMA: PAS DE BOUCLE SUR LES POINTS DE GAUSS C C C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C C ON A BESOIN DES VECTEURS GAMMA(8, ALPHA=1 A 4) = GS(8,4) C CCC DO I=1,24 CCC XELOC(I) = XE(I) CCC END DO C C REMARQUE: ATTENTION, RR, MATRICE DU CORROTATIONNEL EST RANGEE PAR LIGNES! C CCC CALL CUBEROT(XELOC,RR2) CCC CALL VECTROT(RR2,XELOC,1) CCC CALL CUBE_BB(XELOC,B,VOL) C C CALCUL DES FONCTIONS DE FORME GS C XXGB = X * GB C CCC DO J=1,3 CCC DO IA=1,4 CCC XXGB(J,IA) = HOUXGB(XELOC(J),IA) CCC END DO CCC END DO C C GS = (BBB) * XXGB C CCC DO I=1,8 CCC DO J=1,4 CCC GS(I,J)=0.D0 CCC END DO CCC END DO CCC DO J=1,3 CCC DO IA=1,4 CCC DO I=1,8 CCC GS(I,IA)=GS(I,IA) + B(J,I)*XXGB(J,IA) CCC END DO CCC END DO CCC END DO C C GS = GB - GS C CCC DO I=1,4 CCC DO J=1,8 CCC GS(J,I)= (GB(J,I) - GS(J,I))*UNS8 CCC END DO CCC END DO C C CALCUL DE XXVB = X * VB C CCC XXVB(1)= - XELOC(1) + XELOC(4) + XELOC(7) - XELOC(10) CCC & - XELOC(13) + XELOC(16) + XELOC(19) - XELOC(22) CCC XXVB(2)= - XELOC(2) - XELOC(5) + XELOC(8) + XELOC(11) CCC & - XELOC(14) - XELOC(17) + XELOC(20) + XELOC(23) CCC XXVB(3)= - XELOC(3) - XELOC(6) - XELOC(9) - XELOC(12) CCC & + XELOC(15) + XELOC(18) + XELOC(21) + XELOC(24) C C CALCUL DES RELATIONS CONTRAINTES ET DEFORMATIONS GENERALISEES C CCC HIJ(1) = UNS3*XXVB(2)*XXVB(3)/XXVB(1) CCC HIJ(2) = UNS3*XXVB(1)*XXVB(3)/XXVB(2) CCC HIJ(3) = UNS3*XXVB(2)*XXVB(1)/XXVB(3) CCC HIJ(4) = UNS3*XXVB(3) CCC HIJ(5) = UNS3*XXVB(1) CCC HIJ(6) = UNS3*XXVB(2) C C CALCUL DES COEFS A METTRE DANS KIJ POUR COMPOSER KSTAB C C C ICI IL FAUT PRENDRE LA MOYENNE DES CONTRAINTES SUR LES POINTS DE GAUSS C CECI EST DEJA CALCULE PLUS HAUT DANS LA VARIABLE SIGMOY(6) C C C POUR LA PROJECTION 28 ON OBTIENT C CCC XK11_1 = SIGMOY(1)*HIJ(1) CCC XK11_2 = SIGMOY(1)*HIJ(1)*UNS3 CCC XK11_3 = SIGMOY(6)*HIJ(6) CCC XK11_4 = SIGMOY(6)*HIJ(6) C CCC XK22_1 = SIGMOY(2)*HIJ(2) CCC XK22_2 = SIGMOY(2)*HIJ(2)*UNS3 CCC XK22_3 = SIGMOY(5)*HIJ(5) CCC XK22_4 = SIGMOY(5)*HIJ(5) C CCC XK33_1 = ZERO CCC XK33_2 = (SIGMOY(1)*HIJ(1)+SIGMOY(2)*HIJ(2))*UNS3 C CCC XK13_1 = ZERO CCC XK13_2 = ZERO CCC XK23_1 = ZERO CCC XK23_2 = ZERO CCC XK31_1 = ZERO CCC XK31_2 = ZERO CCC XK32_1 = ZERO CCC XK32_2 = ZERO C C CCC CALL ZDANUL(XK11,64) CCC CALL ZDANUL(XK22,64) CCC CALL ZDANUL(XK12,64) CCC CALL ZDANUL(XK21,64) CCC CALL ZDANUL(XK13,64) CCC CALL ZDANUL(XK23,64) CCC CALL ZDANUL(XK31,64) CCC CALL ZDANUL(XK32,64) CCC CALL ZDANUL(XK33,64) C CCC DO J=1,8 CCC DO I=1,8 C C IL FAUT CALCULER K11 K22 K33 K13 K23 K31 K32 MATRICES 8*8 C CCC XK11(I,J)=XK11_1*GS(I,3)*GS(J,3)+XK11_2*GS(I,4)*GS(J,4) C & +XK11_3*GS(I,3)*GS(J,1)+XK11_4*GS(I,1)*GS(J,3) CCC XK22(I,J)=XK22_1*GS(I,3)*GS(J,3)+XK22_2*GS(I,4)*GS(J,4) C & +XK22_3*GS(I,3)*GS(J,2)+XK22_4*GS(I,2)*GS(J,3) CCC XK33(I,J)=XK33_1*GS(I,3)*GS(J,3)+XK33_2*GS(I,4)*GS(J,4) C CCC XK13(I,J)=XK13_1*GS(I,3)*GS(J,1)+XK13_2*GS(I,1)*GS(J,3) CCC XK23(I,J)=XK23_1*GS(I,3)*GS(J,2)+XK23_2*GS(I,2)*GS(J,3) CCC XK31(I,J)=XK31_1*GS(I,3)*GS(J,1)+XK31_2*GS(I,1)*GS(J,3) CCC XK32(I,J)=XK32_1*GS(I,3)*GS(J,2)+XK32_2*GS(I,2)*GS(J,3) C CCC END DO CCC END DO C C ASSEMBLAGE DE KSTAB C CCC CALL ASSE_KSTAB(XKSTAB,XK11,XK22,XK33,XK12,XK21,XK13,XK23, CCC & XK31,XK32) C C C REMISE EN ORDRE DE KSTAB C CCC DO J=1,8 CCC DO I=1,24 CCC TMPKE(I,(J-1)*3+1)=XKSTAB(I,J) CCC TMPKE(I,(J-1)*3+2)=XKSTAB(I,J+8) CCC TMPKE(I,(J-1)*3+3)=XKSTAB(I,J+16) CCC END DO CCC END DO CCC CALL ZDANUL(XKSTAB,576) CCC DO I=1,8 CCC DO J=1,24 CCC XKSTAB((I-1)*3+1,J)=TMPKE(I,J) CCC XKSTAB((I-1)*3+2,J)=TMPKE(I+8,J) CCC XKSTAB((I-1)*3+3,J)=TMPKE(I+16,J) CCC END DO CCC END DO C C IL FAUT REPASSER DANS LE REPERE GLOBAL AVEC RR2^T . KSTAB . RR2 C EN FAIT C'EST RR2. KSTAB .RR2^T CAR RR2 RANGEE PAR LIGNES C CCC CALL ASSE_KSTAB_GL(XKSTAB,RR2) C C AJOUT DE KSTAB A KE C CCC DO J=1,24 CCC DO I=1,24 CCC RE(I,J)=RE(I,J)+XKSTAB(I,J) CCC END DO CCC END DO C C C ENDIF CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC IF(ICLE.EQ.3)THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CALCUL DE LA MATRICE DE MASSE C C ENTREE: PROPEL(1) = MASSE VOLUMIQUE C C XE(1 A 14)= COORDONNEES DES NOEUDS C C SORTIE: RE(24,24) = MATRICE MASSE C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC XRO = PROPEL(1) CALL ZDANUL(RE,576) C C CALCUL DU PLAN MOYEN: C DO J=1,3 DO I=1,4 XCOQ(J,I) = 0.5*XE((I-1)*3+J)+0.5*XE((I-1+4)*3+J) END DO END DO C CALCUL DE LA SURFACE MOYENNE: DO J=1,3 X12(J) = (XCOQ(J,1)+XCOQ(J,2))*0.5 X34(J) = (XCOQ(J,3)+XCOQ(J,4))*0.5 X23(J) = (XCOQ(J,2)+XCOQ(J,3))*0.5 X41(J) = (XCOQ(J,4)+XCOQ(J,1))*0.5 END DO DO J=1,3 V12_34(J) = X34(J)-X12(J) V41_23(J) = X23(J)-X41(J) END DO C C CALCUL DE V41_23 ^ V12_34: C VINF(1) = V41_23(2) * V12_34(3) - V41_23(3) * V12_34(2) VINF(2) = V41_23(3) * V12_34(1) - V41_23(1) * V12_34(3) VINF(3) = V41_23(1) * V12_34(2) - V41_23(2) * V12_34(1) C C CALCUL DU VOLUME DE L'ELEMENT: C CALL CUBKSI(5,XXG5,BKSIP) CALL CUCALB(BKSIP(1,1,3),XE,B,AJAC) XVOLUME = AJAC*8. CALL ZDANUL(XMASSE,64) XXXMASSE = 0. DO I=1,8 XMASSE(I,I) = (8./27.)*AJAC*XRO XXXMASSE = XXXMASSE+XMASSE(I,I) ENDDO DO I=1,8 DO J=1,8 IF(I.NE.J)THEN XMASSE(I,J) = (1./64.)*(2.+VB(I,1)*VB(J,1)*2./3.) & *(2.+VB(I,2)*VB(J,2)*2./3.) & *(2.+VB(I,3)*VB(J,3)*2./3.) & *AJAC*XRO XXXMASSE = XXXMASSE+XMASSE(I,I) ENDIF ENDDO ENDDO CALL ZDANUL(TMPKE,576) DO K=1,3 DO I=1,8 DO J=1,8 TMPKE(I+(K-1)*8,J+(K-1)*8) = XMASSE(I,J) ENDDO ENDDO ENDDO DO J=1,8 DO I=1,24 TMPKE2(I,(J-1)*3+1)=TMPKE(I,J) TMPKE2(I,(J-1)*3+2)=TMPKE(I,J+8) TMPKE2(I,(J-1)*3+3)=TMPKE(I,J+16) END DO END DO CALL ZDANUL(RE,576) DO I=1,8 DO J=1,24 RE((I-1)*3+1,J)=TMPKE2(I,J) RE((I-1)*3+2,J)=TMPKE2(I+8,J) RE((I-1)*3+3,J)=TMPKE2(I+16,J) END DO END DO ENDIF IF(ICLE.EQ.10)THEN CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C C CALCUL DE LA MATRICE KP C C ENTREE: PRESSION DANS PROPEL(1) C C KTAFL DANS PROPEL(2) C C RE=RE+KP C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PRESSI = PROPEL(1) IFACE = nint(PROPEL(2)) ITEST = 0 IF (IFACE.EQ.1) THEN C C DEFINITION DE LA NORMALE N0, FACE INFERIEURE C DO I=1,4 DO J=1,3 XCOQ(J,I) = XE((I-1)*3+J) END DO END DO DO J=1,3 XCOQ14(J) = (XCOQ(J,1)+XCOQ(J,4))/2. XCOQ23(J) = (XCOQ(J,2)+XCOQ(J,3))/2. XCOQ34(J) = (XCOQ(J,3)+XCOQ(J,4))/2. XCOQ12(J) = (XCOQ(J,1)+XCOQ(J,2))/2. ENDDO DO J=1,3 XCOQ12_34(J) = XCOQ34(J)-XCOQ12(J) XCOQ14_23(J) = XCOQ23(J)-XCOQ14(J) ENDDO C XNORMALE(1) = XCOQ14_23(2)*XCOQ12_34(3) & -XCOQ14_23(3)*XCOQ12_34(2) XNORMALE(2) = XCOQ14_23(3)*XCOQ12_34(1) & -XCOQ14_23(1)*XCOQ12_34(3) XNORMALE(3) = XCOQ14_23(1)*XCOQ12_34(2) & -XCOQ14_23(2)*XCOQ12_34(1) XNORMALE(1) = -XNORMALE(1) XNORMALE(2) = -XNORMALE(2) XNORMALE(3) = -XNORMALE(3) C XXG5(1) = -1. XXG5(2) = 0. XXG5(3) = 0. XXG5(4) = 0. XXG5(5) = 0. ITEST = 1 ENDIF IF (IFACE.EQ.2) THEN C C DEFINITION DE LA NORMALE N0, FACE SUPERIEURE C DO I=1,4 DO J=1,3 XCOQ(J,I) = XE(12+(I-1)*3+J) END DO END DO DO J=1,3 XCOQ14(J)=(XCOQ(J,1)+XCOQ(J,4))/2. XCOQ23(J)=(XCOQ(J,2)+XCOQ(J,3))/2. XCOQ34(J)=(XCOQ(J,3)+XCOQ(J,4))/2. XCOQ12(J)=(XCOQ(J,1)+XCOQ(J,2))/2. ENDDO DO J=1,3 XCOQ12_34(J)=XCOQ34(J)-XCOQ12(J) XCOQ14_23(J)=XCOQ23(J)-XCOQ14(J) ENDDO C XNORMALE(1)=XCOQ14_23(2)*XCOQ12_34(3) & -XCOQ14_23(3)*XCOQ12_34(2) XNORMALE(2)=XCOQ14_23(3)*XCOQ12_34(1) & -XCOQ14_23(1)*XCOQ12_34(3) XNORMALE(3)=XCOQ14_23(1)*XCOQ12_34(2) & -XCOQ14_23(2)*XCOQ12_34(1) C XXG5(1)=+1. XXG5(2)=0. XXG5(3)=0. XXG5(4)=0. XXG5(5)=0. ITEST=1 ENDIF IF(ITEST.EQ.0)THEN WRITE(6,*)'ARRET DANS SHB8, CALCUL DE KP, FACE NON DEFINIE' STOP ENDIF C C CALCUL DE B SUR LA FACE INF OU SUP, SUIVANT XXG5 C C DANS UE: ON MET L'EQUIVALENT DE BKSIP, CALCULER C POUR LES 4 NOEUDS DE LA FACE IFACE C CALL CUKSKP(IFACE,UE) C XNORM = SQRT(XNORMALE(1)*XNORMALE(1)+XNORMALE(2)*XNORMALE(2) & +XNORMALE(3)*XNORMALE(3)) DO J=1,3 XNORMALE(J) = XNORMALE(J)/XNORM ENDDO CALL CUCBKP(IFACE,UE,XE,BKP,AJAC) CALL ASSKP2(IFACE,BKP,XNORMALE,XKP) C C REMISE EN ORDRE DE KP C DO I=1,24 DO J=1,8 TMPKE(I,(J-1)*3+1)=XKP(I,J) TMPKE(I,(J-1)*3+2)=XKP(I,J+8) TMPKE(I,(J-1)*3+3)=XKP(I,J+16) END DO END DO CALL ZDANUL(XKP,576) DO I=1,8 DO J=1,24 XKP((I-1)*3+1,J)=TMPKE(I,J) XKP((I-1)*3+2,J)=TMPKE(I+8,J) XKP((I-1)*3+3,J)=TMPKE(I+16,J) END DO END DO DO I=1,24 DO J=1,24 RE(I,J)=RE(I,J)+XKP(I,J)*XNORM*PRESSI/8. ENDDO ENDDO ENDIF RETURN END