C PIOCAF SOURCE PV090527 24/03/14 21:15:02 11870 SUBROUTINE PIOCAF(NBNN,nbsh,IDIM,TAB1,NCOELE,NBPTEL,SHP,XE1, 1 XE2,SH1,TAB,MWRK6,LHOOK, 2 IFOU,KCAS,KERRE) C======================================================================= C C TRANSFORME LES CONTRAINTES DE PIOLA KIRCHHOFF EN CONTRAINTES DE C CAUCHY C ENTREE C ------- C NBNN = NOMBRE DE POINTS PAR ELEMENTS C NBSH = NOMBRE DE fonctions de forme C IDIM = DIMENSION DE L ESPACE SUPPORT C C TAB1(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE PIOLA KIRCHHOFF C C NCOELE = NOMBRE DE COMPOSTS TABLEAU DES CONTRAINTES C C NBPTEL = NOMBRE DE POINTS DE GAUSS C SHP(6,NBNN,NBPTEL)= FONCTIONS DE FORME C C KCAS = 1 SI CONTRAINTES, 2 SI DEFORMATIONS C C TABLEAUX DE TRAVAIL C-------------------- C XE1(3,NBNN) = COORDONNEES CORRESPONDANT A LA CONFIGURATION DEPART C C XE2(3,NBNN) = COORDONNEES CORRESPONDANT A LA CONFIGURATION ACTUEL C C SH1(6,NBNN) = FONCTIONS DE FORME EN UN POINT DE GAUSS C C SORTIES C--------- C TAB(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE CAUCHY C C C AOUT 85 C MODIF PEGON FEV 90 CAS BIDIM C PASSAGE AUX NOUVEAUX CHAMELEMS PAR P.DOWLATYARI 12/4/91 C C======================================================================= IMPLICIT INTEGER(I-N) IMPLICIT REAL*8(A-H,O-Z) -INC CCREEL * SEGMENT MWRK6 INTEGER ITRES1(NBPTEL) REAL*8 PRODDI(NBPTEL,LHOO2),PRODDO(NBPTEL,LHOO2) REAL*8 DDHOOK(LHOOK,LHOOK),DDHOMU(LHOOK,LHOOK) REAL*8 VEC(LHOOK),VEC2(LHOOK) ENDSEGMENT * C dimension r_z(3) DIMENSION TAB1(NBPTEL,*) DIMENSION TAB(NBPTEL,*) *as xfem 2010_01_13 DIMENSION SHP(6,NBSH,*) * DIMENSION SHP(6,NBNN,*) *fin as xfem 2010_01_13 DIMENSION XE1(3,*),XE2(3,*),SH1(6,*) C C TABLEAUX DE TRAVAIL DIMENSIONNES ICI C DIMENSION XJAC(3,3),FAC(6),FACI(6) C C TABLEAUX INDIQUANT LA CORRESPONDANCE ENTRE INDICES I,J ET NUMERO C DE LA COMPOSANTE DE CONTRAINTES OU DE DEFORMATIONS C DIMENSION IN(6),JN(6),ITAB(3,3) C DATA FAC/1.D0,1.D0,1.D0,0.5D0,0.5D0,0.5D0/ DATA FACI/1.D0,1.D0,1.D0,2.D0,2.D0,2.D0/ DATA IN/1,2,3,1,1,2/ DATA JN/1,2,3,2,3,3/ C DATA ITAB(1,1),ITAB(2,1),ITAB(3,1)/1,4,5/ DATA ITAB(1,2),ITAB(2,2),ITAB(3,2)/4,2,6/ DATA ITAB(1,3),ITAB(2,3),ITAB(3,3)/5,6,3/ C KERRE=0 C C MISE A ZERO DES CONTRAINTES OU DES DEFORMATIONS C ** DO 50 IA=1,NBPTEL ** DO 50 IB=1,NCOELE ** TAB(IA,IB)=0.D0 **50 CONTINUE CALL ZERO(TAB,NBPTEL,NCOELE) C C BOUCLE SUR LES POINTS DE GAUSS C DO 130 IC=1,NBPTEL C if (nbsh.ne.nbnn) then CALL HPRIMEX(XE1,NBNN,IDIM,NBSH,SHP,IC,SH1,DJAC) else CALL HPRIME(XE1,NBNN,IDIM,SHP,IC,SH1,DJAC) endif C C CALCUL DE LA MATRICE F (=XJAC) C CALL ZERO(XJAC,3,3) DO IF=1,IDIM JF = IF + 1 r_z1 = 0.D0 r_z2 = 0.D0 r_z3 = 0.D0 do id=1,nbnn r_z1 = r_z1 + SH1(JF,ID)*XE2(1,ID) if (idim.ge.2) r_z2 = r_z2 + SH1(JF,ID)*XE2(2,ID) if (idim.ge.3) r_z3 = r_z3 + SH1(JF,ID)*XE2(3,ID) enddo if (idim.ge.1) XJAC(1,IF) = r_z1 if (idim.ge.2) XJAC(2,IF) = r_z2 if (idim.ge.3) XJAC(3,IF) = r_z3 enddo C CC - MODES DE CALCUL EN DEFORMATIONS "PLANES" GENERALISEES CC CAS 3D : IF (IDIM.EQ.3) THEN C RIEN FAIRE ! CC CAS 2D : ELSE IF(IDIM.EQ.2) THEN CCCC CAS 2D AXISYMETRIQUE IF (IFOU.EQ.0) THEN R1=0.D0 R2=0.D0 DO 150 ID=1,NBNN R1=R1+SH1(1,ID)*XE1(1,ID) R2=R2+SH1(1,ID)*XE2(1,ID) 150 CONTINUE if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then XJAC(3,3)=R2/R1 else xjac(3,3)=xgrand*xzprec endif ELSE CCCC CAS 2D PLAN XJAC(3,3)=1.D0 CCCC CAS 2D PLAN DEFO GENE * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(3,1) (cf. PIOCAP) IF (IFOU.EQ.-3) THEN XJAC(3,3) = XJAC(3,3) + XE2(3,1) ENDIF ENDIF CC CAS 1D : ELSE IF (IDIM.EQ.1) THEN CCCC CAS 1D PLAN IF (IFOU.GE.3 .AND. IFOU.LE.11) THEN XJAC(2,2) = 1.D0 XJAC(3,3) = 1.D0 * Rq : "Deplacement" UY du PTGENE est stocke dans XE2(2,1) (cf. PIOCAP) IF (IFOU.EQ.7 .OR. IFOU.EQ.8 .OR. IFOU.EQ.11) THEN XJAC(2,2) = XJAC(2,2) + XE2(2,1) ENDIF * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(3,1) (cf. PIOCAP) c* IF (IFOU.EQ.9 .OR. IFOU.EQ.10 .OR. IFOU.EQ.11) THEN IF (IFOU.GE.9) THEN XJAC(3,3) = XJAC(3,3) + XE2(3,1) ENDIF CCCC CAS 1D AXIS et SPHE ELSE CALL DISTRR(XE1,SH1,NBNN,R1) CALL DISTRR(XE2,SH1,NBNN,R2) if (r1.lt.-xpetit/xzprec.or.r1.gt.xpetit/xzprec)then FR2R1 = R2 /R1 else FR2R1=XGRAND*XZPREC endif XJAC(3,3) = FR2R1 CCCC CAS 1D SPHE IF (IFOU.EQ.15) THEN XJAC(2,2) = FR2R1 ELSE CCCC CAS 1D AXIS XJAC(2,2) = 1.D0 * Rq : "Deplacement" UZ du PTGENE est stocke dans XE2(2,1) (cf. PIOCAP) IF (IFOU.EQ.14) THEN XJAC(2,2) = XJAC(2,2) + XE2(2,1) ENDIF ENDIF ENDIF ENDIF C GO TO (500,600,700),KCAS C C KCAS=1 CAS DES CONTRAINTES C ---------------------------- C 500 CONTINUE C CCCCCCCCCCCC CALCUL DE L'INVERSE DU DETERMINANT DE F C IF (IDIM.EQ.3) THEN DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3)) DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3)) DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2)) ELSE IF (IDIM.EQ.2) THEN DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3) ELSE IF (IDIM.EQ.1) THEN DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3) ENDIF if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then DETF=1./(DETF) else DETF=XGRAND*xzprec endif C C CALCUL DES CONTRAINTES DE CAUCHY C DO 160 ID=1,NCOELE IND=IN(ID) JND=JN(ID) TABCD = TAB(IC,ID) DO 170 IE=1,IDIM TABAU1=0.D0 DO 171 IF=1,IDIM ICO=ITAB(IF,IE) TABAU1=TAB1(IC,ICO)*XJAC(JND,IF) 1 +TABAU1 171 CONTINUE TABAU1 = TABAU1 * XJAC(IND,IE) TABCD = TABCD + TABAU1 170 CONTINUE TAB(IC,ID)=TABCD*DETF + TAB(IC,ID) 160 CONTINUE C C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE C IF (IDIM.EQ.2) THEN TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF ELSE IF (IDIM.EQ.1) THEN TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF ENDIF GO TO 130 C C KCAS=2 CAS DES DEFORMATIONS C ----------------------------- C 600 CONTINUE C C CCCCCCCCCCCC CALCUL DE L'INVERSE DE F C CALL INVMA1(XJAC,3,3,KERRE) IF(KERRE.NE.0) THEN WRITE(6,77881) ((XJAC(MI,MJ),MJ=1,3),MI=1,3) 77881 FORMAT(2X,' MATRICE SINGULIERE' /(3(1X,1PE12.5)/)) RETURN ENDIF C C CALCUL DES DEFORMATIONS C DO 260 ID=1,NCOELE IND=IN(ID) JND=JN(ID) TABCD=0.D0 DO 271 IE=1,IDIM * unrolling de la boucle 270 pour mettre xjac en registre ** TABAU1=0.D0 ** DO 270 IF=1,IDIM ** ICO=ITAB(IF,IE) ** TABAU1=TABAU1+ ** 1 FAC(ICO)*TAB1(IC,ICO)*XJAC(IF,JND) *270 CONTINUE ICO1=ITAB(1,IE) tabau1= FAC(ICO1)*TAB1(IC,ICO1)*XJAC(1,JND) ICO2=ITAB(2,IE) tabau1= tabau1 + FAC(ICO2)*TAB1(IC,ICO2)*XJAC(2,JND) if(idim.eq.3) then ICO3=ITAB(3,IE) tabau1= tabau1 + FAC(ICO3)*TAB1(IC,ICO3)*XJAC(3,JND) endif TABCD=TABCD+TABAU1*XJAC(IE,IND) 271 CONTINUE TAB(IC,ID)=TAB(IC,ID)+TABCD*FACI(ID) 260 CONTINUE C C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE C IF(IDIM.EQ.2) THEN TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3) ELSE IF(IDIM.EQ.1) THEN TAB(IC,2)=TAB1(IC,2)*XJAC(2,2)*XJAC(2,2)*DETF TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF ENDIF C GO TO 130 C C KCAS=3 CAS DE LA MATRICE DE HOOKE C ----------------------------------- C 700 CONTINUE C CCCCCCCCCCCC CALCUL DE L'INVERSE DU DETERMINANT DE F C IF (IDIM.EQ.3) THEN DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3)) DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3)) DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2)) ELSE IF (IDIM.EQ.2) THEN DETF = ( XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) ) * XJAC(3,3) ELSE IF (IDIM.EQ.1) THEN DETF = XJAC(1,1) * XJAC(2,2) * XJAC(3,3) ENDIF if (detf.lt.-xpetit/xzprec.or.detf.gt.xpetit/xzprec)then DETF=1./(DETF) else DETF=XGRAND*xzprec endif C IJ=1 DO JJ=1,LHOOK DO II=1,LHOOK DDHOOK(II,JJ)=PRODDI(IC,IJ) IJ=IJ+1 enddo enddo * CALL ZERO(DDHOMU,LHOOK,LHOOK) KEY=2 C DO 720 LC=1,LHOOK CALL ZERO (VEC,LHOOK,1) DO 760 ID=1,LHOOK IND=IN(ID) JND=JN(ID) DO IE=1,IDIM DO IF=1,IDIM ICO=ITAB(IF,IE) IF(ICO.EQ.LC) THEN VEC(ID)=VEC(ID)+ & FAC(ICO)*XJAC(IE,IND)*XJAC(IF,JND) ENDIF enddo enddo VEC(ID)=VEC(ID)*FACI(ID) 760 CONTINUE C C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE C IF (IDIM.EQ.2) THEN IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3) * ELSE IF (IDIM.EQ.1) THEN IF(LC.EQ.2) VEC(2)=XJAC(2,2)*XJAC(2,2) IF(LC.EQ.3) VEC(3)=XJAC(3,3)*XJAC(3,3) ENDIF C CALL MATVE1(DDHOOK,VEC,LHOOK,LHOOK,VEC2,KEY) C DO 761 ID=1,LHOOK IND=IN(ID) JND=JN(ID) DO IE=1,IDIM DO IF=1,IDIM ICO=ITAB(IF,IE) DDHOMU(ID,LC)=VEC2(ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF 1 +DDHOMU(ID,LC) enddo enddo 761 CONTINUE C C ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE C IF (IDIM.EQ.2) THEN DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF ELSE IF (IDIM.EQ.1) THEN DDHOMU(2,LC)=VEC2(2)*XJAC(2,2)*XJAC(2,2)*DETF DDHOMU(3,LC)=VEC2(3)*XJAC(3,3)*XJAC(3,3)*DETF ENDIF 720 CONTINUE C IJ=1 DO JJ=1,LHOOK DO II=1,LHOOK PRODDO(IC,IJ)=DDHOMU(II,JJ) IJ=IJ+1 enddo enddo * * 130 CONTINUE RETURN END