C PIOCAX SOURCE AS218369 10/01/29 21:15:18 6623 SUBROUTINE PIOCAX(NBNN,IDIM,TAB1,NCOELE,NBPTEL,IPMINT,XE1,XE2, 1 TABA,MRACC,SH1,TAB,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 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 IPMINT = POINTEUR DES FONCTIONS DE FORME C TABA = pointeur tableau avec ddl de saut c MRACC = pointeur tableau de raccourci pour les C enrichissements elementaires 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) C -INC SMLREEL -INC SMINTE DIMENSION TAB1(NBPTEL,*) DIMENSION TAB(NBPTEL,*) DIMENSION XE1(3,*),XE2(3,*),SH1(6,*) *as xfem 2010_01_13 PARAMETER (NBNNMAX=20) SEGMENT MRACC INTEGER TLREEL(NBNN) ENDSEGMENT * ddl de saut enrichissement SEGMENT TABA REAL*8 TABA1(IDIM,NBNN),TABA2(IDIM,NBNN) ENDSEGMENT * phi et H aux nnoeuds ; phi et H au point de Gauss courant DIMENSION xphi1(NBNNMAX),xh1(NBNNMAX),TABPHG(2),tabdh(nbnnMAX) *fin as xfem 2010_01_13 C C TABLEAUX DE TRAVAIL DIMENSIONNES ICI C DIMENSION XJAC(3,3),FAC(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 IN/1,2,3,1,1,2/ DATA JN/1,2,3,2,3,3/ C DATA ITAB(1,1),ITAB(1,2),ITAB(1,3)/1,4,5/ DATA ITAB(2,1),ITAB(2,2),ITAB(2,3)/4,2,6/ DATA ITAB(3,1),ITAB(3,2),ITAB(3,3)/5,6,3/ C KERRE=0 MINTE = IPMINT NBSH = SHPTOT(/2) IDIM1=IDIM+1 C C MISE A ZERO DES CONTRAINTES OU DES DEFORMATIONS C DO 50 IB=1,NCOELE DO 50 IA=1,NBPTEL TAB(IA,IB)=0.D0 50 CONTINUE C C BOUCLE SUR LES POINTS DE GAUSS C DO 130 IC=1,NBPTEL tabphg(2)=0.D0 *as xfem 2010_01_13 * Initialisation de SH1(Ni, Ni,x, Ni,y) do i3=1,nbnn xh1(i3)=0.D0 do i4=1,IDIM1 SH1(i4,i3)=SHPTOT(i4,i3,IC) enddo enddo * Calcul de H et phi aux noeuds et au point de Gauss IC do 131 i1=1,nbnn mlree1=tlreel(i1) if(mlree1.eq.0) goto 131 tabphg(1)=0.D0 do i2=1,nbnn xphi1(i2)=mlree1.PROG(i2) if (abs(xphi1(i2)).lt.1.d-7) then xh1(i2)=0.D0 else xh1(i2)=sign(1.d0,xphi1(i2)) endif tabphg(1)=tabphg(1)+SH1(1,i2)*xphi1(i2) enddo if (abs(tabphg(1)).lt.1.d-7) then tabphg(2)=0.D0 else tabphg(2)=sign(1.d0,tabphg(1)) endif 131 continue * Calcul des H(x)-H(xi) : do i3=1,nbnn tabdh(i3)=tabphg(2)-xh1(i3) enddo * Calcul de SH1 : call jacobix(XE1,TABA1,TABDH,SH1,IDIM,NBNN,DJac) C C CALCUL DE LA MATRICE F C CALL ZERO(XJAC,3,3) DO 140 ID=1,NBNN DO 140 IE=1,IDIM * r_z = XE2(IE,ID) r_z = XE2(IE,ID)+(tabdh(ID)*TABA2(IE,ID)) DO 140 IF=1,IDIM XJAC(IE,IF)=XJAC(IE,IF) + SH1(IF+1,ID)*r_z 140 CONTINUE *fin as xfem 2010_01_13 IF(IDIM.EQ.2) THEN XJAC(3,3)=1.D0 IF(IFOU.EQ.0) THEN C CCCCCCCCCCCCC CAS AXISYMETRIQUE C R1=0. R2=0. DO 150 ID=1,NBNN R1=R1+SH1(1,ID)*XE1(1,ID) R2=R2+SH1(1,ID)*XE2(1,ID) 150 CONTINUE XJAC(3,3)=R2/(R1+1.E-20) ENDIF ENDIF C C GO TO (500,600),KCAS C C KCAS=1 CAS DES CONTRAINTES C ---------------------------- C 500 CONTINUE C C CCCCCCCCCCCC CALCUL DE DETERMINANT DE F C IF(IDIM.EQ.2) THEN DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) DETF = DETF * XJAC (3,3) ENDIF 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)) ENDIF DETF=1./(DETF+1.E-20) C C CALCUL DES CONTRAINTES DE CAUCHY C DO 160 ID=1,NCOELE IND=IN(ID) JND=JN(ID) DO 170 IE=1,IDIM DO 170 IF=1,IDIM ICO=ITAB(IE,IF) TAB(IC,ID)=TAB1(IC,ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF 1 +TAB(IC,ID) 170 CONTINUE 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 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 INVALM(XJAC,3,3,KERRE,1.D-9) 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) DO 270 IE=1,IDIM DO 270 IF=1,IDIM ICO=ITAB(IE,IF) TAB(IC,ID)=TAB(IC,ID) + 1 FAC(ICO)*TAB1(IC,ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID) 270 CONTINUE 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) ENDIF C 130 CONTINUE C RETURN END