1 : C PRO2 SOURCE BP208322 13/09/02 21:15:00 7816 2 : SUBROUTINE PRO2(IPGEA,IPCHEL,isort,IPOUT,ilphmo) 3 : C===================================================================== 4 : C 5 : C SOUS-PROGRAMME DE L'OPERATEUR PROIET 6 : C 7 : C ENTREES : 8 : C 9 : C IPGEA = POINTEUR DE UN OBJET MAILLAGE. 10 : C IPCHEL = POINTEUR DE UN OBJET CHAMELEM 11 : C isort 0 on veut un chpo brutal a 1 sous zone non repartitionné 12 : C en sortie pointuer su un segment ipca 13 : C 1 on veut un champ conforme partitionne 14 : c ilphmo pointeur sur un listmot phases du modele cible 15 : C SORTIES : 16 : C 17 : C IPOUT = POINTEUR SUR UN OBJET DE TYPE CHPOINT 18 : C 19 : C===================================================================== 20 : C 21 : C TRES INSPIRE DE LA FONCTION ACCRO DEVELOPPEE PAR JMB 22 : C 23 : C===================================================================== 24 : 25 : IMPLICIT INTEGER(I-N) 26 : IMPLICIT REAL*8 (A-H,O-Z) 27 : *- 28 : -INC CCREEL 29 : -INC CCOPTIO 30 : -INC SMCOORD 31 : -INC SMELEME 32 : -INC SMCHPOI 33 : -INC SMLCHPO 34 : -INC SMLMOTS 35 : -INC SMCHAML 36 : -INC TMTRAV 37 : -INC CCGEOME 38 : -INC SMMODEL 39 : -INC SMLENTI 40 : SEGMENT nomiii 41 : character*4 NOMINC(0) 42 : endsegment 43 : DATA EPSI,GRAND/1.D-9,1.D20/ 44 : SEGMENT MDATA 45 : INTEGER IPCA(NBSZ) 46 : ENDSEGMENT 47 : SEGMENT IPEX(NODES) 48 : SEGMENT ISOM(NBSOM(IPT3.ITYPEL)) 49 : SEGMENT MTR1 50 : CHARACTER*4 IPCOM(0) 51 : ENDSEGMENT 52 : SEGMENT SEPHAS 53 : CHARACTER*8 LIPHAS(NPHAS) 54 : ENDSEGMENT 55 : SEGMENT LLIZA(NLIZA6) 56 : SEGMENT/MTR4/(IPHAR(0)) 57 : segment inomil(0) 58 : SEGMENT/MTRA/(NOPOIN(XCOOR(/1)/(IDIM+1))) 59 : * 60 : LOGICAL drphas 61 : C 62 : segment info 63 : integer infell(JG) 64 : endsegment 65 : C 66 : SEGMENT ICOUNT 67 : INTEGER IEINT(2,NODES),IDEJVU(NODES),IACCRO 68 : * ient (1,i)= j, ient(2,i)=k veut dire que le ieme point de ipt2 est dans 69 : * le k ieme element de la jeme sous zone. 70 : * idejvu(i)=1 veut dire que le ieme noeud de ipt2 est deja traité 71 : REAL*8 QQQ(3,NODES),CKRIT(NODES) 72 : * CKRIT est le meilleur des criteres déja rencontré 73 : ENDSEGMENT 74 : SEGMENT ICOSOM 75 : INTEGER IPTT(NBSZ),ICENV(NBSZ),IVOL(NBSZ) 76 : ENDSEGMENT 77 : SEGMENT WORK1 78 : REAL*8 CORDEL(3,NBNN),QSI(3),XPO(3),SHPP(6,NBNN),V1(3),V2(3) 79 : REAL*8 XEL(3,NBNN),BPSS(3,3),XPT(3),cored(3,NBNN),V3(3) 80 : REAL*8 R(3),U(3),W(3,3) 81 : INTEGER II(3) 82 : ENDSEGMENT 83 : segment siele 84 : integer iiele(2,nsous) 85 : endsegment 86 : external shape 87 : 88 : C lecture optionnelle du modele associe au MCHAML origine 89 : IPMOD=0 90 : SIELE= 0 90 : call lirobj('MMODEL',IPMOD,0,iret) 91 : * write(6,*) ' dans pro2 lecture de ipmod ',ipmod 92 : if(ipmod.ne.0) then 93 : mmodel = ipmod 94 : segact mmodel 95 : NSOUS=KMODEL(/1) 96 : 97 : segini siele 98 : 99 : DO 30 ISOUS = 1, NSOUS 100 : IMODEL = KMODEL(ISOUS) 101 : SEGACT,IMODEL 102 : MELEME = IMAMOD 103 : SEGACT MELEME 104 : MELE = NEFMOD 105 : C COQUE INTEGREE OU PAS ? 106 : IF(INFMOD(/1).NE.0)THEN 107 : NPINT=INFMOD(1) 108 : ELSE 109 : NPINT=0 110 : ENDIF 111 : C 111 : call elquoi(mele,npint,6,ipinf,imodel) 112 : if (ierr.ne.0) then 113 : segdes imodel,mmodel 114 : C#MC segsup mchelm 115 : iret=0 116 : return 117 : endif 118 : C 119 : info=ipinf 120 : iiele(1,isous)=IMAMOD 121 : iiele(2,isous)= infelL(14) 122 : segdes meleme,imodel 123 : segsup info 124 : 30 CONTINUE 125 : endif 126 : C lecture optionnelle du critere de rattrapage de points exterieurs 127 : ccrit=1.D-5 127 : call lirree(ccrit,0,iretou) 128 : ccrit = -ABS(ccrit) 129 : * write(6,*) ' dans pro2 lecture de ccrit', ccrit 130 : C 131 : C ON VA CHANGER LE MAILLAGE RECEPTEUR EN ELEMENTS POI1 132 : C SI IL NE L EST PAS DEJA 133 : 134 : IPT2 = IPGEA 135 : SEGACT IPT2 135 : IF(IPT2.ITYPEL.NE.1) CALL CHANGE(IPGEA,1) 136 : SEGDES IPT2 137 : IF(IERR.NE.0) RETURN 138 : C 139 : IPT2 = IPGEA 140 : SEGACT IPT2 141 : NODES = IPT2.NUM(/2) 142 : SEGDES IPT2 143 : SEGINI ICOUNT 144 : IACCRO = 0 145 : do 3259 iu=1,nodes 146 : 3259 ckrit(iu)=-1.d10 147 : C 148 : C LE MAILLAGE SUPPORT DU CHAMELEM DOIT ETRE MASSIF 149 : C 150 : MCHEL1 = IPCHEL 151 : SEGACT MCHEL1 152 : NBSZ = MCHEL1.IMACHE(/1) 153 : * write(6,*) ' nbsz',nbsz 154 : C 155 : nbnnma=0 156 : DO 10 IOB=1, NBSZ 157 : C 158 : IPT1=MCHEL1.IMACHE(IOB) 159 : SEGACT IPT1 160 : C 161 : ITY = IPT1.ITYPEL 162 : C En dimension 1 : EF massifs s'appuient sur SEG2 ou SEG3 163 : IF (IDIM.EQ.1) THEN 164 : IF (ITY.NE.2.AND.ITY.NE.3) THEN 164 : CALL ERREUR(16) 165 : SEGDES IPT1 166 : RETURN 167 : ENDIF 168 : ELSE 169 : C En dimension 2,3 : 170 : cbp IF (ITY.LT.4) THEN 171 : IF (ITY.LT.4.OR.(IDIM.eq.3.and.ITY.LT.14)) THEN 172 : write(ioimp,*) 'IDIM,ITYPEL=',IDIM,ITY 172 : CALL ERREUR(16) 173 : SEGDES IPT1 174 : RETURN 175 : ENDIF 176 : ENDIF 177 : nbnnma=max(nbnnma,ipt1.num(/1)) 178 : C 179 : 10 CONTINUE 180 : C SEGACT MCOORD 181 : 182 : SEGINI MDATA 183 : Ckich on regarde les phases du MCHAML a projeter 184 : nphas = nbsz 185 : segini sephas 186 : liphas(1) = mchel1.conche(1)(17:24) 187 : nphas = 1 188 : if ( mchel1.conche(/2).gt.1) then 189 : do 11 jpha =2, mchel1.conche(/2) 190 : do k = 1,nphas 191 : if (mchel1.conche(jpha)(17:24).eq.liphas(k)) goto 11 192 : enddo 193 : nphas = nphas + 1 194 : liphas(nphas) = mchel1.conche(jpha)(17:24) 195 : 11 continue 196 : segadj sephas 197 : endif 198 : Ckich on fabrique un maillage contenant tous les maillages du champ 199 : c de pointeurs distincts 200 : nbsous=nbsz 201 : nbelem=0 202 : nbnn=0 203 : nbref=0 204 : segini ipt6 205 : 206 : nbsous = 1 207 : ipt6.lisous(1) = mchel1.imache(1) 208 : DO 20 IOB=2,NBSZ 209 : do is= 1, nbsous 210 : if (ipt6.lisous(is).eq.mchel1.imache(iob)) goto 20 211 : enddo 212 : c verifie que les maillages sont disjoints : avertit 213 : do 22 is= 1, nbsous 214 : ip1 = ipt6.lisous(is) 215 : ip2 = mchel1.imache(iob) 215 : call interb(ip1,ip2,iret,iob1) 216 : if (iret.eq.1) goto 22 217 : write(ioimp,*) 'maillages non disjoints dans le MCHAML' 218 : * kich : on peut faire mieux, pardon. 219 : goto 23 220 : 22 continue 221 : 222 : 23 nbsous = nbsous + 1 223 : IPT6.lisous(nbsous)=MCHEL1.IMACHE(IOB) 224 : 20 continue 225 : segadj ipt6 226 : C 227 : C ICI LE ZONAGE 228 : C============================================ 228 : CALL ZONEG2(IPT6,IPT2,ICount,1) 229 : c write(6,*) ' dans pro2 sortie de zoneg2' 230 : C============================================ 231 : C 232 : C------------------------- 233 : C TABLEAU des POINT PAS VUS 234 : C en previson de rectifications ulterieures 235 : NBOUB = NODES - IACCRO 236 : * write(6,*) ' nboubl ' , nboub 237 : SEGACT IPT2 238 : 239 : * SEGINI IOUBLI 240 : * SEGACT IEXCLU 241 : IUOO = 0 242 : DO 21 I=1,NODES 243 : IP = IPT2.NUM(1,I) 244 : IF(IDEJVU(I).EQ.1) THEN 245 : go to 210 246 : ELSE 247 : IF(CKRIT(I).GT.CCRIT) THEN 248 : IDEJVU(I)=1 249 : IUOO= IUOO + 1 250 : * IPTOU(IUO,1)=IP 251 : ELSE 252 : * IP = IPT2.NUM(1,I) 253 : * write(6,*) ' ip ckrit (i)' 254 : * write(6,*) ip, ckrit(i), qqq(1,i),qqq(2,i),qqq(3,i) 255 : GO TO 21 256 : ENDIF 257 : endif 258 : 210 IPCA(IEINT(1,I))=1 259 : 21 CONTINUE 260 : iaccro=iaccro+iuoo 261 : * write(6,*) ' nodes iaccro',nodes,iaccro 262 : CC on donne le nombre de points non projetés 263 : 264 : if (NBOUB.ne.iuoO) then 265 : INTERR(1)=NBOUB-IUOO 266 : * write(6,*) ' nombre de points non accrochés ',interr(1) 266 : CALL ERREUR(-322) 267 : endif 268 : 269 : segact ipt6 270 : Ckich : on fabrique 1 chpoint par phase presente dans le MMODEL 271 : c stocke dans un listchpo, dans l ordre de MOTS (posibilite terme nul) 272 : if (ilphmo.lt.0) then 273 : jgn = 8 274 : jgm = 1 275 : segini mlmots 276 : mots(1) = liphas(1) 277 : else 278 : mlmots = ilphmo 279 : c n1 = nphas 280 : n1 = mots(/2) 281 : segini mlchpo 282 : endif 283 : 284 : ipout = 0 285 : DO 700 IMP = 1,MOTS(/2) 286 : drphas = .false. 287 : do iphas = 1,NPHAS 288 : if(mots(imp).eq.liphas(iphas)) then 289 : drphas = .true. 290 : goto 703 291 : endif 292 : enddo 293 : 703 if(.not.drphas) then 294 : write(6,*) ' pas de projection phase ', mots(imp) 295 : goto 700 296 : endif 297 : 298 : C** recherche du noms de composantes existantes dans le chpoint 299 : segini nomiii 300 : NA=0 301 : DO 220 i=1,nbsz 302 : ckich maintenant maillage disjoints if(ipca(i).eq.0) go to 220 303 : if (mchel1.conche(i)(17:24).ne.liphas(iphas)) goto 220 304 : mchaml = mchel1.ichaml(i) 305 : segact mchaml 306 : do 221 j=1,nomche(/2) 307 : IF( TYPCHE(J).NE.'REAL*8') GO TO 221 308 : IF(NA.NE.0) THEN 309 : do 222 k=1,nominc(/2) 310 : if( nominc(k).eq.nomche(j)(1:4)) go to 221 311 : 222 continue 312 : nominc(**)=nomche(j) 313 : ELSE 314 : nominc(**)=nomche(j) 315 : na=1 316 : ENDIF 317 : 221 continue 318 : 220 continue 319 : nnin=nominc(/2) 320 : nnnoe=iaccro 321 : * write(6,*) ' nnin nnoe ',nnin,nnnoe,iphas 322 : segini mtrav 323 : C 324 : C**** remplissage de mtrav en calculant les valeurs 325 : C 326 : do 311 iu = 1 , nominc(/2) 327 : inco(iu)=nominc(iu) 328 : 311 continue 329 : segsup nomiii 330 : 331 : c recense les zones du mchel1 compatibles avec le maillage de ipt6 et la phase 332 : nliza6 = ipt6.lisous(/1) 333 : segini lliza 334 : do iza6 = 1,nliza6 335 : do ips = 1,nbsz 336 : IF (ipt6.lisous(iza6).eq.mchel1.imache(ips).AND. 337 : & liphas(iphas).eq.mchel1.conche(ips)(17:24)) THEN 338 : if (lliza(iza6).eq.0) then 339 : jg = 1 340 : segini mlenti 341 : LLIZA(IZA6) = mlenti 342 : lect(1) = ips 343 : else 344 : mlenti = lliza(iza6) 345 : jg = lect(/1)+1 346 : segadj mlenti 347 : lect(jg) = ips 348 : endif 349 : ENDIF 350 : enddo 351 : mlenti = lliza(iza6) 352 : enddo 353 : 354 : segact ipt2 355 : ipla=0 356 : nbnn = nbnnma 357 : SEGINI ,WORK1 358 : do 224 I=1,nodes 359 : if(idejvu(i).eq.1) then 360 : ipla=ipla+1 361 : igeo(ipla)=ipt2.num(1,i) 362 : ip=ipt2.num(1,i) 363 : igeo(ipla)=ipt2.num(1,i) 364 : iza6=ieint(1,i) 365 : iel=ieint(2,i) 366 : * 367 : IF (LLIZA(iza6).EQ.0) goto 302 368 : mlenti = lliza(iza6) 369 : nza = lect(/1) 370 : 371 : do 300 iz = 1,nza 372 : IZA=lect(iz) 373 : 374 : 301 mchaml=mchel1.ichaml(iza) 375 : ipt1=mchel1.imache(iza) 376 : segact ipt1, mchaml 377 : 378 : 379 : if(ipmod.ne.0) then 380 : C recherche element fini pour fonction de forme dans qsijs 381 : do lk=1,nsous 382 : if(ipt1.eq.iiele(1,lk)) then 383 : kel=iiele(2,lk) 384 : goto 31 385 : endif 386 : enddo 387 : 31 continue 388 : else 389 : C par defaut on utilise les fonctions de transformation geometriques 390 : kel =IPT1.ITYPEL 391 : endif 392 : C 393 : C Coordonnées des sommets de l'élements IEL 394 : C 395 : nbnn=ipt1.num(/1) 395 : CALL DOXE(XCOOR,IDIM,NBNN,IPT1.NUM,IEL,CORDEL) 396 : C 397 : IDI2 = IDIM 398 : do ih=1,idim 399 : qsi(ih)=qqq(ih,i) 400 : enddo 401 : C 402 : C Coordonnées du point dans l'element de référence 403 : C 403 : CALL SHAPE(QSI(1),QSI(2),QSI(3),KEL,SHPP,IRET) 404 : * if ( qsi(1).lt.-1.00001.or.qsi(1).gt.100001) then 405 : * if ( qsi(2).lt.-1.00001.or.qsi(2).gt.100001) then 406 : * if ( qsi(3).lt.-1.00001.or.qsi(3).gt.100001) then 407 : * write(6,*) ip,iza,iel,kel,(qsi(ih),ih=1,idim) 408 : * write(6,*) ( shpp(1,ioup),ioup=1,nbnn) 409 : * endif 410 : * endif 411 : * endif 412 : C 413 : C Calcul des valeurs pour chaque type de composantes 414 : C 415 : do 200 il=1,ielval(/1) 416 : IF (TYPCHE(il).EQ.'REAL*8') THEN 417 : C 418 : do 225 kcol=1,nnin 419 : * write(6,2323)il,nomche(il)(1:4),inco(1) 420 : * 2323 format('IL ', I2, 'nomche', A4, ' inco' , A4) 421 : if(nomche(il)(1:4).eq.inco(kcol)) go to 226 422 : 225 continue 423 : 226 continue 424 : * if(kcol.ne.1) then 425 : * write(6,*) 'kcol nnin',kcol,nnin 426 : * write(6,*)' iza il ipla',iza,il,ipla 427 : * return 428 : * endif 429 : MELVAL = IELVAL(il) 430 : SEGACT MELVAL 431 : IEMN = MIN(IEL,VELCHE(/2)) 432 : NHAR(kcol) = MCHEL1.INFCHE(iza,3) 433 : DO 150, NP = 1, NBNN 434 : IBMN = MIN(NP,velche(/1)) 435 : bb(kcol,ipla)=SHPP(1,NP)*VELCHE(IBMN,IEMN)+bb(kcol,ipla) 436 : C 437 : 150 CONTINUE 438 : ibin(kcol,ipla)=ibin(kcol,ipla)+1 439 : * SEGDES MELVAL 440 : C 441 : ENDIF 442 : 200 CONTINUE 443 : 300 continue 444 : 302 continue 445 : endif 446 : 224 CONTINUE 447 : 448 : do iza6 = 1,nliza6 449 : mlenti = lliza(iza6) 450 : segsup mlenti 451 : enddo 452 : segsup lliza 453 : 454 : ITRAV=MTRAV 455 : DO 4370 IB=1,NNNOE 456 : DO 4370 IA=1,NNIN 457 : * write(6,fmt='('' ia ib bb ibin '',2i4,e12.5,i5)')ia,ib,bb(ia,ib), 458 : * * ibin(ia,ib) 459 : IF (IBIN(IA,IB).GT.1) THEN 460 : * write(6,*) ' ibin gt 1 ib ,ia' , ib ,ia, ibin(ia,ib) 461 : BB(IA,IB)=BB(IA,IB) / IBIN(IA,IB) 462 : IBIN(IA,IB)=1 463 : ENDIF 464 : 4370 CONTINUE 465 : * SEGDES MTRAV 466 : C 467 : C RECONSTUCTION DE LA PARTITION 468 : C 468 : CALL CRECHP(ITRAV,ICHPOI1) 469 : C 470 : IPOUT=ICHPOI1 471 : MCHPOI=ICHPOI1 472 : SEGACT MCHPOI*MOD 473 : * le champ est forcement diffus 474 : JATTRI(1)=1 475 : MTYPOI='PROJECTI' 476 : IFOPOI=IFOUR 477 : SEGDES MCHPOI 478 : SEGSUP MTRAV 479 : ckich dans le listchpoint 480 : * write(6,*) 'pro2',ilphmo,mchpoi,ipout 481 : if (ilphmo.gt.0) then 482 : * write(6,*) 'pro2',mchpoi,ipout,'iphas',liphas(iphas),iphas,imp 483 : ichpoi(imp) = MCHPOI 484 : endif 485 : 700 CONTINUE 486 : if (ipout.le.0) then 486 : call erreur(-79) 487 : return 488 : endif 489 : if (ilphmo.gt.0) ipout = mlchpo 490 : SEGSUP ICOUNT,MDATA 491 : 492 : RETURN 493 : END 494 : 495 : 496 :
© Cast3M 2003 - All rights reserved.
Disclaimer