1 : C PHAPRO SOURCE CHAT 11/03/16 21:28:55 6902 2 : subroutine phapro 3 : IMPLICIT INTEGER(I-N) 4 : IMPLICIT REAL*8(A-H,O-Z) 5 : * 6 : * Calcule la nouvelle proportions de phase 7 : * en entree 8 : * : le modele et son chamelem de materiau associé 9 : * : un champ de temperature proposée avec les mult de 10 : * lagrange 11 : * en sortie chamelem de proportion de phase mis a jour 12 : * 13 : -INC CCOPTIO 14 : -INC SMMODEL 15 : -INC SMCHAML 16 : -INC SMCHPOI 17 : -INC SMELEME 18 : -INC SMRIGID 19 : -INC SMCOORD 20 : -INC CCREEL 21 : 22 : segment icp 23 : integer icpr1(npt) 24 : endsegment 25 : pointeur icpr.icp 26 : segment icon(imax) 27 : segment ipos(ino) 28 : segment ipo2(ino) 29 : segment imoy(npart) 30 : segment idej(npart) 31 : segment qrest(ino) 32 : segment qtot(ino) 33 : segment ivoi(nvoi) 34 : segment qmoye 35 : real*8 qmoy(nbelem) 36 : endsegment 37 : xpre=xpetit 38 : umxpre=1.d0 - xpre 39 : * lecture du chamelem de source de proportions de phase 39 : call lirobj('MCHAML',mchelm,1,iretou) 40 : if(ierr.ne.0) return 41 : * lecture du chamelem de source de chaleur latente 41 : call lirobj('MCHAML',mchel1,1,iretou) 42 : if(ierr.ne.0) return 43 : * lecture des reactions 43 : call lirobj('CHPOINT',mchpo1,1,iretou) 44 : if(ierr.ne.0) return 45 : * recopie du mcheml entré 45 : call copie8(mchelm,ipoin2) 46 : mchelm=ipoin2 47 : * ON regarde d'abord le chpoint des reactions. s'il est vide il n'y a 48 : * rien a faire 49 : segact mchpo1 50 : if(mchpo1.ipchp(/1).eq.0) then 50 : call ecrobj('MCHAML',ipoin2) 51 : segdes mchpo1 52 : return 53 : endif 54 : * test du chpoint en entree 55 : if(mchpo1.ipchp(/1) .ne.1) then 55 : call erreur(180) 56 : return 57 : endif 58 : msoupo=mchpo1.ipchp(1) 59 : segact msoupo 60 : meleme=igeoc 61 : segact meleme 62 : mpoval=ipoval 63 : segact mpoval 64 : * creation du tableau des connexions a partir du chamelem de chaleur 65 : * latente 66 : segact mchel1 67 : npt = xcoor(/1)/(idim+1) 68 : segini icp 69 : ngrand=0 70 : do 520 ima=1,mchel1.imache(/1) 71 : meleme=mchel1.imache(ima) 72 : segact meleme 73 : ngrand = max(ngrand,num(/2)) 74 : do 521 nel=1,num(/2) 75 : do 521 nbn=1,num(/1) 76 : ia = num(nbn,nel) 77 : icpr1(ia)=icpr1(ia)+1 78 : 521 continue 79 : 520 continue 80 : imax=0 81 : ino=1 82 : ngrand=ngrand+1 83 : do 522 ic=1,icpr1(/1) 84 : if(icpr1(ic).ne.0) then 85 : ino=ino+1 86 : imax=imax+icpr1(ic) 87 : endif 88 : 522 continue 89 : segini icon,ipos,ipo2 90 : ia=1 91 : nvoi=0 92 : do 525 ic=1,icpr1(/1) 93 : if(icpr1(ic).ne.0) then 94 : ia = ia + 1 95 : nvoi = max(icpr1(ic),nvoi) 96 : ipos(ia)=ipos(ia-1) + icpr1(ic) 97 : icpr1(ic)=ia-1 98 : endif 99 : 525 continue 100 : do 523 ima = 1,mchel1.imache(/1) 101 : meleme= mchel1.imache(ima) 102 : nplu= ngrand*(ima-1) 103 : do 524 nel=1,num(/2) 104 : do 524 nbn=1,num(/1) 105 : ia = num(nbn,nel) 106 : ib=icpr1(ia) 107 : ipo2(ib)=ipo2(ib)+1 108 : id=ipos(ib)+ipo2(ib) 109 : icon(id)=nel+nplu 110 : 524 continue 111 : 523 continue 112 : nvoi=nvoi+1 113 : segini ivoi 114 : segsup ipo2 115 : * activation du chamelem de chaleur latente et de proportion de phase 116 : segact mchelm 117 : do 560 n1=1,imache(/1) 118 : meleme=imache(n1) 119 : mchaml=ichaml(n1) 120 : mcham1=mchel1.ichaml(n1) 121 : segact meleme,mchaml,mcham1 122 : do 561 n2=1,ielval(/1) 123 : melval=ielval(n2) 124 : melva1=mcham1.ielval(n2) 125 : segact melval*MOD,melva1*MOD 126 : 561 continue 127 : 560 continue 128 : * fabrication de la proportion de phase deja mangee par element 129 : ifois=0 130 : segini,icpr=icp 131 : 1000 continue 132 : * write(6,*) ' icpr 1 ' 133 : * write(6,*) ( icpr.icpr1(iop),iop=1,10) 134 : ifois=ifois+1 135 : npart=imache(/1) 136 : if(ifois.eq.2) then 137 : ino=qrest(/1) 138 : segsup qrest,qtot 139 : do 1005 IMO=1,imoy(/1) 140 : qmoye = imoy(imo) 141 : nbelem=qmoy(/1) 142 : segsup qmoye 143 : 1005 continue 144 : segsup imoy 145 : endif 146 : segini imoy,qrest,qtot 147 : do 570 n1=1,imache(/1) 148 : meleme=imache(n1) 149 : do 580 ico=1,mchel1.imache(/1) 150 : if(meleme.eq.mchel1.imache(ico)) go to 581 151 : 580 continue 151 : call erreur(21) 152 : 581 continue 153 : mcham1=mchel1.ichaml(ico) 154 : mchaml=ichaml(n1) 155 : nbelem = num(/2) 156 : segini qmoye 157 : imoy(n1)=qmoye 158 : do 571 inom=1,nomche(/2) 159 : if( nomche(inom).eq.'PPHA') go to 573 160 : 571 continue 161 : moterr(1:4)='PPHA' 161 : call erreur(236) 162 : return 163 : 573 continue 164 : do 583 inom1=1,mcham1.nomche(/2) 165 : if(mcham1.nomche(inom1).eq.'Q') go to 584 166 : 583 continue 167 : moterr(1:4)='Q' 167 : call erreur(236) 168 : return 169 : 584 continue 170 : melva1=mcham1.ielval(inom1) 171 : melval=ielval(inom) 172 : nbnn=num(/1) 173 : do 575 nbe=1,num(/2) 174 : qs=0.d0 175 : nbnn=num(/1) 176 : * ipr=1 177 : do 574 nbn=1,num(/1) 178 : is=icpr1(num(nbn,nbe)) 179 : qtot(is) =melva1.velche(nbn,nbe)+qtot(is) 180 : qrest(is)=qrest(is)+ 181 : $ melva1.velche(nbn,nbe)*(1.d0-velche(nbn,nbe)) 182 : qs=qs + velche(nbn,nbe) 183 : 574 continue 184 : qmoy(nbe)=qs/nbnn 185 : 575 continue 186 : * write(6,*) ' qmoy' 187 : * write(6,*) ( qmoy(ip),ip=1,qmoy(/1)) 188 : 570 continue 189 : * boucle sur les noeuds concernés on y passe deux fois, la premiere 190 : * pour traiter les noeuds qui basculent completement 191 : do 600 ipc=1,mchpo1.ipchp(/1) 192 : msoupo=mchpo1.ipchp(ipc) 193 : segact msoupo 194 : meleme=igeoc 195 : segact meleme 196 : mpoval=ipoval 197 : segact mpoval 198 : do 602 noc=1,nocomp(/2) 199 : if( nocomp(noc).eq.'Q') go to 603 200 : 602 continue 201 : 603 continue 202 : * write(6,*)'num' 203 : * write(6,*) ( num(1,iop),iop=1,num(/2)) 204 : do 601 nel=1,num(/2) 205 : * on commence par comparer qreac avec qlat pour voir si tout est mangé 206 : qreac=vpocha(nel,noc) 207 : ia = num(1,nel) 208 : ib = icpr1(ia) 209 : * write(6,*)' nel,ia,ib,icpr.icpr1(ia)', 210 : * $ nel,ia,ib,icpr.icpr1(ia) 211 : if(icpr.icpr1(ia).EQ.0) go to 601 212 : qres= qrest(ib) 213 : qto=qtot(ib) 214 : ic=ipos(ib)+1 215 : id=ipos(ib+1) 216 : * write(6,*)'ia,ib,qreac,qres,qto', 217 : * $ ia,ib,qreac,qres,qto,icpr.icpr1(ia) 218 : * write(6,*) 'ifois' , ifois 219 : if(ifois.eq.2) go to 620 220 : if( qreac.lt.0.D0) then 221 : if(-qreac.lt.qres*umxpre) go to 601 222 : xval=1.d0 223 : else 224 : if(qreac .lt.(qto-qres)*umxpre) go to 601 225 : xval=0.d0 226 : endif 227 : * tous les elements sont manges 228 : ipt2= imache(1) 229 : * write(6,987) ((ipt2.num(iu,ip),iu=1,4),ip=1,92) 230 : * 987 format(8I8) 231 : icpr.icpr1(ia)=0 232 : * write(6,*) 'traitement complet de ia ib', ia,ib 233 : do 610 ie=ic,id 234 : il=icon(ie) 235 : ilv=mod(il,ngrand) 236 : ieme= (il - ilv)/ngrand + 1 237 : ipt1= imache(ieme) 238 : mchaml=ichaml(ieme) 239 : melval=ielval(1) 240 : do 611 ino=1,ipt1.num(/1) 241 : if( ipt1.num(ino,ilv).eq.ia) then 242 : * write(6,*)'ie,il,ilv,ieme,mchaml,melval',ie,il,ilv, 243 : * $ ieme,mchaml,melval 244 : * write(6,*) ' ia,ie,ipt1.num(ino,ilv)',ia,ie,ipt1.num(ino,ilv) 245 : velche(ino,ilv)=xval 246 : go to 610 247 : endif 248 : 611 continue 249 : 610 continue 250 : go to 601 251 : 620 continue 252 : * il va falloir en tartiner sur certains elements 253 : * on prend d'abord l'élement dont la plus gande proportion est dans le 254 : * sens qui nous interesse 255 : * si ifon =1 on fonds sinon on solidifie 256 : do 631 ivo=ic,id 257 : ivoi(ivo-ic+1)=0 258 : 631 continue 259 : 630 continue 260 : ifon=0 261 : iiem=-1 262 : qma= -1.5 263 : qmi=1.5 264 : if(qreac.lt.0.D0) then 265 : ifon=1 266 : qreac=-qreac 267 : endif 268 : * write(6,*) ' ifon ' , ifon 269 : do 502 ie=ic,id 270 : if(ivoi(ie-ic+1).eq.1) go to 502 271 : il=icon(ie) 272 : ilv = mod(il,ngrand) 273 : ieme= (il - ilv)/ngrand + 1 274 : ipt1= imache(ieme) 275 : qmoye=imoy(ieme) 276 : * write(6,*) ' qmoy(ilv) ilv ieme ',ilv, qmoy(ilv) 277 : if( ifon.eq.1) then 278 : if(qmoy(ilv).lt.umxpre. and. qmoy(ilv).GT.qma) then 279 : iiem=ieme 280 : ieg=ie 281 : iilv=ilv 282 : qma=qmoy(ilv) 283 : endif 284 : else 285 : if(qmoy(ilv).ge.xpre.and.qmoy(ilv).lt.qmi) then 286 : iiem=ieme 287 : iilv=ilv 288 : ieg=ie 289 : qmi=qmoy(ilv) 290 : endif 291 : endif 292 : 502 continue 293 : * write(6,*) ' qma qmi iiem iilv', qma, qmi,iiem,iilv 294 : if ( iiem.eq.-1) go to 601 295 : ivoi(ieg-ic+1)=1 296 : * mise a jour du point en question dans l'element trouvé 297 : mchaml=ichaml(iiem) 298 : mcham1=mchel1.ichaml(iiem) 299 : do 506 n2=1,mcham1.nomche(/2) 300 : if(mcham1.nomche(n2).eq.'Q') GO TO 507 301 : 506 continue 302 : 507 continue 303 : melva1=mcham1.ielval(n2) 304 : melval=mchaml.ielval(1) 305 : ipt1=imache(iiem) 306 : do 505 ino=1,ipt1.num(/1) 307 : ir= ipt1.num(ino,iilv) 308 : * write(6,*) ' ir, ia ifon', ir,ia,ifon 309 : if( ir.eq.ia) then 310 : if(ifon.eq.1) then 311 : qmax=melva1.velche(ino,iilv)*(1.d0 - velche(ino,iilv)) 312 : * write(6,*) ' qmax qreac melva1. ',qmax,qreac, 313 : * $ melva1.velche(ino,iilv) 314 : if(qmax.le.qreac) then 315 : velche(ino,iilv)=1.d0 316 : if( (qreac-qmax)/melva1.velche(ino,iilv).le.xpre) 317 : $ go to 601 318 : qreac=qreac-qmax 319 : qreac = - qreac 320 : * write(6,*) ' on repart en 630 qreac ', qreac 321 : go to 630 322 : else 323 : xpro= qreac / melva1.velche(ino,iilv) 324 : velche(ino,iilv)=velche(ino,iilv) + xpro 325 : endif 326 : else 327 : qmin=melva1.velche(ino,iilv)*velche(ino,iilv) 328 : if(qmin.le.qreac) then 329 : velche(ino,iilv)=0.d0 330 : if(abs((qmin-qreac)/melva1.velche(ino,iilv)).le.xpre) 331 : $ go to 601 332 : qreac = qreac - qmin 333 : go to 630 334 : else 335 : xpro = qreac /melva1.velche(ino,iilv) 336 : velche(ino,iilv)=velche(ino,iilv) - xpro 337 : endif 338 : endif 339 : go to 601 340 : endif 341 : 505 continue 342 : 601 continue 343 : 600 continue 344 : if(ifois.eq.1) go to 1000 345 : * travail terminé il faut desactiver 346 : segsup icon,ipos,ipo2,icp,ivoi,icpr,imoy,qrest,qtot,qmoye 347 : * desactivation mchelm,mchel1 348 : do 961 ich=1,mchelm.imache(/1) 349 : meleme=imache(ich) 350 : segdes meleme 351 : mchaml=ichaml(ich) 352 : do 962 ive=1,ielval(/1) 353 : melval=ielval(ive) 354 : * test de validité des valeur 0.d0 <= val <= 1.d0 355 : do 976 iu=1,velche(/2) 356 : do 976 io=1,velche(/1) 357 : if( velche(Io,iu). lT.0.D0) then 358 : * write(6,*) ' phapro: proportion < 0. : ' ,velche(Io,iu) 359 : velche(Io,iu)=0.d0 360 : elseif (velche(Io,iu). gT.1.D0) then 361 : * write(6,*) ' phapro: proportion > 1. : ',velche(Io,iu) 362 : velche(Io,iu)=1.d0 363 : endif 364 : 976 continue 365 : segdes melval 366 : 962 continue 367 : segdes mchaml 368 : 961 continue 369 : segdes mchelm 369 : call ecrobj ('MCHAML',mchelm) 370 : mchelm=mchel1 371 : do 963 ich=1,mchelm.imache(/1) 372 : meleme=imache(ich) 373 : segdes meleme 374 : mchaml=ichaml(ich) 375 : do 964 ive=1,ielval(/1) 376 : melval=ielval(ive) 377 : segdes melval 378 : 964 continue 379 : segdes mchaml 380 : 963 continue 381 : do 965 ipc=1,mchpo1.ipchp(/1) 382 : msoupo = mchpo1.ipchp(ipc) 383 : meleme=igeoc 384 : segdes meleme 385 : mpoval=ipoval 386 : segdes mpoval 387 : segdes msoupo 388 : 965 continue 389 : segdes mchpo1 390 : return 391 : end 392 : 393 : 394 : 395 : 396 : 397 : 398 : 399 : 400 : 401 : 402 : 403 : 404 : 405 : 406 : 407 : 408 : 409 : 410 :
© Cast3M 2003 - All rights reserved.
Disclaimer