1 : C REDUAF SOURCE FANDEUR 13/01/30 21:15:00 7686 2 : subroutine reduaf(jchelm,mmodtm,mchel2,istri,iret,KERRE) 3 : implicit real*8(a-h,o-z) 4 : implicit integer (I-N) 5 : * on essaye de faire la reduction du champ par elemnt jchelm 6 : * sur le modele mmodel 7 : * resultat mchel2 iret=1 ok iret=0 rate 8 : *² INFCHE(.,6): LE SMINTE POINTE SUR UN SEGMENT D'INTEGRATION : 9 : * =1 AUX NOEUDS 10 : * =2 AU CENTRE DE GRAVITE 11 : * =3 AUX POINTS DE GAUSS POUR LA RIGIDITE 12 : * =4 AUX POINTS DE GAUSS POUR LA MASSE 13 : * =5 AUX POINTS DE GAUSS POUR LES CONTRAINTES 14 : * 15 : -INC CCOPTIO 16 : -INC SMCHAML 17 : -INC SMMODEL 18 : -INC SMCOORD 19 : c*-INC SMINTE 20 : -INC SMELEME 21 : SEGMENT INFO 22 : INTEGER INFELL(JG) 23 : ENDSEGMENT 24 : segment travg 25 : integer icpr(xcoor(/1)/(idim+1)) 26 : integer ityp (nbele) 27 : integer imin(nbele) 28 : integer ipart(nbele) 29 : integer ipos(nbele) 30 : endsegment 31 : segment inum(ino+1) 32 : segment ielno(ima) 33 : segment izone(nz) 34 : segment infc(nz) 35 : segment locglo(nmacom,nz) 36 : segment nokkom 37 : character*(8) nokomp(nbco) 38 : integer maxptv(nbco),maxpti(nbco) 39 : endsegment 40 : segment idimma 41 : integer imaxdi(n1) 42 : endsegment 43 : 44 : iret=0 45 : KERRE = 0 46 : 47 : mchelm=jchelm 48 : mmodel=mmodtm 49 : segact mchelm,mmodel 50 : segini,mchel2=mchelm 51 : nz=imache(/1) 52 : 53 : * Cas tres particulier : 54 : IF (nz.EQ.0) THEN 55 : * write(ioimp,*) 'CAS PARTICULIER nz = ',nz 56 : segdes,mchelm,mmodel 57 : * segdes,mchel2 58 : iret = 1 59 : RETURN 60 : ENDIF 61 : 62 : if (istri.eq.1) goto 1954 63 : * test formulation melange: on ajoute les modeles de phases 64 : nvim=0 65 : do i=1,kmodel(/1) 66 : imodel=kmodel(i) 67 : segact imodel 68 : if (formod(1).eq.'MELANGE') nvim=nvim+ivamod(/1) 69 : enddo 70 : if (nvim.ne.0) then 71 : nc=kmodel(/1) 72 : n1=nc+nvim 73 : segini mmode1 74 : do i=1,kmodel(/1) 75 : mmode1.kmodel(i)=kmodel(i) 76 : enddo 77 : do 600 i=1,kmodel(/1) 78 : imodel=kmodel(i) 79 : if (formod(1).eq.'MELANGE') then 80 : do 601 j=1,ivamod(/1) 80 : if (tymode(j).ne.'IMODEL') call erreur(21) 81 : nc=nc+1 82 : mmode1.kmodel(nc)=ivamod(j) 83 : 601 continue 84 : endif 85 : 600 continue 86 : mmodel=mmode1 87 : endif 88 : 89 : 1954 continue 90 : * 91 : * test rapide sur les pointeurs des maillages on prend tout ce que l'on peut 92 : * 93 : n133=0 94 : ia=0 95 : segini izone 96 : do 500 i=1,kmodel(/1) 97 : imodel=kmodel(i) 98 : segact imodel 99 : * write (6,*) ' reduad model i imodel nefmod conmod ', 100 : * > i,imodel,nefmod,conmod 101 : if(nefmod.eq.22) then 102 : segdes imodel 103 : go to 500 104 : endif 105 : do 501 j=1,nz 106 : * on parcourt tous les NZ chamelem elementaires 107 : * write (6,*) ' reduaf j imamod imache conmod conche ', 108 : * > j,imamod,imache(j),conmod,conche(j) 109 : if( imamod.eq.imache(j)) then 110 : * correspondance du maillage chamelem j avec maillage modele i 111 : if( conche(j).eq.conmod.or.conche(j).eq.' ') then 112 : ia=ia+1 113 : if (ia.gt.mchel2.imache(/1)) then 114 : n1=ia+10 115 : n133=ia 116 : l1=titche(/1) 117 : n3=infche(/2) 118 : segadj mchel2 119 : endif 120 : * write (6,*) ' reduaf boucle 501 ia j conche ',ia,j,conche(j) 121 : mchel2.conche(ia)=conche(j) 122 : mchel2.imache(ia)=imamod 123 : mcham3=ichaml(j) 124 : segini,mcham2=mcham3 125 : *? do 505 k=1,mcham2.ielval(/1) 126 : *? melva3=mcham2.ielval(k) 127 : *?** segini,melva2=melva3 128 : *? melva2=melva3 129 : *? mcham2.ielval(k)=melva2 130 : *? segdes melva2,melva3 131 : *? 505 continue 132 : mchel2.ichaml(ia)=mcham2 133 : *? segdes mcham2,mcham3 134 : segdes,mcham2 135 : do 504 k=1,infche(/2) 136 : mchel2.infche(ia,k)=infche(j,k) 137 : 504 continue 138 : izone(j)=1 139 : endif 140 : endif 141 : 501 continue 142 : 500 continue 143 : * les a-t-on tous pris? si oui, le resultat est le chamelem d'entree 144 : *** if(ia.eq.nz) then 145 : *** iret = 1 146 : * desactivation des segments 147 : *** do 509 i=1,kmodel(/1) 148 : *** imodel=kmodel(i) 149 : *** segdes imodel 150 : **509 continue 151 : *** segdes mchelm,mmodel,mchel2 152 : *** segsup izone 153 : *** return 154 : *** else 155 : * sinon, le resultat est un chamelem partiel 156 : if(ia.ne.0) then 157 : n1=ia 158 : n133=ia 159 : l1=titche(/1) 160 : n3=infche(/2) 161 : segadj mchel2 162 : mchel3=mchel2 163 : segini,mchel2=mchelm 164 : ib=0 165 : do 506 j=1,nz 166 : if( izone(j).eq.1) go to 506 167 : ib=ib+1 168 : * write (6,*) ' reduaf boucle 506 ib j conche ',ib,j,conche(j) 169 : mchel2.conche(ib)=conche(j) 170 : mchel2.imache(ib)=imache(j) 171 : mchel2.ichaml(ib)=ichaml(j) 172 : do 507 k=1,infche(/2) 173 : mchel2.infche(ib,k)=infche(j,k) 174 : 507 continue 175 : 506 continue 176 : n1=ib 177 : 178 : l1=titche(/1) 179 : n3=infche(/2) 180 : segadj mchel2 181 : segdes mchelm 182 : mchelm=mchel2 183 : mchel2=0 184 : else 185 : mchel3=0 186 : endif 187 : *** endif 188 : * on doit tester element par element 189 : * travail sur les maillages support du mchaml 190 : * 191 : segsup izone 192 : * 193 : * le nouveau chmelem ne peut avoir plus que kmodel(/1) * 194 : *(nbinfele+conmod+nophas) 195 : *nbinfele est le nombre de infele different dans le chamelem a reduire 196 : * on commence par compter nbinfele 197 : nz=imache(/1) 198 : segini infc 199 : izo=0 200 : n3=infche(/2) 201 : do 200 i=1,nz 202 : if(infc(i).ne.0) go to 200 203 : izo=izo+1 204 : infc(i)=izo 205 : do 201 j=i+1,nz 206 : do 202 k=1,n3 207 : if(k.eq.4) go to 202 208 : if(infche(i,k).ne.infche(j,k)) go to 201 209 : 202 continue 210 : * la zone j a le meme infche que la zone i 211 : * on teste aussi si meme constituant 212 : if(conche(i).ne.conche(j))go to 201 213 : infc(j)=izo 214 : 201 continue 215 : 200 continue 216 : * 217 : * on va dire que chaque sous zone du modele a izo sous partie 218 : * donc nouveau chamelem a izo*kmodel(/1) sous partie : creation. 219 : * Au debut chaque sous mchel contient toutes les composantes 220 : n1=izo * kmodel(/1) 221 : segini,idimma 222 : l1=titche(/1) 223 : segini mchel2 224 : mchel2.titche=titche 225 : mchel2.ifoche=ifoche 226 : do 205 i=1,kmodel(/1) 227 : imodel=kmodel(i) 228 : segact imodel 229 : MELE1=NEFMOD 230 : if( nefmod.eq.22) go to 205 231 : lesupp=6 232 : if(infele(6).ne.0) then 233 : infe6 = infele(6) 234 : else 234 : CALL ELQUOI(MELE1,0,lesupp,IPINF,IMODEL) 235 : info=ipinf 236 : segact info 237 : infe6 = infell(6) 238 : segsup info 239 : endif 240 : * write(6,*) 'infell',( infell( iou),iou=1,infell(/1)) 241 : na=(i-1)*izo 242 : do 206 j=1,izo 243 : *** imaxdi(na+j)=max(infell(6),infell(2)) 244 : * imaxdi(na+j)= infell(6) 245 : imaxdi(na+j)= infe6 246 : mchel2.imache(na+j)=imamod 247 : mchel2.conche(na+j)=conmod 248 : * write (6,*) ' reduad boucle 206 j na conmod ',j,na,conmod 249 : 206 continue 250 : * segsup info 251 : 205 continue 252 : 253 : *on cree locglo(j,i)=k qui dit que la j ieme inconnue 254 : * de la i ieme zone du mchelm est en k ieme position globale 255 : * d'abord activation du mchelm et max de n2=nomche(/2)(=ielval(/1)) 256 : nz=imache(/1) 257 : nmacom=0 258 : do 150 i=1,nz 259 : ipt1=imache(i) 260 : segact ipt1 261 : mchaml=ichaml(i) 262 : segact mchaml 263 : nmacom=max(nmacom,ielval(/1)) 264 : do 151 j=1,ielval(/1) 265 : melval=ielval(j) 266 : * write(6,*) 'reduaf activation de melval ith',ith,melval 267 : segact melval 268 : 151 continue 269 : 150 continue 270 : nbco=20 271 : nbcv=0 272 : segini nokkom,locglo 273 : do 152 i=1,nz 274 : mchaml=ichaml(i) 275 : do 153 j=1,nomche(/2) 276 : do 154 k=1,nbcv 277 : if(nokomp(k).eq.nomche(j)) go to 155 278 : 154 continue 279 : nbcv=nbcv+1 280 : if(nbcv.gt.nbco) then 281 : nbco=nbco+20 282 : segadj nokkom 283 : endif 284 : nokomp(nbcv)=nomche(j) 285 : k=nbcv 286 : 155 continue 287 : melval=ielval(j) 288 : if( typche(j).eq.'REAL*8') then 289 : * write(6,*) 'i j k nokomp maxpt velche ', 290 : * $ i,j, k , nokomp(k),maxptv(k),velche(/1) 291 : maxptv(k)=max(maxptv(k),velche(/1)) 292 : else 293 : maxpti(k)=max(maxpti(k),ielche(/1)) 294 : endif 295 : locglo(j,i)=k 296 : 153 continue 297 : 152 continue 298 : nbco=nbcv 299 : segadj nokkom 300 : n2=nbco 301 : * 302 : * il reste a prendre chaque sous chamelem est a le dispachtcher 303 : * dans les differents bouts. Pour faire des tests rapides on cherche pour 304 : * chaque element dans les sous zones du model quel est le point de 305 : * plus petit numero et on construit un tableau qui dit que les elements 306 : * ayant un certain point mini i sont dans les zones j 307 : * 308 : nbsous=kmodel(/1) 309 : nbele=0 310 : * on compte le nombre d'elements 311 : do 1 i=1,nbsous 312 : imodel=kmodel(i) 313 : if(nefmod.eq.22) go to 1 314 : meleme=imamod 315 : segact meleme 316 : nbele= num(/2)+ nbele 317 : 1 continue 318 : segact mcoord 319 : segini travg 320 : 321 : * fabrication d'une numerotation locale icpr 322 : * fabrication de imin(k) qui donne le plus petit noeud de l'element 323 : * fabrication de ipart (k) qui dit dans quel sous paquet est l'element 324 : * fabrication de ityp qui donne le type de l element 325 : * fabrique de ipos(k): numero d'ordre de l'element dans le sous paquet 326 : idej=0 327 : ino=0 328 : do 2 i=1,nbsous 329 : imodel = kmodel(i) 330 : if( nefmod.eq.22) go to 2 331 : meleme=imamod 332 : do 3 j=1,num(/2) 333 : ia = num(1,j) 334 : do 4 m=1,num(/1) 335 : ib=num(m,j) 336 : ia = min ( ia,ib) 337 : if(icpr(ib).eq.0) then 338 : ino=ino+1 339 : icpr(ib)=ino 340 : endif 341 : 4 continue 342 : jidej= j + idej 343 : imin(jidej)=ia 344 : ipart(jidej)=i 345 : ityp(jidej)=itypel 346 : ipos(jidej)=j 347 : 3 continue 348 : idej = idej + num(/2) 349 : 2 continue 350 : segini inum 351 : * on va compter combien (au max)d'elements touchent un noeud mini 352 : do 5 i=1,nbele 353 : ia=imin(i) 354 : ib = icpr(ia) 355 : inum(ib)=inum(ib)+1 356 : 5 continue 357 : * pour chaque noeud, on stocke la fin de rangement dans ielno 358 : do i=1,ino 359 : inum(i+1)=inum(i+1)+inum(i) 360 : enddo 361 : ima=inum(ino+1) 362 : 363 : * write (6,*) ' reduaf ima',ima 364 : segini ielno 365 : * dans le tableau ielno on dira quels sont les elements ayant 366 : * le ino iem noeud local comme plus petit noeud 367 : do 7 i=1,nbele 368 : ia = imin(i) 369 : iloc=icpr(ia) 370 : ib=inum(iloc) 371 : inum(iloc)=inum(iloc)-1 372 : ielno(ib)=i 373 : 7 continue 374 : * on a maintenant dans inum le debut de rangement dans ielno 375 : * 376 : * on travaille sous chaml par sou-chaml 377 : * 378 : do 8 i=1,imache(/1) 379 : ipt1=imache(i) 380 : segact ipt1 381 : mchaml=ichaml(i) 382 : segact mchaml 383 : do 218 j=1,ielval(/1) 384 : melval=ielval(j) 385 : segact melval 386 : 218 continue 387 : ityref=ipt1.itypel 388 : c* indinf=infc(i) 389 : c* iconce= (i-1)*izo 390 : do 9 j=1,ipt1.num(/2) 391 : imi=ipt1.num(1,j) 392 : c* ina=0 393 : do 10 k=2,ipt1.num(/1) 394 : imi= min(imi,ipt1.num(k,j)) 395 : 10 continue 396 : iloc = icpr(imi) 397 : if( iloc.eq.0) go to 9 398 : do 11 k=inum(iloc)+1,inum(iloc+1) 399 : ie = ielno(k) 400 : if(ie.eq.0) write (6,*) ' remplissage incorrect de ielno ' 401 : * l element ie est un bon candidat on teste le itypel 402 : * puis le constituant 403 : if(ityref.ne.ityp(ie)) go to 11 404 : nbs= ipart(ie) 405 : imodel=kmodel(nbs) 406 : meleme=imamod 407 : if(nefmod.eq.22) go to 11 408 : iloo=ipos(ie) 409 : if(conche(i).ne.conmod.and.conche(i).ne.' ' ) go to 11 410 : * on regarde si chaque noeud de j appartient aussi a ie 411 : do 12 m=1,ipt1.num(/1) 412 : if(num(m,iloo).ne.ipt1.num(m,j)) go to 11 413 : 12 continue 414 : * on a trouve une correspondance on peut dispatcher les valeurs 415 : ipla= ipos(ie) 416 : ires=(nbs-1)*izo+infc(i) 417 : if( mchel2.ichaml(ires) .eq.0) then 418 : n2=nbco 419 : segini mcham2 420 : mchel2.ichaml(ires)=mcham2 421 : do 21 m=1,n3 422 : mchel2.infche(ires,m)=infche(i,m) 423 : 21 continue 424 : else 425 : mcham2=mchel2.ichaml(ires) 426 : endif 427 : do 16 m=1,ielval(/1) 428 : inoco=locglo(m,i) 429 : melval=ielval(m) 430 : if(mcham2.ielval(inoco).eq.0) then 431 : * nb de point support du chamelem 432 : iou=infche(i,6) 433 : if( iou.eq.1) then 434 : n2ptel=ipt1.num(/1) 435 : n1ptel=n2ptel 436 : elseif(iou.eq.2) then 437 : n2ptel=1 438 : n1ptel=1 439 : else 440 : n1ptel=imaxdi(ires) 441 : n2ptel=imaxdi(ires) 442 : endif 443 : ipt3=imamod 444 : * n1ptel=velche(/1) 445 : * n1ptel=imaxdi(ires) 446 : * n2ptel=ielche(/1) 447 : * n2ptel=imaxdi(ires) 448 : * write(6,*) ' n1ptel ', n1ptel 449 : if(typche(m).eq.'REAL*8') then 450 : n1el=ipt3.num(/2) 451 : n2el=0 452 : n2ptel=0 453 : else 454 : n2el=ipt3.num(/2) 455 : n1el=0 456 : n1ptel=0 457 : endif 458 : segini melva2 459 : mcham2.ielval(inoco)=melva2 460 : mcham2.typche(inoco)=typche(m) 461 : mcham2.nomche(inoco)=nokomp(inoco) 462 : else 463 : * verification du type 464 : if( mcham2.typche(inoco).ne.typche(m)) then 465 : * incompatibilite du type de composante entre champs 466 : KERRE = 472 467 : endif 468 : melva2= mcham2.ielval(inoco) 469 : endif 470 : * write(6,*) ' typche(m)' ,typche(m) 471 : if(typche(m).eq.'REAL*8' ) then 472 : ** if(ipla.gt.melva2.velche(/2) ) then 473 : * write(6,*) ' erreur ipla dim2',ipla,melva2.velche(/2) 474 : * endif 475 : * if(velche(/1) . gt.melva2.velche(/1)) then 476 : * write(6,*) 'bizarre ***********',velche(/1),melva2.velche(/1) 477 : * endif 478 : ibj=min (velche(/2),j) 479 : * write(6,*) 'melva2.velche(/1)',melva2.velche(/1) 480 : do 17 n=1,melva2.velche(/1) 481 : ibn=min(velche(/1),n) 482 : melva2.velche(n,ipla)=velche(ibn,ibj) 483 : 17 continue 484 : else 485 : 486 : ibj=min (ielche(/2),j) 487 : do 18 n=1,melva2.ielche(/1) 488 : ibn=min(ielche(/1),n) 489 : melva2.ielche(n,ipla)=ielche(ibn,ibj) 490 : 18 continue 491 : endif 492 : 16 continue 493 : 11 continue 494 : 9 continue 495 : 8 continue 496 : * 497 : * suppression des segments de trvail 498 : * 499 : segsup locglo,infc,ielno 500 : segsup inum 501 : segsup travg 502 : segsup inum 503 : segsup nokkom 504 : * 505 : * compactage du resultat 506 : * 507 : n1=0 508 : do 30 i=1,mchel2.imache(/1) 509 : if( mchel2.ichaml(i).ne.0) then 510 : n1 = n1+1 511 : * write (6,*) ' reduad boucle 30 i n1 conche ',i,n1, 512 : * > mchel2.conche(i) 513 : mchel2.conche(n1)=mchel2.conche(i) 514 : mchel2.ichaml(n1)=mchel2.ichaml(i) 515 : mchel2.imache(n1)=mchel2.imache(i) 516 : do 31 j=1,n3 517 : mchel2.infche(n1,j)=mchel2.infche(i,j) 518 : 31 continue 519 : mcham2=mchel2.ichaml(i) 520 : * on compacte les composantes 521 : n2=0 522 : do 32 j=1,mcham2.ielval(/1) 523 : if( mcham2.ielval(j).ne.0) then 524 : n2=n2+1 525 : mcham2.ielval(n2)=mcham2.ielval(j) 526 : mcham2.typche(n2)=mcham2.typche(j) 527 : mcham2.nomche(n2)=mcham2.nomche(j) 528 : melva2=mcham2.ielval(j) 529 : segdes melva2 530 : endif 531 : 32 continue 532 : segadj mcham2 533 : segdes mcham2 534 : endif 535 : 30 continue 536 : * 537 : * desactivation des entrees :modele et chamelem 538 : * 539 : do 33 i=1,kmodel(/1) 540 : imodel=kmodel(i) 541 : ipt1=imamod 542 : segdes ipt1 543 : segdes imodel 544 : 33 continue 545 : segdes mmodel 546 : do 34 i=1,imache (/1) 547 : ipt1=imache(i) 548 : segdes ipt1 549 : mchaml=ichaml(i) 550 : do 35 j=1,ielval(/1) 551 : melval=ielval(j) 552 : *pv segdes melval 553 : 35 continue 554 : segdes mchaml 555 : 34 continue 556 : segdes mchelm 557 : * 558 : * ajout de mchel3 559 : * 560 : if( n1.eq.0) then 561 : segsup mchel2 562 : if(mchel3.eq.0) then 563 : * write(6,*) ' on tombe ici' 564 : KERRE = 21 565 : return 566 : else 567 : mchel2=mchel3 568 : go to 1111 569 : endif 570 : endif 571 : n1dep=n1 572 : n1=n1+n133 573 : segadj mchel2 574 : if(mchel3.ne.0) then 575 : segact mchel3 576 : do 40 j=1,n133 577 : mchel2.ichaml(n1dep+j)=mchel3.ichaml(j) 578 : mchel2.imache(n1dep+j)=mchel3.imache(j) 579 : mchel2.conche(n1dep+j)=mchel3.conche(j) 580 : * write (6,*) ' reduaf boucle 40 j n1dep conch', 581 : * > j,n1dep,mchel3.conche(j) 582 : do 41 k=1,mchel3.infche(/2) 583 : mchel2.infche(n1dep+j,k)=mchel3.infche(j,k) 584 : 41 continue 585 : 40 continue 586 : segsup mchel3 587 : endif 588 : segdes mchel2 589 : * 590 : * sortie 591 : * 592 : 1111 iret=1 593 : * 594 : * on compacte le resultat 595 : * 596 : segact mchel2 597 : do i=1,mchel2.ichaml(/1) 598 : mchaml=mchel2.ichaml(i) 599 : segact mchaml*mod 600 : * write (6,*) ' reduaf boucle finale i conche ', 601 : * > i,mchel2.conche(i) 602 : do j=1,ielval(/1) 603 : ielin=ielval(j) 603 : call comred(ielin) 604 : c* ielval(j)=ielin 605 : *pv melva1=ielin 606 : *pv segdes melva1 607 : enddo 608 : segdes mchaml 609 : enddo 610 : *pv segdes mchel2 611 : segsup idimma 612 : 613 : return 614 : end 615 : 616 : 617 :
© Cast3M 2003 - All rights reserved.
Disclaimer