1 : C PIOCAX SOURCE AS218369 10/01/29 21:15:18 6623 2 : SUBROUTINE PIOCAX(NBNN,IDIM,TAB1,NCOELE,NBPTEL,IPMINT,XE1,XE2, 3 : 1 TABA,MRACC,SH1,TAB,IFOU,KCAS,KERRE) 4 : C======================================================================= 5 : C 6 : C TRANSFORME LES CONTRAINTES DE PIOLA KIRCHHOFF EN CONTRAINTES DE 7 : C CAUCHY 8 : C ENTREE 9 : C ------- 10 : C NBNN = NOMBRE DE POINTS PAR ELEMENTS 11 : C IDIM = DIMENSION DE L ESPACE SUPPORT 12 : C 13 : C TAB1(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE PIOLA KIRCHHOFF 14 : C 15 : C NCOELE = NOMBRE DE COMPOSTS TABLEAU DES CONTRAINTES 16 : C 17 : C NBPTEL = NOMBRE DE POINTS DE GAUSS 18 : C IPMINT = POINTEUR DES FONCTIONS DE FORME 19 : C TABA = pointeur tableau avec ddl de saut 20 : c MRACC = pointeur tableau de raccourci pour les 21 : C enrichissements elementaires 22 : C 23 : C KCAS = 1 SI CONTRAINTES, 2 SI DEFORMATIONS 24 : C 25 : C TABLEAUX DE TRAVAIL 26 : C-------------------- 27 : C XE1(3,NBNN) = COORDONNEES CORRESPONDANT A LA CONFIGURATION DEPART 28 : C 29 : C XE2(3,NBNN) = COORDONNEES CORRESPONDANT A LA CONFIGURATION ACTUEL 30 : C 31 : C SH1(6,NBNN) = FONCTIONS DE FORME EN UN POINT DE GAUSS 32 : C 33 : C SORTIES 34 : C--------- 35 : C TAB(NBPTEL,NCOELE) =TABLEAU DES CONTRAINTES DE CAUCHY 36 : C 37 : C 38 : C AOUT 85 39 : C MODIF PEGON FEV 90 CAS BIDIM 40 : C PASSAGE AUX NOUVEAUX CHAMELEMS PAR P.DOWLATYARI 12/4/91 41 : C 42 : C======================================================================= 43 : IMPLICIT INTEGER(I-N) 44 : IMPLICIT REAL*8(A-H,O-Z) 45 : C 46 : -INC SMLREEL 47 : -INC SMINTE 48 : 49 : DIMENSION TAB1(NBPTEL,*) 50 : DIMENSION TAB(NBPTEL,*) 51 : DIMENSION XE1(3,*),XE2(3,*),SH1(6,*) 52 : *as xfem 2010_01_13 53 : PARAMETER (NBNNMAX=20) 54 : SEGMENT MRACC 55 : INTEGER TLREEL(NBNN) 56 : ENDSEGMENT 57 : * ddl de saut enrichissement 58 : SEGMENT TABA 59 : REAL*8 TABA1(IDIM,NBNN),TABA2(IDIM,NBNN) 60 : ENDSEGMENT 61 : 62 : * phi et H aux nnoeuds ; phi et H au point de Gauss courant 63 : DIMENSION xphi1(NBNNMAX),xh1(NBNNMAX),TABPHG(2),tabdh(nbnnMAX) 64 : *fin as xfem 2010_01_13 65 : C 66 : C TABLEAUX DE TRAVAIL DIMENSIONNES ICI 67 : C 68 : DIMENSION XJAC(3,3),FAC(6) 69 : C 70 : C TABLEAUX INDIQUANT LA CORRESPONDANCE ENTRE INDICES I,J ET NUMERO 71 : C DE LA COMPOSANTE DE CONTRAINTES OU DE DEFORMATIONS 72 : C 73 : DIMENSION IN(6),JN(6),ITAB(3,3) 74 : C 75 : DATA FAC/1.D0,1.D0,1.D0,0.5D0,0.5D0,0.5D0/ 76 : DATA IN/1,2,3,1,1,2/ 77 : DATA JN/1,2,3,2,3,3/ 78 : C 79 : DATA ITAB(1,1),ITAB(1,2),ITAB(1,3)/1,4,5/ 80 : DATA ITAB(2,1),ITAB(2,2),ITAB(2,3)/4,2,6/ 81 : DATA ITAB(3,1),ITAB(3,2),ITAB(3,3)/5,6,3/ 82 : C 83 : KERRE=0 84 : 85 : MINTE = IPMINT 86 : NBSH = SHPTOT(/2) 87 : 88 : IDIM1=IDIM+1 89 : C 90 : C MISE A ZERO DES CONTRAINTES OU DES DEFORMATIONS 91 : C 92 : DO 50 IB=1,NCOELE 93 : DO 50 IA=1,NBPTEL 94 : TAB(IA,IB)=0.D0 95 : 50 CONTINUE 96 : C 97 : C BOUCLE SUR LES POINTS DE GAUSS 98 : C 99 : DO 130 IC=1,NBPTEL 100 : 101 : tabphg(2)=0.D0 102 : *as xfem 2010_01_13 103 : * Initialisation de SH1(Ni, Ni,x, Ni,y) 104 : do i3=1,nbnn 105 : xh1(i3)=0.D0 106 : do i4=1,IDIM1 107 : SH1(i4,i3)=SHPTOT(i4,i3,IC) 108 : enddo 109 : enddo 110 : * Calcul de H et phi aux noeuds et au point de Gauss IC 111 : do 131 i1=1,nbnn 112 : mlree1=tlreel(i1) 113 : if(mlree1.eq.0) goto 131 114 : tabphg(1)=0.D0 115 : do i2=1,nbnn 116 : xphi1(i2)=mlree1.PROG(i2) 117 : if (abs(xphi1(i2)).lt.1.d-7) then 118 : xh1(i2)=0.D0 119 : else 120 : xh1(i2)=sign(1.d0,xphi1(i2)) 121 : endif 122 : tabphg(1)=tabphg(1)+SH1(1,i2)*xphi1(i2) 123 : enddo 124 : if (abs(tabphg(1)).lt.1.d-7) then 125 : tabphg(2)=0.D0 126 : else 127 : tabphg(2)=sign(1.d0,tabphg(1)) 128 : endif 129 : 130 : 131 continue 131 : * Calcul des H(x)-H(xi) : 132 : do i3=1,nbnn 133 : tabdh(i3)=tabphg(2)-xh1(i3) 134 : enddo 135 : * Calcul de SH1 : 135 : call jacobix(XE1,TABA1,TABDH,SH1,IDIM,NBNN,DJac) 136 : C 137 : C CALCUL DE LA MATRICE F 138 : C 138 : CALL ZERO(XJAC,3,3) 139 : DO 140 ID=1,NBNN 140 : DO 140 IE=1,IDIM 141 : * r_z = XE2(IE,ID) 142 : r_z = XE2(IE,ID)+(tabdh(ID)*TABA2(IE,ID)) 143 : DO 140 IF=1,IDIM 144 : XJAC(IE,IF)=XJAC(IE,IF) + SH1(IF+1,ID)*r_z 145 : 140 CONTINUE 146 : *fin as xfem 2010_01_13 147 : IF(IDIM.EQ.2) THEN 148 : XJAC(3,3)=1.D0 149 : IF(IFOU.EQ.0) THEN 150 : C 151 : CCCCCCCCCCCCC CAS AXISYMETRIQUE 152 : C 153 : R1=0. 154 : R2=0. 155 : DO 150 ID=1,NBNN 156 : R1=R1+SH1(1,ID)*XE1(1,ID) 157 : R2=R2+SH1(1,ID)*XE2(1,ID) 158 : 150 CONTINUE 159 : XJAC(3,3)=R2/(R1+1.E-20) 160 : ENDIF 161 : ENDIF 162 : C 163 : C 164 : GO TO (500,600),KCAS 165 : C 166 : C KCAS=1 CAS DES CONTRAINTES 167 : C ---------------------------- 168 : C 169 : 500 CONTINUE 170 : C 171 : C 172 : CCCCCCCCCCCC CALCUL DE DETERMINANT DE F 173 : C 174 : IF(IDIM.EQ.2) THEN 175 : DETF=XJAC(1,1)*XJAC(2,2)-XJAC(1,2)*XJAC(2,1) 176 : DETF = DETF * XJAC (3,3) 177 : ENDIF 178 : IF(IDIM.EQ.3) THEN 179 : DETF=XJAC(1,1)*(XJAC(2,2)*XJAC(3,3)-XJAC(3,2)*XJAC(2,3)) 180 : DETF=DETF-XJAC(2,1)*(XJAC(1,2)*XJAC(3,3)-XJAC(3,2)*XJAC(1,3)) 181 : DETF=DETF+XJAC(3,1)*(XJAC(1,2)*XJAC(2,3)-XJAC(1,3)*XJAC(2,2)) 182 : ENDIF 183 : DETF=1./(DETF+1.E-20) 184 : C 185 : C CALCUL DES CONTRAINTES DE CAUCHY 186 : C 187 : DO 160 ID=1,NCOELE 188 : IND=IN(ID) 189 : JND=JN(ID) 190 : DO 170 IE=1,IDIM 191 : DO 170 IF=1,IDIM 192 : ICO=ITAB(IE,IF) 193 : TAB(IC,ID)=TAB1(IC,ICO)*XJAC(IND,IE)*XJAC(JND,IF)*DETF 194 : 1 +TAB(IC,ID) 195 : 170 CONTINUE 196 : 160 CONTINUE 197 : C 198 : C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE 199 : C 200 : IF(IDIM.EQ.2) THEN 201 : TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3)*DETF 202 : ENDIF 203 : GO TO 130 204 : C 205 : C KCAS=2 CAS DES DEFORMATIONS 206 : C ----------------------------- 207 : C 208 : 600 CONTINUE 209 : C 210 : C 211 : CCCCCCCCCCCC CALCUL DE L'INVERSE DE F 212 : C 212 : CALL INVALM(XJAC,3,3,KERRE,1.D-9) 213 : IF(KERRE.NE.0) THEN 214 : WRITE(6,77881) ((XJAC(MI,MJ),MJ=1,3),MI=1,3) 215 : 77881 FORMAT(2X,' MATRICE SINGULIERE' /(3(1X,1PE12.5)/)) 216 : RETURN 217 : ENDIF 218 : C 219 : C CALCUL DES DEFORMATIONS 220 : C 221 : DO 260 ID=1,NCOELE 222 : IND=IN(ID) 223 : JND=JN(ID) 224 : DO 270 IE=1,IDIM 225 : DO 270 IF=1,IDIM 226 : ICO=ITAB(IE,IF) 227 : TAB(IC,ID)=TAB(IC,ID) + 228 : 1 FAC(ICO)*TAB1(IC,ICO)*XJAC(IE,IND)*XJAC(IF,JND)/FAC(ID) 229 : 270 CONTINUE 230 : 260 CONTINUE 231 : C 232 : C PEGON : ON NE FAIT PAS LA TRANSFORMATION SUR LA 3-EME COMPOSANTE 233 : C 234 : IF(IDIM.EQ.2) THEN 235 : TAB(IC,3)=TAB1(IC,3)*XJAC(3,3)*XJAC(3,3) 236 : ENDIF 237 : C 238 : 130 CONTINUE 239 : C 240 : RETURN 241 : END 242 : 243 :
© Cast3M 2003 - All rights reserved.
Disclaimer