1 : C PROM SOURCE FD218221 14/01/31 21:15:12 7936 2 : 3 : SUBROUTINE PROM(IPMODE,IPCHE,IPCHT,ISUP,IPOUT) 4 : IMPLICIT INTEGER(I-N) 5 : IMPLICIT REAL*8 (A-H,O-Z) 6 : EXTERNAL SHAPE 7 : *-------------------------------------------------------------------- * 8 : * * 9 : * Sous-programme appelé par l'opérateur PROI * 10 : * Projection d'un chamelem sur le support d un modele par * 11 : * minimisation de l integrale (Ue-Us)**2 ds sur chaque element * 12 : * Ue evalue a partir des fonctions de formes du modele associe * 13 : * au champ d origine si il est fourni ( sinon a partir des * 14 : * fonctions de transformation geometriques de l element) * 15 : * * 16 : * L integrale est exprime à partir des fonctions de forme associees* 17 : * au modèle support recepteur * 18 : *---------------------------------------------------------------------* 19 : * * 20 : * ENTREE : * 21 : * IPMODE modele sur lequel on va projeter * 22 : * IPCHE MCHAML de caracteristique (si le maodele comprends * 23 : * des coques) * 24 : * IPCHT MCHAML a projeter (obligatoirement defini aux noeuds) 25 : * ISUP remis a 1 pour le moment (on travaille aux noeuds)* 26 : * * 27 : * SORTIE : * 28 : * IPOUT MCHAML resultat * 29 : * * 30 : *---------------------------------------------------------------------* 31 : * 32 : -INC CCOPTIO 33 : -INC SMCOORD 34 : -INC SMMODEL 35 : -INC SMCHAML 36 : -INC SMELEME 37 : -INC SMCHPOI 38 : -INC SMLCHPO 39 : -INC SMLMOTS 40 : -INC TMTRAV 41 : -INC SMINTE 42 : -INC CCREEL 43 : * 44 : PARAMETER ( O1=1.D0,O2=2.D0) 45 : character*4 nomm 46 : segment snomip 47 : character*4 nomip(0) 48 : endsegment 49 : * SEGMENT MDATA 50 : * INTEGER IPCA(*),IPEXI(*) 51 : * ENDSEGMENT 52 : segment ipcb(nsschp) 53 : segment icpr(2,xcoor(/1)/idim+1) 54 : segment icpr1(xcoor(/1)/idim+1) 55 : C 56 : segment info 57 : integer infell(JG) 58 : endsegment 59 : C 60 : segment sicoq 61 : integer icocoq(nbsous),imfr(nbsous),iiele(nbsous) 62 : integer iminte(nbsous),imint1(nbsous),ibnap(nbsous) 63 : integer inbg(nbsous),ibbas(nbsous) 64 : endsegment 65 : C 66 : segment WRK5 67 : REAL*8 W(200),A(NBN1,NBN1),B(NBN1,NBN1),XPU(3),RHS(NBN1) 68 : REAL*8 SHPXXX(6,NBNN,200),XX(3,200),VAL(NBN1),POIDS(200) 69 : INTEGER NDD 70 : endsegment 71 : C 72 : SEGMENT WRK4 73 : REAL*8 BPSS(3,3),XEL(3,NBN1),XU(3),SHP(6,NBN1),S(6) 74 : REAL*8 TXR(3,3,NBN1),XE(3,NBG) 75 : ENDSEGMENT 76 : SEGMENT VECT 77 : REAL*8 VEC1(IDIM) 78 : REAL*8 VEC2(IDIM) 79 : REAL*8 VECN(IDIM,ibbas(isous)) 80 : ENDSEGMENT 81 : 82 : C= Logique a TRUE si mode de calcul est axisymetrique 1D/2D ou Fourier 83 : LOGICAL LOGAXI 84 : C= Quelques constantes (2.Pi et 4.Pi) 85 : PARAMETER (X2Pi =6.283185307179586476925286766559D0) 86 : PARAMETER (X4Pi=12.566370614359172953850573533118D0) 87 : 88 : C on remet isup a 1 pour le moment 89 : ISUP=1 90 : MMODEL=IPMODE 91 : SEGACT,MMODEL 92 : NSOUS=KMODEL(/1) 93 : C 94 : ndec=0 94 : call lirent(ndec,0,iret) 95 : if(ndec.eq.0) ndec = 3 96 : * 97 : NBNOE = XCOOR(/1)/(IDIM+1) 98 : 99 : if(ipche.NE.0) then 100 : MCHELM=IPCHE 101 : segact MCHELM 102 : NCHAM=ICHAML(/1) 103 : endif 104 : * 105 : * Création du maillage provisoire support pour interpolation 106 : * 107 : NBSOUS = NSOUS 108 : NBREF = 0 109 : NBELEM = 0 110 : NBNN = 0 111 : C ipt2 contiendra les super couche ou points de gauss 112 : SEGINI IPT2 113 : SEGINI SICOQ 114 : c listmots des phases 115 : ilphmo = -1 116 : jgn = 8 117 : jgm = nsous 118 : segini mlmots 119 : ilphmo = mlmots 120 : jgm = 1 121 : * 122 : * Boucle sur l'ensemble des sous zones du modeles sur lequel on va 123 : * projetter pour initialisations 124 : nponou=0 125 : npcoq =0 126 : DO 30 ISOUS = 1, NSOUS 127 : * 128 : IMODEL = KMODEL(ISOUS) 129 : SEGACT,IMODEL 130 : MELEME = IMAMOD 131 : SEGACT MELEME 132 : MELE = NEFMOD 133 : NBELE1= NUM(/2) 134 : NBN1 = NUM(/1) 135 : ibnap(isous)=1 136 : C 137 : C COQUE INTEGREE OU PAS ? 138 : C 139 : IF(INFMOD(/1).NE.0)THEN 140 : NPINT=INFMOD(1) 141 : if( npint.ne.0)IBNAP(isous) =NPINT 142 : ELSE 143 : NPINT=0 144 : C ibnap(isous)=1 145 : ENDIF 146 : 147 : if( infmod(/1).lt.2+isup) then 147 : call elquoi(mele,npint,isup,ipinf,imodel) 148 : if (ierr.ne.0) then 149 : segdes imodel,mmodel 150 : segsup mchelm 151 : iret=0 152 : return 153 : endif 154 : info=ipinf 155 : mfr =infell(13) 156 : iele =infell(14) 157 : minte=infell(11) 158 : minte1=infell(12) 159 : nbg=infell(6) 160 : segsup info 161 : else 162 : minte=infmod(isup+2) 163 : minte1=infMOD(12) 164 : nbg=infele(6) 165 : mfr =infele(13) 166 : iele =infele(14) 167 : endif 168 : imfr(isous)=mfr 169 : iiele(isous)=iele 170 : iminte(isous)=minte 171 : imint1(isous)=minte1 172 : inbg(isous)=nbg 173 : * write(6,*) 'mfr iele mele minte minte1 nbg',mfr,iele, 174 : * & mele,minte,minte1,nbg 175 : NBSOUS = 0 176 : NBREF = 0 177 : NBELEM = NBELE1 178 : 179 : IF(mfr.eq.3.or.mfr.eq.5.or.mfr.eq.9) THEN 180 : C------------------coques 181 : IF(IPCHE.eq.0) THEN 182 : C pour le coques il faut les caracteristiques 182 : call erreur(404) 183 : segdes meleme ,imodel,mmodel 184 : * segsup info 185 : return 186 : ENDIF 187 : C 188 : icocoq(isous)=1 189 : nbnn = ndec*ndec 190 : if(iele.lt.4) nbnn=ndec 191 : ibbas(isous)= nbnn 192 : if(idim.eq.3.or.(iele.lt.4.and.idim.eq.2)) then 193 : ibnap(isous)=ndec 194 : nbnn= nbnn*ndec 195 : inbg(isous)=nbnn 196 : endif 197 : 198 : ELSE 199 : C--------------- massifs 200 : 201 : IF (mfr.EQ.1) THEN 202 : IF (IDIM.EQ.3) THEN 203 : NBNN=ndec*ndec*ndec 204 : ELSE IF (IDIM.EQ.2) THEN 205 : NBNN=ndec*ndec 206 : C* ELSE IF (IDIM.EQ.1) THEN 207 : ELSE 208 : NBNN=ndec 209 : ENDIF 210 : ELSE 211 : NBNN=ndec*ndec 212 : ENDIF 213 : 214 : C petite rectif utile pour les elements speciaux 215 : if(iele.eq.32) then 216 : C les polygones 217 : nbnn =nbnn*nbn1 218 : elseif(iele.eq.25.or.iele.eq.26) then 219 : C les pyramides 220 : nbnn = nbnn*4 221 : endif 222 : C 223 : inbg(isous)= NBNN 224 : ENDIF 225 : * 226 : SEGINI IPT1 227 : IPT1.ITYPEL = 28 228 : ipt2.lisous(isous) = ipt1 229 : segdes ipt1 230 : nponou=nponou+nbnn*nbelem 231 : 232 : if (isous.eq.1) then 233 : mots(1) = conmod(17:24) 234 : else 235 : do ipl = 1,jgm 236 : if (mots(ipl).eq.conmod(17:24)) goto 27 237 : enddo 238 : jgm = jgm + 1 239 : mots(jgm) = conmod(17:24) 240 : 27 continue 241 : endif 242 : 243 : segdes meleme,imodel 244 : * segsup info 245 : 30 CONTINUE 246 : segadj mlmots 247 : C----------------------------------------------------------- 248 : C 249 : C 250 : C on va fabriquer des points support provisoires a l interieur 251 : C de chaque element des maillages du modele recepteur 252 : C 253 : C 254 : C---------------------------------------------------------- 255 : C write(6,*) ' nombre de points crees ' ,nponou 256 : C (DIP) ajustement de MCOORD pour les noeuds provisoires 257 : C il sera remis a sa taille initiale a la fin du programme 258 : NBPTS = NBNOE + nponou 259 : SEGADJ MCOORD 260 : *----------------------------------------------------------- 261 : * Boucle sur l'ensemble des sous zones du modeles 262 : *------------------------------------------------------------ 263 : inupo =NBNOE 264 : DO 100 ISOUS = 1,NSOUS 265 : IMODEL = KMODEL(ISOUS) 266 : segact imodel 267 : MELEME = IMAMOD 268 : segact meleme 269 : nbele1= num(/2) 270 : nbn1 = num(/1) 271 : ipt1=ipt2.lisous(isous) 272 : segact ipt1*mod 273 : NNB = ibbas(isous) 274 : NBG = inbg(isous) 275 : SEGINI WRK4 276 : MELVA1 = 0 277 : MELVA2 = 0 278 : IF(icocoq(isous).ne.0) THEN 279 : C -------------------------------la sz est une coque ------ 280 : 281 : NUCHA = 0 282 : DO 15, NUCH = 1, NCHAM 283 : * 284 : IF ( CONCHE(NUCH).EQ.CONMOD.AND. 285 : & IMACHE(NUCH).EQ.IMAMOD) NUCHA = NUCH 286 : * 287 : 15 CONTINUE 288 : * 289 : IF (NUCHA.NE.0) THEN 290 : MCHAML=ICHAML(NUCHA) 291 : SEGACT,MCHAML 292 : * 293 : MELVA1 = 0 294 : MELVA2 = 0 295 : NCOMP = IELVAL(/1) 296 : DO 20, I = 1, NCOMP 297 : IF (NOMCHE(I)(1:4).EQ.'EPAI') THEN 298 : MELVA1 = IELVAL(I) 299 : SEGACT, MELVA1 300 : ELSEIF (NOMCHE(I)(1:4).EQ.'EXCE') THEN 301 : MELVA2 = IELVAL(I) 302 : SEGACT, MELVA2 303 : ENDIF 304 : 20 CONTINUE 305 : ENDIF 306 : ENDIF 307 : C 308 : MELE = NEFMOD 309 : c coque integree ou pas ? 310 : if(infmod(/1).ne.0)then 311 : npint=infmod(1) 312 : else 313 : npint=0 314 : endif 315 : C 316 : nbnap=ibnap(isous) 317 : mfr =imfr(isous) 318 : iele =iiele(isous) 319 : minte=iminte(isous) 320 : minte1=imint1(isous) 321 : nbg = inbg(isous) 322 : C 323 : segact minte 324 : if(mfr.eq.5) segact minte1 325 : 326 : nnb = ibbas(isous) 327 : C 328 : ipo=0 329 : *------------------------------------------------------------------- 330 : * Pavage de l element de reference 331 : *------------------------------------------------------------------- 332 : NBNN = inbg(isous) 333 : segini WRK5 334 : iwkr= wrk5 335 : iele=iiele(isous) 335 : CALL PAVPOI(WRK5,iele,ndec,nbn1) 336 : segini vect 337 : *------------------------------------------------------------------- 338 : * Boucle sur les elements de la sous zone du modele 339 : *------------------------------------------------------------------- 340 : 341 : DO IEL=1,NBELE1 341 : CALL DOXE(XCOOR,IDIM,NBN1,NUM,IEL,XEL) 342 : 343 : C si c est une coque pouvant etre excentree 344 : C recherche des normales aux noeuds 345 : if(icocoq(isous).ne.0) then 346 : C 347 : IF(IELE.EQ.4) THEN 348 : C --------coq3 dkt dst ---------------- 349 : do 98 i=1,idim 350 : VEC1(i)=xel(i,2)-xel(i,1) 351 : 98 VEC2(i)=xel(i,3)-xel(i,1) 352 : vecn(1,1) = vec1(2)*vec2(3)-vec1(3)*vec2(2) 353 : vecn(2,1) = vec1(3)*vec2(1)-vec1(1)*vec2(3) 354 : vecn(3,1) = vec1(1)*vec2(2)-vec1(2)*vec2(1) 355 : vnor =sqrt(vecn(1,1)*vecn(1,1)+vecn(2,1)*vecn(2,1) 356 : & +vecn(3,1)*vecn(3,1)) 357 : do 91 i=1,idim 358 : 91 vecn(i,1)=vecn(i,1)/vnor 359 : do 99 i=1,idim 360 : do 99 j=2,ibbas(isous) 361 : 99 vecn(i,j)=vecn(i,1) 362 : ELSEIF (IELE.EQ.8) THEN 363 : C --------------- coq4 ---------------- 363 : call cq4loc(xel,xe,BPSS,irrt,1) 364 : C ici c est une estimation du plan moyen 365 : * write(6,fmt='(3E12.5)') ((bpss(i,j),i=1,3),j=1,3) 366 : do i=1,3 367 : do ip=1,ibbas(isous) 368 : vecn(i,ip)=bpss(3,i) 369 : enddo 370 : enddo 371 : ELSEIF(IELE.eq.2) THEN 372 : C ---------------- coq2 -------------- 373 : vnor=0.D0 374 : do 92 i=1,idim 375 : VEC1(i)=xel(i,2)-xel(i,1) 376 : 92 vnor=vnor + vec1(i)*vec1(i) 377 : vnor = sqrt(vnor) 378 : vecn(1,1)=-vec1(2)/vnor 379 : vecn(2,1)=vec1(1)/vnor 380 : vecn(1,2)=vecn(1,1) 381 : vecn(2,2)=vecn(2,1) 382 : ELSEIF (iele.eq.5.or.iele.eq.10) THEN 383 : C preparation pour coques integrees 384 : SURFX=XZERO 385 : SURFY=XZERO 386 : SURFZ=XZERO 387 : ENDIF 388 : endif 389 : C ------------------- boucle sur les points d integration 390 : inuo=inupo 391 : C write(6,*) ' ndd et inupo ' ,ndd,inupo 392 : do igau=1,ndd 393 : C 394 : do j=1,idim 395 : xpu(j)=xx(j,igau) 396 : enddo 397 : C write(6,*)(xpu(J),j=1,idim) 398 : C 399 : if(mele.GE.111.and.mele.le.122) then 400 : C les plolygones 400 : call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret) 401 : else 401 : CALL SHAPE(xpu(1),xpu(2),xpu(3),iele,shp,iret) 402 : endif 403 : 404 : C PROJECTION DANS L ELEMENT REEL 405 : DO 501 IK = 1,IDIM 406 : XU(IK)= XZERO 407 : DO 501 JK = 1,nbn1 408 : XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK)) 409 : 501 CONTINUE 410 : C 411 : if ( mfr.eq.5) then 412 : C coques integrees ds l epaisseur 413 : DO I=1,6 414 : S(I)=XZERO 415 : CONTINUE 416 : ENDDO 417 : DO INOE=1,NBN1 418 : S(1)=S(1)+SHP(2,INOE)*XEL(2,INOE) 419 : S(2)=S(2)+SHP(3,INOE)*XEL(3,INOE) 420 : S(3)=S(3)+SHP(3,INOE)*XEL(2,INOE) 421 : S(4)=S(4)+SHP(2,INOE)*XEL(3,INOE) 422 : S(5)=S(5)+SHP(3,INOE)*XEL(1,INOE) 423 : S(6)=S(6)+SHP(2,INOE)*XEL(1,INOE) 424 : ENDDO 425 : SURFX=S(1)*S(2)-S(3)*S(4) 426 : SURFY=S(4)*S(5)-S(2)*S(6) 427 : SURFZ=S(6)*S(3)-S(5)*S(1) 428 : SURF=SQRT(SURFX**2+SURFY**2+SURFZ**2) 429 : C write(6,*) ' normales au pt igau du plan de base ' 430 : vecn(1,igau)=surfx/surf 431 : vecn(2,igau)=surfy/surf 432 : vecn(3,igau)=surfz/surf 433 : endif 434 : C 435 : IF(ICOCOQ(ISOUS).EQ.0) THEN 436 : C elements standards 437 : inupo=inupo+1 438 : ipt1.num(igau,iel)=inupo 439 : do i=1,idim 440 : xcoor((inupo-1)*(idim+1)+i)=XU(i) 441 : enddo 442 : C write(6,3335) (xu(i),i=1,idim) 443 : C3335 format(3e12.5) 444 : ELSE 445 : C traitement special des coques 446 : if(melva1.ne.0) then 447 : ibmn =MIN(iel,melva1.velche(/2)) 448 : igmn =MIN(igau,melva1.velche(/1)) 449 : epai = melva1.velche(igmn,ibmn) 450 : else 451 : epai =xzero 452 : endif 453 : C 454 : if(melva2.ne.0) then 455 : ibmn =MIN(iel,melva2.velche(/2)) 456 : igmn =MIN(igau,melva2.velche(/1)) 457 : exce = melva2.velche(igmn,ibmn) 458 : else 459 : exce =xzero 460 : endif 461 : 462 : dp = epai/ibnap(isous) 463 : dnor = exce -(epai + dp)/2.D0 464 : C 465 : do inap=1,ibnap(isous) 466 : dnor=dnor+dp 467 : ipos = (inap-1)*ibbas(isous)+igau 468 : ipt1.num(ipos,iel)=inuo+ipos 469 : do i=1,idim 470 : xcoor((inuo+ipos-1)*(idim+1)+i)= XU(i)+dnor*vecn(i,igau) 471 : enddo 472 : enddo 473 : inupo=inupo+ibnap(isous) 474 : ENDIF 475 : 476 : C ---------------fin de la boucle sur les points de d integration 477 : enddo 478 : 479 : C fin de la boucle sur les elements 480 : enddo 481 : segsup wrk5,vect 482 : segdes minte 483 : if(mfr.eq.5) segdes minte1 484 : if(icocoq(isous).ne.0) SEGDES,MCHAML 485 : if (melva1.ne.0) segdes melva1 486 : if (melva2.ne.0) segdes melva2 487 : segdes meleme,imodel 488 : segdes ipt1 489 : SEGSUP WRK4 490 : * segsup info 491 : * 492 : 100 CONTINUE 493 : C fin de creation des noeuds supports provisoires 494 : 2003 format(3I4,2X,4e12.5,I4) 495 : C--------------//////////////////////////////////------------------ 496 : C maintenant in faut un meleme poi1 de tous les points sur lesquels 497 : C on doit interpoler une valeur 498 : segini icpr1 499 : ipt5=ipt2 500 : nbelem=0 501 : do io=1,max(1,ipt2.lisous(/1)) 502 : if (ipt2.lisous(/1).ne.0) ipt5=ipt2.lisous(io) 503 : segact ipt5 504 : do ip=1,ipt5.num(/1) 505 : do il=1,ipt5.num(/2) 506 : ipt=ipt5.num(ip,il) 507 : if (icpr1(ipt).eq.0) then 508 : nbelem=nbelem+1 509 : icpr1(ipt)=nbelem 510 : endif 511 : enddo 512 : enddo 513 : ** segdes ipt5 514 : enddo 515 : ** segdes ipt2 516 : nbnn=1 517 : nbsous=0 518 : nbref=0 519 : ** write (6,*) ' nbelem ',nbelem 520 : segini ipt4 521 : ipgeom=ipt4 522 : nbelem=0 523 : do i=1,icpr1(/1) 524 : if (icpr1(i).ne.0) then 525 : nbelem=nbelem+1 526 : ** write (6,*) ' i nbelem icpr ',i,nbelem,icpr1(i) 527 : ipt4.num(1,nbelem)=i 528 : endif 529 : enddo 530 : segsup icpr1 531 : C---------------------------------------------------------------- 532 : C---------------------------------------------------------------- 533 : *-------------on est pret a faire l interpolation sur ipgeom ---- 534 : C 535 : *write(6,*) ' maillage apres changer poi1 ' ,ipgeom 535 : C call ecmail(ipgeom) 536 : isort=0 536 : CALL PRO2(IPGEOM,IPCHT,isort,IPOUT,ilphmo) 537 : if (ierr.ne.0) return 538 : C---------------------------------------------------------------- 539 : C---------------------------------------------------------------- 540 : C 541 : C Definition de "constantes" selon le mode de calcul 542 : C Modes de calcul 1D/2D axisymetrique 543 : IF (IFOMOD.EQ.0.OR.IFOMOD.EQ.4) THEN 544 : LOGAXI=.TRUE. 545 : CteAxi=X2Pi 546 : C Mode de calcul 2D Fourier 547 : ELSE IF (IFOMOD.EQ.1) THEN 548 : LOGAXI=.TRUE. 549 : IF (NIFOUR.EQ.0) THEN 550 : CteAxi=X2Pi 551 : ELSE 552 : CteAxi=XPi 553 : ENDIF 554 : ELSE 555 : LOGAXI=.FALSE. 556 : ENDIF 557 : C---------------------------------------------------------------- 558 : C --- il faut maintenant reconstituer les MCHAML 559 : C --- a partir du chpo construit sur ipgeom ---------------- 560 : mlchpo = ipout 561 : segact mlchpo 562 : mchel5= ipcht 563 : segact mchel5 564 : N1=NSOUS 565 : L1=12 566 : N3=6 567 : SEGINI MCHEL1 568 : * MCHEL1.TITCHE='SCALAIRE' 569 : MCHEL1.TITCHE = mchel5.titche 570 : MCHEL1.IFOCHE=IFOUR 571 : segdes mchel5 572 : *////////////////////////////// 573 : 574 : * boucle phases modele 575 : DO 7000 IPHAS = 1,JGM 576 : * write(6,*) 'iphas ',iphas 577 : segini icpr,snomip 578 : * mdata=isort 579 : * segact mdata 580 : mchpoi=ichpoi(iphas) 581 : segact mchpoi 582 : nsschp=ipchp(/1) 583 : segini ipcb 584 : 585 : C ip=0 586 : do i=1,nsschp 587 : * mchpoi=ipca(i) 587 : C call ecchpo(mchpoi) 588 : * segact mchpoi 589 : msoupo=ipchp(i) 590 : ipcb(i)=msoupo 591 : segact msoupo 592 : mpoval=ipoval 593 : segact mpoval 594 : nc = nocomp(/2) 595 : C 596 : C tableau des general des composantes 597 : do ic = 1,nc 598 : nomm= nocomp(ic) 599 : if(nomip(/2).eq.0) then 600 : nomip(**)=nomm 601 : else 602 : do k=1,nomip(/2) 603 : if(nomm.eq.nomip(k)) goto 4322 604 : enddo 605 : nomip(**)=nomm 606 : 4322 continue 607 : endif 608 : enddo 609 : * write(6,*) (nomip(l),l=1,nomip(/2)) 610 : * reperage de la position ses points dans le chpo 611 : ipt5=igeoc 612 : segact ipt5 613 : * ip=0 614 : do j=1,ipt5.num(/2) 615 : * ip=ip+1 616 : icpr(1,ipt5.num(1,j))=j 617 : icpr(2,ipt5.num(1,j))= i 618 : enddo 619 : enddo 620 : C i=icpr(2,k) point k venant du msoupo du msoupo i dans ipcb 621 : 622 : inupo=0 623 : segact ipt2 624 : DO 2100 isous = 1,nsous 625 : imodel= KMODEL(isous) 626 : segact imodel 627 : if (conmod(17:24).ne.mots(iphas)) goto 2100 628 : meleme = imamod 629 : segact meleme 630 : idi2=idim 631 : ipt5=ipt2.lisous(isous) 632 : segact ipt5 633 : if(icocoq(isous).ne.0) idi2 = idim-1 634 : mfr=imfr(isous) 635 : * création du nouveau segment MCHAML 636 : N2=nomip(/2) 637 : SEGINI MCHAML 638 : MCHEL1.IMACHE(isous)=MELEME 639 : MCHEL1.ICHAML(isous)=MCHAML 640 : MCHEL1.CONCHE(isous)=CONMOD 641 : MCHEL1.INFCHE(isous,1)=0 642 : MCHEL1.INFCHE(isous,2)=0 643 : MCHEL1.INFCHE(isous,3)=NIFOUR 644 : MCHEL1.INFCHE(isous,4)=0 645 : MCHEL1.INFCHE(isous,5)=0 646 : MCHEL1.INFCHE(isous,6)=1 647 : * 648 : N1EL=NUM(/2) 649 : N2PTEL=0 650 : N2EL=0 651 : C-------------------------------------------------------- 652 : C boucle sur les composantes 653 : C-------------------------------------------------------- 654 : NBNN=num(/1) 655 : N1PTEL= NBNN 656 : NBN1=NBNN 657 : 658 : C 659 : SEGINI WRK5,wrk4 660 : iele = iiele(isous) 661 : iwrk= wrk5 661 : CALL PAVPOI(WRK5,iele,ndec,nbnn) 662 : C ----- remplissage de shpxxx ( shptot) 663 : DO IGAU=1,NDD 664 : do j=1,idim 665 : xpu(j)=xx(j,igau) 666 : enddo 667 : if(mele.GE.111.and.mele.le.122) then 667 : call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret) 668 : else 668 : call shape(xpu(1),xpu(2),xpu(3),iele,shp,iret) 669 : endif 670 : do id=1,6 671 : do no=1,nbnn 672 : shpxxx(id,no,igau)=shp(id,no) 673 : enddo 674 : enddo 675 : ENDDO 676 : N1PTEL=NBNN 677 : C ----------- creation des melvals recepteurs ----- 678 : DO icomp=1,NOMIP(/2) 679 : SEGINI MELVAL 680 : NOMCHE(icomp)=nomip(icomp) 681 : TYPCHE(icomp)='REAL*8' 682 : IELVAL(icomp)=MELVAL 683 : ENDDO 684 : C-------------------------------------------------- 685 : C-------- boucle sur les elements ---------- 686 : AXIS=O1 687 : DO 2210 iel=1,num(/2) 687 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IEL,XEL) 687 : CALL ZERO(A,NBNN,NBNN) 688 : C SSOM1=XZERO 689 : C 690 : C 691 : DO inap=1,ibnap(isous) 692 : DO 40 IGAU=1,ndd 693 : do i=1,idim 694 : xpu(i)=xx(i,igau) 695 : enddo 696 : if(mele.GE.111.and.mele.le.122) then 697 : C les plolygones 697 : call shpoly(xpu(1),xpu(2),xpu(3),nbn1,mele,shp,iret) 698 : else 698 : call shape(xpu(1),xpu(2),xpu(3),iele,shp,iret) 699 : endif 699 : CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC) 700 : C Modes de calcul 1D/2D axisymetrique et 2D Fourier 701 : IF (LOGAXI) THEN 701 : CALL DISTRR(XEL,SHP,NBNN,RR) 702 : AXIS=CteAxi*RR 703 : C Mode de calcul 1D spherique 704 : ELSE IF (IFOMOD.EQ.5) THEN 704 : CALL DISTRR(XEL,SHP,NBNN,RR) 705 : AXIS=X4Pi*RR*RR 706 : ENDIF 707 : C 708 : C on veut calculer la fonction au point d integration projete 709 : C PROJECTION 710 : DO IK = 1,IDIM 711 : XU(IK)= XZERO 712 : DO JK = 1,NBNN 713 : XU(IK) = XU(IK) + SHP(1,JK)*(XEL(IK,JK)) 714 : ENDDO 715 : ENDDO 716 : C 717 : CPZA=poids(igau)*DJAC*AXIS 718 : DO 50 INO1=1,NBNN 719 : DO 60 INO2=1,NBNN 720 : Q = SHPXXX(1,INO1,IGAU)*SHPXXX(1,INO2,IGAU)*CPZA 721 : CCC Q = SHPXXX(1,INO1,IGAU)*SHPXXX(1,INO2,IGAU) 722 : CCC 1 *poids(igau)*DJAC*AXIS 723 : A(INO1,INO2)=A(INO1,INO2)+Q 724 : 60 CONTINUE 725 : 50 CONTINUE 726 : 40 CONTINUE 727 : enddo 728 : C 729 : C ---- fin preparation matrice 'masse ' sur element courant 730 : C do i=1,NBNN 731 : C do j=1,NBNN garder eventuellement pour calculer 732 : C B(i,j)=A(i,j) la difference sur les integrales 733 : C enddo 734 : C enddo 735 : C 735 : call invalm(A,NBNN,NBNN,KERRE,1.D-10) 736 : C ------------------ boucle sur les composantes 737 : DO 2220 icomp=1,nomip(/2) 738 : idebco = inupo 739 : do i=1,nbnn 740 : RHS(i) =XZERO 741 : enddo 742 : C --------- boucle sur les points d integration 743 : C 744 : AXIS = O1 745 : DO 42 inap=1,ibnap(isous) 746 : DO 41 IGAU=1,ndd 747 : do i=1,idim 748 : xpu(i)=xx(i,igau) 749 : enddo 749 : CALL SHAPE(xpu(1),xpu(2),xpu(3),iele,shp,iret) 749 : CALL JACOBI(XEL,SHP,IDI2,NBNN,DJAC) 750 : C Modes de calcul 1D/2D axisymetrique et 2D Fourier 751 : IF (LOGAXI) THEN 751 : CALL DISTRR(XEL,SHP,NBNN,RR) 752 : AXIS=CteAxi*RR 753 : C Mode de calcul 1D spherique 754 : ELSE IF (IFOMOD.EQ.5) THEN 754 : CALL DISTRR(XEL,SHP,NBNN,RR) 755 : AXIS=X4Pi*RR*RR 756 : ENDIF 757 : 758 : ipop=(inap-1)*ndd+igau 759 : inupo=ipt5.num(ipop,iel) 760 : 761 : jh=icpr (1,inupo) 762 : ipa=icpr(2,inupo) 763 : C il faut verifier si vpocha existe pour ce point 764 : msoupo = ipcb(ipa) 765 : mpoval=ipoval 766 : do l=1,nocomp(/2) 767 : if(nocomp(l).eq.nomip(icomp)) then 768 : vvv=vpocha(jh,icomp) 769 : goto 4555 770 : endif 771 : enddo 772 : vvv=XZERO 773 : 4555 continue 774 : C write(6,7777) inupo,jh,poids(igau),vvv 775 : C 776 : melval=ielval(icomp) 777 : CPZA=vvv*poids(igau)*DJAC*AXIS 778 : DO INO1 = 1,NBNN 779 : RHS(ino1)=RHS(ino1)+CPZA*SHPXXX(1,ino1,igau) 780 : CCC RHS(ino1)=RHS(ino1)+vvv*poids(igau)* 781 : CCC 1 DJAC*SHPXXX(1,INO1,IGAU)*AXIS 782 : ENDDO 783 : 41 continue 784 : 42 continue 785 : 786 : 7777 format(2I4,2X,2E12.5) 787 : C SSOM1 =SSOM1+SURF*poids(igau)*vvv 788 : C 789 : do i=1,NBNN 790 : val(i)=XZERO 791 : do j=1,NBNN 792 : val(i)=val(i)+A(i,j)*RHS(J) 793 : enddo 794 : enddo 795 : C 796 : do I=1,NBNN 797 : velche(i,iel)=val(i) 798 : enddo 799 : C 800 : * if(icomp.lt.nomip(/2)) inupo =idebco 801 : C ---- fin de la boucle sur les composantes 802 : 2220 continue 803 : C fin de la boucle sur les elements 804 : 2210 continue 805 : C 806 : do i=1,ielval(/1) 807 : ipc1 = ielval(i) 807 : call comred(ipc1) 808 : melval = ipc1 809 : segdes melval 810 : enddo 811 : segdes meleme,mchaml 812 : 813 : C---------------- fin du traitement des massifs 814 : segsup wrk5,wrk4,ipt5 815 : 816 : C segsup ipt3 817 : segdes imodel,meleme 818 : C fin de la boucle sur les sous zones du modele 819 : 2100 continue 820 : C destruction des chpo intermediaires 821 : do i=1,ipcb(/1) 822 : msoupo=ipcb(i) 823 : cgf ipt5= igeoc (correction 7284) 824 : mpoval= ipoval 825 : * mchpoi=ipca(i) 826 : segsup mpoval,msoupo 827 : cgf segsup ipt5 (correction 7284) 828 : enddo 829 : segsup mchpoi 830 : segsup ipcb,snomip 831 : segsup icpr 832 : 7000 CONTINUE 833 : C (fdp) re-ajustement de MCOORD a sa taille initiale 834 : NBPTS = NBNOE 835 : SEGADJ MCOORD 836 : C retrait des maillages temporaires du pre-conditionnement 837 : c (leurs numero de noeuds depasse la taille de MCOOR) 837 : call crech1b 838 : segdes mchel1,mmodel 839 : C (fdp) suppression du maillage temporaire IPT1 840 : ipt1=ipgeom 841 : segsup ipt1 842 : segact ipt2 843 : do io=1,ipt2.lisous(/1) 844 : ipt5=ipt2.lisous(io) 845 : * segsup ipt5 846 : enddo 847 : segsup ipt2 848 : segsup ipt4 849 : if(ipche.ne.0) then 850 : mchelm=ipche 851 : segdes mchelm 852 : endif 853 : IPOUT=MCHEL1 854 : segsup sicoq 855 : 856 : return 857 : END 858 : 859 : 860 : 861 : 862 :
© Cast3M 2003 - All rights reserved.
Disclaimer