Télécharger zoneg2.eso

Retour à la liste

Numérotation des lignes :

zoneg2
  1. C ZONEG2 SOURCE OF166741 25/07/28 21:15:16 12336
  2. SUBROUTINE ZONEG2(IP1,IP2,ICOUNT,IFL)
  3. C====================================================================
  4. C localisation de points dans des elements d une zone elementaire
  5. C zonage et construction d un tableau de correspondance
  6. C
  7. C Entrees :
  8. C ip1 pointeur sur le maillage massif
  9. C ip2 pointeur sur le maillage poi1
  10. C iexx pointeur sur un segment contenant les points deja vus
  11. C ( eviter comptage des point aux frontieres de 2 sous zones)
  12. C IFL = 0 noeuds communs a IP1 et IP2 non repertories( accro)
  13. C IFL = 1 "" comptabilisés ( proi)
  14. C Sorties :
  15. C iccoun pointeur sur le segment ICOUNT contenant le tableau de
  16. C correspondance voir commentaire pro2
  17. C====================================================================
  18. IMPLICIT INTEGER(I-N)
  19. IMPLICIT REAL*8 (A-H,O-Z)
  20.  
  21. -INC PPARAM
  22. -INC CCOPTIO
  23. -INC CCREEL
  24. -INC CCGEOME
  25.  
  26. -INC SMCOORD
  27. -INC SMELEME
  28. SEGMENT /STRAV/(NP1(ITE),NP2(ITE),NP3(ITE),NPI(ITE))
  29. dimension xpu(3)
  30.  
  31. PARAMETER (us3=1.D0/3.D0)
  32. PARAMETER (Un=1.D0)
  33.  
  34. segment inomil(0)
  35. segment wrk4
  36. real*8 xel(3,nbnn),shpp(6,nbnn),qsi(3)
  37. endsegment
  38. segment wrkl
  39. real*8 xell(3,nbnnl),shpl(6,nbnnl)
  40. endsegment
  41.  
  42. SEGMENT ICOUNT
  43. INTEGER IEINT(2,NODES),IDEJVU(NODES),IACCRO
  44. REAL*8 QQQ(3,NODES),CKRIT(NODES)
  45. ENDSEGMENT
  46. EXTERNAL SHAPE
  47. REAL*8 ZDIST
  48.  
  49. ZDIST=0.D0
  50. ip=1
  51. ipo=1
  52. * WRITE(IOIMP,*) ' entree dans zoneg2'
  53. itr=0
  54. irate=0
  55. IDIM1 = IDIM + 1
  56. IPT2 = IP2
  57. ipt1=ip1
  58.  
  59. * on fait le zonage sur les points
  60. SEGACT IPT2
  61. Nodes = IPT2.NUM(/2)
  62. * WRITE(IOIMP,*) ' nodes ' , nodes
  63.  
  64. numnp=nodes
  65. ITE=numnp+3+1
  66. numnp3=ite-1
  67. * write(6,*) ' numnp ite numnp3 ' , numnp , ite , numnp3
  68. NF=100
  69. NMULT=100
  70. segini strav
  71. C
  72. C Cas 1D : idem qu'en 2D ou 3D mais pour eviter multiplications
  73. C -------- des tests (IF) qui alourdissent, on duplique ZONEG2
  74. C pour le 1D en 600
  75. IF (IDIM.EQ.1) GOTO 600
  76.  
  77. C recherche des min et max
  78. C
  79. iref = ipt2.NUM(1,1)*idim1-idim
  80. XMIN=xcoor(iref)
  81. YMIN=xcoor(iref+1)
  82. ZMIN=xcoor(iref+2)
  83. YMAx=ymin
  84. XMAx=xmin
  85. ZMAx=zmin
  86. DO 201 Iyr = 2,nodes
  87. IB=ipt2.NUM(1,Iyr)
  88. IA=IB*IDIM1-idim
  89. IF(XCOOR(IA).LT.XMIN) XMIn=XCOOR(IA)
  90. IF(XCOOR(IA).GT.XMAx) XMAx=XCOOR(IA)
  91. IF(XCOOR(IA+1).LT.YMIn) YMIn=XCOOR(IA+1)
  92. IF(XCOOR(IA+1).GT.YMAx) YMAx=XCOOR(IA+1)
  93. IF( IDIM.EQ.3 ) THEN
  94. IF(XCOOR(IA+2).LT.ZMIn) ZMIn=XCOOR(IA+2)
  95. IF(XCOOR(IA+2).GT.ZMAx) ZMAx=XCOOR(IA+2)
  96. ENDIF
  97. 201 CONTINUE
  98. if(idim.eq.2) then
  99. ZMIN=0.D0
  100. ZMAX=0.D0
  101. ZDIST=1.D0
  102. ZDISTV=1.D0
  103. endif
  104. XDIST=XMAx-XMIn
  105. YDIST=YMAx-YMIn
  106. xdistv=xdist
  107. ydistv=ydist
  108. IF(XDIST.EQ.0.D0) THEN
  109. NDEX=1
  110. XDIST=1.D0
  111. ELSE
  112. NDEX=0
  113. ENDIF
  114. IF(YDIST.EQ.0.D0) THEN
  115. NDEY=1
  116. YDIST=1.D0
  117. ELSE
  118. NDEY=0
  119. ENDIF
  120. * SMO=1E-5*XDIST*YDIST/NUMNP
  121. IF( IDIM.EQ.3 ) then
  122. zdist = zmax-zmin
  123. zdistv=zdist
  124. IF(ZDIST.EQ.0.D0) THEN
  125. NDEZ=1
  126. ZDIST=1.D0
  127. * SMO=1E-5*XDIST*YDIST*ZDIST/NUMNP
  128. ELSE
  129. NDEZ=0
  130. ENDIF
  131. endif
  132. * write(ioimp,*) 'Maillage POI1'
  133. * write(ioimp,*) 'xmin,xmax,ymin,ymax,zmin,zmax=',xmin,xmax,ymin
  134. * $ ,ymax,zmin,zmax
  135. * write(ioimp,*) 'xdistv,ydistv,zdistv=',xdistv,ydistv,zdistv
  136. * write(ioimp,*) 'xdist,ydist,zdist=',xdist,ydist,zdist
  137.  
  138. C
  139. C*** recherche sur les elements de la plus petite taille en x y et z
  140. C
  141. ipt1=ip1
  142. segact ipt1
  143. dminx=1.D10
  144. dminy=1.D10
  145. dminz=1.D10
  146. do 453 isous=1,ipt1.lisous(/1)
  147. meleme=ipt1.lisous(isous)
  148. segact meleme
  149. do 454 iel=1,num(/2)
  150. dmix=1.D10
  151. dmax=-1.D10
  152. dmiy=1.D10
  153. dmay=-1.D10
  154. dmiz=1.D10
  155. dmaz=-1.D10
  156. do 455 ip=1,num(/1)
  157. ire=num(ip,iel)*idim1-idim
  158. if(xcoor(ire).le.dmix) dmix=xcoor(ire)
  159. if(xcoor(ire).gt.dmax) dmax=xcoor(ire)
  160. if(xcoor(ire+1).le.dmiy) dmiy=xcoor(ire+1)
  161. if(xcoor(ire+1).gt.dmay) dmay=xcoor(ire+1)
  162. if(xcoor(ire+2).le.dmiz) dmiz=xcoor(ire+2)
  163. if(xcoor(ire+2).gt.dmaz) dmaz=xcoor(ire+2)
  164. 455 continue
  165. if(dminx.gt.(dmax-dmix) ) dminx= dmax-dmix
  166. if(dminy.gt.(dmay-dmiy) ) dminy= dmay-dmiy
  167. if(dminz.gt.(dmaz-dmiz) ) dminz= dmaz-dmiz
  168. 454 continue
  169. 453 continue
  170. * write(ioimp,*) 'Maillage Massif'
  171. * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz
  172. * WRITE(IOIMP,*) ' xdist,ydist,zdist',xdist,ydist,zdist
  173.  
  174. if( dminx.eq.0.D0) then
  175. if(xdistv.eq.0.D0) then
  176. dminx=1.D0
  177. else
  178. dminx = xdistv/1000.D0
  179. endif
  180. endif
  181. if( dminy.eq.0.D0) then
  182. if(ydistv.eq.0.D0) then
  183. dminy=1.D0
  184. else
  185. dminy = ydistv/1000.D0
  186. endif
  187. endif
  188. if (idim.eq.3) then
  189. if( dminz.eq.0.D0) then
  190. if(zdistv.eq.0.D0) then
  191. dminz=1.D0
  192. else
  193. dminz = zdistv/1000.D0
  194. endif
  195. endif
  196. endif
  197. *
  198. * write(ioimp,*) 'Apres correction 1'
  199. * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz
  200.  
  201. dminx = max( dminx , xdist/10000.D0)
  202. dminy = max(dminy,ydist/10000.D0)
  203. if( idim.eq.3) dminz = max(dminz,zdist/10000.D0)
  204.  
  205. * write(ioimp,*) 'Apres correction 2'
  206. * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz
  207.  
  208. if (ndex.eq.0) ndex=xdist/dminx
  209. ndex=max(1,ndex)
  210. if (ndey.eq.0) ndey=ydist/dminy
  211. ndey=max(1,ndey)
  212. if(idim.eq.3) then
  213. if (ndez.eq.0) ndez=zdist/dminz
  214. ndez=max(1,ndez)
  215. else
  216. ndez=1
  217. zdist=1.D0
  218. dminz=1.01D0
  219. endif
  220. * WRITE(IOIMP,*) ' ndex ndey ndez',ndex,ndey,ndez
  221.  
  222. C CALCUL DU DECOUPAGE EN ZONE
  223. NBZONE=NUMNP*NF
  224. * if( ndex*ndey*ndez.lt.nbzone) then
  225. * WRITE(IOIMP,*) ' ndex ndey ndez nbzone=',ndex,ndey,ndez,nbzone
  226. xnprop = real(ndex)*real(ndey)*real(ndez)
  227. * dimension effective : là où les découpages sont plus grands que 1
  228. idimf=0
  229. if (ndex.gt.1) idimf=idimf+1
  230. if (ndey.gt.1) idimf=idimf+1
  231. if (ndez.gt.1) idimf=idimf+1
  232. * write(ioimp,*) ' idimf=',idimf
  233. * if(idim.eq.3) then
  234. * rap = (real(nbzone)/xnprop)**0.333333333D0
  235. * else
  236. * rap= (real(nbzone)/xnprop)** 0.5D0
  237. * endif
  238. if (idimf.gt.0) then
  239. rap = (real(nbzone)/xnprop)**(1.D0/real(idimf))
  240. * write(ioimp,*) ' rap=',rap
  241. if(rap.gt.1.D0) rap=1.d0
  242. else
  243. rap=1.d0
  244. endif
  245. * write(ioimp,*) ' rap=',rap
  246. if (ndex.gt.1) then
  247. ndec=ndex*rap
  248. ndec=max(1,ndec)
  249. else
  250. ndec=1
  251. endif
  252. *
  253. if (ndey.gt.1) then
  254. mdec=ndey*rap
  255. mdec = max(1,mdec)
  256. else
  257. mdec=1
  258. endif
  259. *
  260. if(idim.eq.3) then
  261. if (ndez.gt.1) then
  262. ldec=ndez*rap
  263. ldec = max(1,ldec)
  264. else
  265. ldec=1
  266. endif
  267. else
  268. ldec=1
  269. endif
  270. *
  271. nbzone=ndec*mdec*ldec
  272. nf= nbzone / numnp3 +1
  273. nmult=nf
  274. * write(ioimp,*) 'Apres correction 3'
  275. * WRITE(IOIMP,*) ' ndec mdec ldec nbzone=',ndec,mdec,ldec,nbzone
  276. * else
  277. * RAPY=YDIST/XDIST
  278. * NDEC=MAX(INT(SQRT(REAL(NBZONE)/RAPY)),1)
  279. * MDEC=NBZONE/NDEC
  280. * LDEC=1
  281. * IF(IDIM.EQ.3) THEN
  282. * NDEC=max(int((real(nbzone)*xdist*xdist/(ydist*zdist))
  283. * $ **0.3333333D0),1)
  284. * MDEC= max(int((real(nbzone)*ydist*ydist/(xdist*zdist))
  285. * $ **0.3333333D0),1)
  286. * LDEC=NBZONE/(MDEC*NDEC)
  287. * ENDIF
  288. NBZONE=NDEC*MDEC*LDEC
  289. * endif
  290. * WRITE (IOIMP,8933) NDEC,MDEC,ldec,NBZONE,nmult
  291. * 8933 FORMAT(' DECOUPAGE EN X ',I6,' EN Y ',I6,' EN Z ',i6,
  292. * $ ' SOIT ',I8,' ZONES ')
  293. * write(6,*) 'diviseur des zones pour aller aux super zones',nmult
  294. * write(6,*) ' nombre max de super zones ', nbzone/nmult+1
  295. XDEC=XDIST/NDEC*1.0000001D0
  296. YDEC=YDIST/MDEC*1.0000001D0
  297. ZDEC=ZDIST/LDEC*1.0000001D0
  298. * WRITE (IOIMP,8935) XMIN,XMAX,YMIN,YMAX,zmin, zmax,XDEC,YDEC,zdec
  299. * 8935 FORMAT (' XMIN,XMAX,YMIN,YMAX,XDEC,YDEC ',2(6G12.5))
  300. NMDEC=NDEC*MDEC
  301. if(idim.eq.2) THEN
  302. DO 2 I=1,NUMNP
  303. IB=IPT2.num(1,i)
  304. IA=IB*IDIM1-idim
  305. IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+
  306. & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC
  307. IPOINT=IZONE/NMULT+1
  308. * WRITE(IOIMP,*) ' point zone ipoint' ,ib,izone,ipoint
  309. np1(ipoint)=np1(ipoint)+1
  310. 2 CONTINUE
  311. ELSE
  312. DO 3 I=1,NUMNP
  313. IB=IPT2.num(1,i)
  314. IA=IB*IDIM1 -idim
  315. IZONE=INT(real((XCOOR(IA)-XMIn)/XDEC))+1+
  316. & INT(real((XCOOR(IA+1)-YMIn)/YDEC))*NDEC +
  317. & INT(real((XCOOR(IA+2)-ZMIn)/ZDEC))*NMDEC
  318. IPOINT=IZONE/NMULT+1
  319. * WRITE(IOIMP,*) ' point zone ipoint' ,ib,izone,ipoint
  320. np1(ipoint)=np1(ipoint)+1
  321. 3 CONTINUE
  322. endif
  323. C np1= NB DE POINTS PAR SUPER-ZONE DE CLASSEMENT
  324. DO 4 I=2,ite
  325. np1(i)=np1(i-1)+np1(i)
  326. 4 CONTINUE
  327. * WRITE(IOIMP,*) ' np1 avant descente'
  328. * WRITE(IOIMP,*) ( np1( iyou),iyou=1,ite)
  329. if(idim.eq.2) then
  330. do 35 i=1,numnp
  331. ia = ipt2.num(1,i)*idim1-idim
  332. IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+
  333. & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC
  334. IPOINT=IZONE/NMULT+1
  335. J=np1(ipoint)
  336. np3(j)=izone
  337. np2(j)=i
  338. np1(ipoint)=np1(ipoint)-1
  339. 35 continue
  340. else
  341. do 36 i=1,numnp
  342. ia = ipt2.num(1,i)*idim1-idim
  343. IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+
  344. & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC +
  345. & INT(real((XCOOR(IA+2)-ZMIn)/ZDEC))*NMDEC
  346. IPOINT=IZONE/NMULT+1
  347. J=np1(ipoint)
  348. np3(j)=izone
  349. np2(j)=i
  350. np1(ipoint)=np1(ipoint)-1
  351. 36 continue
  352. endif
  353. C
  354. * WRITE(IOIMP,*) ' np1 ite avant permutation' ,ite
  355. * WRITE(IOIMP,*) ( np1( iyuo),iyuo=1,ite)
  356. * WRITE(IOIMP,*) ' np2 ite avant permutation' ,ite
  357. * WRITE(IOIMP,*)( np2( iyuo),iyuo=1,ite)
  358. * WRITE(IOIMP,*) ' np3 ite avant permutation' ,ite
  359. * WRITE(IOIMP,*)( np3( iyuo),iyuo=1,ite)
  360. DO 40 I=1,NUMNP3
  361. JD=NP1(I)+2
  362. JF=NP1(I+1)
  363. * write(6,*) ' i jd jf ' ,i,jd,jf
  364. IF (JD.LE.JF) THEN
  365. 45 IPERM=0
  366. DO 50 J=JD,JF
  367. IF(NP3(J-1).GT.NP3(J)) THEN
  368. L1=NP2(J-1)
  369. NP2(J-1)=NP2(J)
  370. NP2(J)=L1
  371. L1=NP3(J-1)
  372. NP3(J-1)=NP3(J)
  373. NP3(J)=L1
  374. IPERM=1
  375. ENDIF
  376. 50 CONTINUE
  377. IF (IPERM.EQ.1) GO TO 45
  378. ENDIF
  379. 40 CONTINUE
  380.  
  381. * WRITE(IOIMP,*) ' dans zoneg2 apres boucle 40'
  382. C
  383. C a ce stade np1(i)=j le debut de la super zone i est en j-ieme position
  384. C +1
  385. C np2(j)=k le premier point dans la super zone i est le k-ieme
  386. C np3(j)=m sa vrai zone est la m-ieme
  387. C
  388. *
  389. * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT POUR CHAQUE element
  390. * on va boucler sur les points qui sont dans les zones en dessous
  391. *
  392. * WRITE(IOIMP,*) ' np1 ite apres peerm ',ite
  393. * WRITE(IOIMP,*) ( np1( iyuo),iyuo=1,ite)
  394. * WRITE(IOIMP,*) ' np2 ite apres peerm ',ite
  395. * WRITE(IOIMP,*)( np2( iyuo),iyuo=1,ite)
  396. * WRITE(IOIMP,*) ' np3 ite apres peerm ',ite
  397. * WRITE(IOIMP,*)( np3( iyuo),iyuo=1,ite)
  398.  
  399. segact ipt1
  400. xpu(3)=0.d0
  401. do 200 isous=1,ipt1.lisous(/1)
  402. * WRITE(IOIMP,*) ' boucle 200 isous' , isous
  403. meleme = ipt1.lisous(isous)
  404. segact meleme
  405. KEL =itypel
  406. nbnn=num(/1)
  407. C
  408. kd = kdegre(kel)
  409. nso = nbsom(kel)
  410. iad=nspos(kel)-1
  411. segini wrk4
  412. segini inomil
  413. C element quaf ? si oui kell=numero lineaire correspondant
  414. C si non kell=0
  415. CALL NQF2NL(kel,kell)
  416. IF (IERR.NE.0) RETURN
  417. IF (KELL.NE.0) THEN
  418. nbnnl=NBNNE(KELL)
  419. segini wrkl
  420. do i=1,nbnnl
  421. inomil(**)=i
  422. enddo
  423. ELSE
  424. if(kd.eq.3) then
  425. C element quadratiques
  426. do 762 i=1,nbnn
  427. do 763 j=1,nso
  428. iso = ibsom(iad+j)
  429. if(i.eq.iso) goto 762
  430. 763 continue
  431. inomil(**)= i
  432. 762 continue
  433. elseif(kd.eq.2) then
  434. C elements lineaires
  435. do i=1,nbnn
  436. inomil(**)=i
  437. enddo
  438. endif
  439. ENDIF
  440. itest=1
  441. do 101 iel=1,num(/2)
  442. call doxe(xcoor,idim,nbnn,num,iel,xel)
  443. * cas quaf
  444. if (kell.ne.0) then
  445. nsom=nbsom(kel)
  446. if(nsom.ne.nbnnl) then
  447. call erreur(5)
  448. return
  449. endif
  450. idx=nspos(kel)-1
  451. do innl=1,nbnnl
  452. do iid=1,idim
  453. xell(iid,innl)=xel(iid,ibsom(idx+innl))
  454. enddo
  455. enddo
  456. endif
  457.  
  458. iref = num(1,iel)*idim1 - idim
  459. xmi=xcoor(iref)
  460. ymi=xcoor(iref+1)
  461. zmi=xcoor(iref+2)
  462. xma=xmi
  463. yma=ymi
  464. zma=zmi
  465. do 102 ipo=2,num(/1)
  466. iref=num(ipo,iel)*idim1-idim
  467. if(xcoor(iref).lt.xmi) xmi = xcoor(iref)
  468. if(xcoor(iref).gt.xma) xma = xcoor(iref)
  469. if(xcoor(iref+1).lt.ymi) ymi = xcoor(iref+1)
  470. if(xcoor(iref+1).gt.yma) yma = xcoor(iref+1)
  471. if(xcoor(iref+2).lt.zmi) zmi = xcoor(iref+2)
  472. if(xcoor(iref+2).gt.zma) zma = xcoor(iref+2)
  473. 102 continue
  474. * if(iel.lt.2.and.isous.eq.1) then
  475. * WRITE(IOIMP,*) ' element ',iel
  476. * WRITE(IOIMP,*) ( num(iuo,iel),iuo=1,num(/1))
  477. * WRITE(IOIMP,*) ' xmi ymi zmi xma yma zma'
  478. * WRITE(IOIMP,*) ' boite reelle'
  479. * WRITE(IOIMP,*) xmi,ymi,zmi,xma,yma,zma
  480. * endif
  481. xdec1=(xma-xmi)/3.D0
  482. ydec1=(yma-ymi)/3.D0
  483. zdec1=(zma-zmi)/3.D0
  484. xmi=xmi-xdec1
  485. xma=xma+xdec1
  486. ymi=ymi-ydec1
  487. yma=yma+ydec1
  488. zmi=zmi-zdec1
  489. zma=zma+zdec1
  490. if(xmi.lt.xmin) xmi=xmin
  491. if(xmi.gt.xmax) xmi=xmax
  492. if(xma.lt.xmin) xma=xmin
  493. if(xma.gt.xmax) xma=xmax
  494. if(ymi.lt.ymin) ymi=ymin
  495. if(ymi.gt.ymax) ymi=ymax
  496. if(yma.lt.ymin) yma=ymin
  497. if(yma.gt.ymax) yma=ymax
  498. if(zmi.lt.zmin) zmi=zmin
  499. if(zmi.gt.zmax) zmi=zmax
  500. if(zma.lt.zmin) zma=zmin
  501. if(zma.gt.zmax) zma=zmax
  502. * if(iel.lt.2.and.isous.eq.1) then
  503. *
  504. * WRITE(IOIMP,*) ' boite apres modif'
  505. * WRITE(IOIMP,*) xmi,ymi,zmi,xma,yma,zma
  506. * endif
  507. * WRITE(IOIMP,*) ' xmi xma ymi yma zmi zma',xmi,xma,ymi,yma,zmi,zma
  508. * if(idim.eq.2) then
  509. * zmi=0.d0
  510. * zma=0.d0
  511. * endif
  512. C
  513. C
  514. C recherche des points interne à l'element
  515. C pour cela on balaye les super zones
  516. C
  517. if(idim.eq.2) then
  518.  
  519. IXDEP=INT(real((Xmi-XMIN)/XDEC))+1+
  520. & INT(real((ymi-YMIN)/YDEC))*NDEC
  521.  
  522.  
  523. ixpas=INT(real((Xma-XMIN)/XDEC))+1+
  524. & INT(real((ymi-YMIN)/YDEC))*NDEC - ixdep
  525.  
  526. iny=INT(real((yma-YMIn)/YDEC))-INT(real((ymi-YMIn)
  527. $ /YDEC))+1
  528.  
  529. do 120 ily=1,iny
  530. ile=ixdep +( ily -1 )* ndec
  531. ilf=ile + ixpas
  532. izv= ile/nf+1
  533. * if (izv.gt.ite) WRITE(IOIMP,*) ' prob ',izv,ixpas,ixdep,
  534. * > ndec,mdec,xdec,ydec,ymi,ymin
  535. jd=np1(izv)+1
  536. do 122 j=jd,numnp3
  537. izd=np3(j)
  538. if(izd.lt.ile) go to 122
  539. if(izd.gt.ilf) go to 122
  540. ipo= np2(j)
  541. C recherche si le point est interne
  542. if(idejvu(ipo).eq.1) go to 122
  543. ip=ipt2.num(1,ipo)
  544. iref=ip*idim1-idim
  545. xpu(1)=xcoor(iref)
  546. xpu(2)=xcoor(iref+1)
  547. if(xpu(1).lt.xmi.or.xpu(1).gt.xma) go to 122
  548. if(xpu(2).lt.ymi.or.xpu(2).gt.yma) go to 122
  549. * cas quaf
  550. qsi(1) = 0.D0
  551. qsi(2) = 0.D0
  552. qsi(3) = 0.D0
  553. if (kell.ne.0) then
  554. call qsijs(xell,kell,nbnnl,idim,xpu,
  555. $ SHPL,qsi,iret)
  556. call qsijs(xel ,kel ,nbnn ,idim,xpu,
  557. $ SHPP,qsi,iret)
  558. else
  559. call qsijs(xel,kel,nbnn,idim,xpu,SHPP,qsi,iret)
  560. endif
  561. IF (IRET.NE.0) then
  562. irate=irate+1
  563. GOTO 122
  564. endif
  565. IF (KELL.NE.0) THEN
  566. CALL SHAPE(qsi(1),qsi(2),qsi(3),kell,shpl
  567. $ ,iret2)
  568. if (iret2.eq.0) then
  569. MOTERR(1:4)=NOMS(KELL)
  570. CALL ERREUR(68)
  571. RETURN
  572. endif
  573. ENDIF
  574.  
  575. * WRITE(IOIMP,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn)
  576.  
  577. do i=1,inomil(/1)
  578. ilp= inomil(i)
  579. if (kell.ne.0) then
  580. if( shpl(1,ilp).lt.0.D0) goto 110
  581. else
  582. if( shpp(1,ilp).lt.0.D0) goto 110
  583. endif
  584. enddo
  585. C
  586. goto 100
  587. 110 xmap=1.D10
  588. do i=1,inomil(/1)
  589. ilp=inomil(i)
  590. if (kell.ne.0) then
  591. xmap = min ( xmap,shpl(1,ilp))
  592. else
  593. xmap = min ( xmap,shpp(1,ilp))
  594. endif
  595. enddo
  596. if( xmap . gt . ckrit(ipo) ) then
  597. do 180 jl =1,NUM(/1)
  598. if(IP.NE.NUM(jl,iel)) go to 180
  599. if(IFL.EQ.0) go to 190
  600. idejvu(ipo)=1
  601. itr=itr+1
  602. go to 181
  603. 180 continue
  604. 181 continue
  605. ckrit(ipo)=xmap
  606. ieint(1,ipo)=isous
  607. ieint(2,ipo)=iel
  608. QQQ(1,ipo) = qsi(1)
  609. QQQ(2,ipo) = qsi(2)
  610. endif
  611. 190 continue
  612. go to 122
  613. 100 CONTINUE
  614.  
  615. C on a trouve l element contenant le point ipo attention
  616. C si le point appartient aux deux maillages on ne le stocke pas
  617. C pour accro mais on le stocke pour proi
  618. do 18 jl =1,NUM(/1)
  619. if((IP.EQ.NUM(jl,iel)).AND.(IFL.EQ.0)) GOTO
  620. $ 19
  621. 18 continue
  622. if(idejvu(ipo).EQ.0) then
  623. idejvu(ipo) = 1
  624. itr = itr + 1
  625. IEINT(1,ipo) = isous
  626. IEINT(2,ipo) = iel
  627. QQQ(1,ipo) = qsi(1)
  628. QQQ(2,ipo) = qsi(2)
  629. * WRITE(IOIMP,*) 'trouve point' ,IP, 'Element' ,J
  630. * WRITE(IOIMP,2375)ip,j,(xpu(i1),i1=1,idim),(qsi(i2),i2=1,idim)
  631. endif
  632. 19 continue
  633. * 2375 format(2i4,6e12.5)
  634. 122 continue
  635. 120 continue
  636. endif
  637. if(idim.eq.3) then
  638.  
  639.  
  640.  
  641. IXDEP=INT(real((Xmi-XMIN)/XDEC))+1+
  642. & INT(real((ymi-YMIN)/YDEC))*NDEC +
  643. & INT(real((zmi-zmin)/ZDEC))*NMDEC
  644.  
  645.  
  646. ixpas=INT(real((Xma-XMIN)/XDEC))+1+
  647. & INT(real((ymi-YMIN)/YDEC))*NDEC +
  648. & INT(real((zmi-zmin)/ZDEC))*NMDEC - ixdep
  649.  
  650.  
  651. iny=INT(real((yma-YMIn)/YDEC))-INT(real((ymi-YMIn)
  652. $ /YDEC))+1
  653. inz= INT(real((zma-zmin)/ZDEC))-INT(real((zmi-zmin)
  654. $ /ZDEC))+1
  655.  
  656. * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) ' inz iny',inz,iny
  657. do 240 itz=1,inz
  658. idexx= ixdep + (itz-1)*nmdec
  659. do 220 ily=1,iny
  660. ile=idexx +( ily -1 )* ndec
  661. ilf=ile + ixpas
  662. * if(iel.lt.2.and.isous.eq.1)WRITE(IOIMP,*)' ile,ilf' , ile,ilf
  663. izv= ile/nf + 1
  664. jd=np1(izv)+1
  665. * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) ' jd,jf ',jd,jf
  666. do 222 j=jd,numnp3
  667. izd=np3(j)
  668. * WRITE(IOIMP,*) ' izd ', izd
  669. if(izd.lt.ile) go to 222
  670. if(izd.gt.ilf) go to 223
  671. ipo= np2(j)
  672. C recherche si le point est interne
  673. if(idejvu(ipo).eq.1) go to 222
  674. ip=ipt2.num(1,ipo)
  675. * if(iel.lt.2.and.isous.eq.1)WRITE(IOIMP,*) ' ipo ip ',ipo,ip
  676. iref=ip*idim1-idim
  677. xpu(1)=xcoor(iref)
  678. xpu(2)=xcoor(iref+1)
  679. xpu(3)=xcoor(iref+2)
  680.  
  681. * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) 'cooor '
  682. * $ ,xpu(1),xpu(2),xpu(3)
  683. if(xpu(1).lt.xmi.or.xpu(1).gt.xma) go to 222
  684. if(xpu(2).lt.ymi.or.xpu(2).gt.yma) go to 222
  685. if(xpu(3).lt.zmi.or.xpu(3).gt.zma) go to 222
  686. qsi(1) = 0.D0
  687. qsi(2) = 0.D0
  688. qsi(3) = 0.D0
  689. call qsijs(xel,kel,nbnn,idim,xpu,SHPP,qsi,iret)
  690. IF (IRET.NE.0) THEN
  691. irate=irate+1
  692. if(irate.eq.itest) then
  693. * WRITE(IOIMP,*) 'solution pas trouvee', irate
  694. itest=itest *2
  695. endif
  696. GOTO 222
  697. endif
  698. IF (KELL.NE.0) THEN
  699. CALL SHAPE(qsi(1),qsi(2),qsi(3),kell,shpl
  700. $ ,iret2)
  701. ENDIF
  702. C
  703. * WRITE(IOIMP,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn)
  704. do i=1,inomil(/1)
  705. ilp= inomil(i)
  706. if (kell.ne.0) then
  707. if( shpl(1,ilp).lt.0.D0) goto 210
  708. else
  709. if( shpp(1,ilp).lt.0.D0) goto 210
  710. endif
  711. enddo
  712. C
  713. goto 300
  714. 210 xmap=1.D10
  715. do i=1,inomil(/1)
  716. ilp=inomil(i)
  717. if (kell.ne.0) then
  718. xmap = min ( xmap,shpl(1,ilp))
  719. else
  720. xmap = min ( xmap,shpp(1,ilp))
  721. endif
  722. enddo
  723. if( xmap . gt . ckrit(ipo) ) then
  724. do 280 jl =1,NUM(/1)
  725. if(IP.NE.NUM(jl,iel)) go to 280
  726. if(IFL.EQ.0) go to 290
  727. idejvu(ipo)=1
  728. itr=itr+1
  729. * if(iel.lt.3.and.isous.eq.1) WRITE(IOIMP,*)'egalite de point'
  730. go to 281
  731. 280 continue
  732. 281 continue
  733. ckrit(ipo)=xmap
  734. ieint(1,ipo)=isous
  735. ieint(2,ipo)=iel
  736. QQQ(1,ipo) = qsi(1)
  737. QQQ(2,ipo) = qsi(2)
  738. QQQ(3,ipo) = qsi(3)
  739. endif
  740. 290 continue
  741. go to 222
  742. 300 CONTINUE
  743.  
  744. C on a trouve l element contenant le point ipo attention
  745. C si le point appartient aux deux maillages on ne le stocke pas
  746. C pour accro mais on le stocke pour proi
  747. do 28 jl =1,NUM(/1)
  748. if((IP.EQ.NUM(jl,iel)).AND.(IFL.EQ.0))
  749. $ GOTO 29
  750. 28 continue
  751. if(idejvu(ipo).EQ.0) then
  752.  
  753. idejvu(ipo) = 1
  754. * if(iel.lt.3.and.isous.eq.1) WRITE(IOIMP,*)'point interne'
  755. itr = itr + 1
  756. IEINT(1,ipo) = isous
  757. IEINT(2,ipo) = iel
  758. QQQ(1,ipo) = qsi(1)
  759. QQQ(2,ipo) = qsi(2)
  760. QQQ(3,ipo) = qsi(3)
  761. * if(iel.lt.3.and.isous.eq.1) then
  762. * WRITE(IOIMP,*) 'trouve point' ,IP, 'Element' ,iel
  763. * WRITE(IOIMP,2375)ip,iel,(xpu(i1),i1=1,idim),(qsi(i2),i2=1,idim
  764. * endif
  765. endif
  766. 29 continue
  767. 222 continue
  768. 223 continue
  769. 220 continue
  770. 240 continue
  771. endif
  772. 101 continue
  773. segsup inomil,wrk4
  774. IF (KELL.NE.0) THEN
  775. SEGSUP wrkl
  776. ENDIF
  777. segdes meleme
  778. 200 continue
  779. segdes ipt1
  780. GOTO 800
  781.  
  782. C==========================
  783. C= CAS 1D - DIMENSION 1 =
  784. C==========================
  785. 600 CONTINUE
  786.  
  787. C Recherche des min et max
  788. XMIN=XGRAND
  789. XMAX=-XGRAND
  790. DO i=1,nodes
  791. j=IPT2.NUM(1,i)*IDIM1-IDIM
  792. xx=XCOOR(j)
  793. XMIN=MIN(xx,XMIN)
  794. XMAX=MAX(xx,XMAX)
  795. ENDDO
  796. XDIST=XMAX-XMIN
  797. xdistv=XDIST
  798. IF (XDIST.EQ.0.D0) XDIST=1.D0
  799.  
  800. C Recherche sur les elements de la plus petite taille en x
  801. SEGACT,ipt1
  802. nbsous=ipt1.lisous(/1)
  803. dminx=XGRAND
  804. DO isous=1,nbsous
  805. meleme=ipt1.lisous(isous)
  806. SEGACT,meleme
  807. nbnn=num(/1)
  808. DO iel=1,num(/2)
  809. dmix=XGRAND
  810. dmax=-XGRAND
  811. DO i=1,nbnn
  812. j=num(i,iel)*IDIM1-IDIM
  813. xx=XCOOR(j)
  814. dmix=MIN(dmix,xx)
  815. dmax=MAX(dmax,xx)
  816. ENDDO
  817. xx=dmax-dmix
  818. dminx=MIN(dminx,xx)
  819. ENDDO
  820. ENDDO
  821. C- WRITE(IOIMP,*) 'dminx',dminx,' xdist',xdist
  822. IF (dminx.EQ.0.D0) THEN
  823. IF (xdistv.EQ.0.D0) THEN
  824. dminx=1.D0
  825. ELSE
  826. dminx=xdistv/1000.D0
  827. ENDIF
  828. ENDIF
  829. dminx=MAX(dminx,xdist/10000.D0)
  830.  
  831. ndex=xdist/dminx
  832. ndex=MAX(1,ndex)
  833. C- WRITE(IOIMP,*) ' ndex=',ndex
  834.  
  835. C Calcul du decoupage en zones
  836. nbzone=NUMNP*NF
  837. rap=REAL(nbzone)/ndex
  838. rap=MIN(Un,rap)
  839. ndec=ndex*rap
  840. ndec = max(1,ndec)
  841. nbzone=ndec
  842. XDEC=xdist/ndec*1.0000001D0
  843. nf=(nbzone/numnp3)+1
  844. nmult=nf
  845. C- WRITE(IOIMP,601) ndec,nbzone
  846. C- 601 FORMAT(' DECOUPAGE EN X ',I6,' SOIT ',I8,' ZONES ')
  847. C- WRITE(IOIMP,602) XMIN,XMAX,XDEC
  848. C- 602 FORMAT(' XMIN,XMAX,XDEC ',3G12.5)
  849.  
  850. DO i=1,numnp
  851. j=IPT2.num(1,i)*IDIM1-IDIM
  852. izone=INT((XCOOR(j)-XMIN)/XDEC)+1
  853. ipoint=izone/nmult+1
  854. np1(ipoint)=np1(ipoint)+1
  855. C- WRITE(IOIMP,*) 'point,izone,ipoint =',i,j,izone,ipoint
  856. ENDDO
  857. C np1 = nb de points par super-zone de classement
  858. DO i=2,NUMNP3
  859. np1(i)=np1(i)+np1(i-1)
  860. ENDDO
  861. C- WRITE(IOIMP,*) ' np1 avant descente'
  862. C- WRITE(IOIMP,*) (np1(i),i=1,ite)
  863. DO i=1,numnp
  864. j=IPT2.num(1,i)*IDIM1-IDIM
  865. izone=INT((XCOOR(j)-XMIN)/XDEC)+1
  866. ipoint=izone/nmult+1
  867. j=np1(ipoint)
  868. np3(j)=izone
  869. np2(j)=i
  870. np1(ipoint)=j-1
  871. ENDDO
  872.  
  873. C- WRITE(IOIMP,*) ' np1 ite avant permutation' ,ite
  874. C- WRITE(IOIMP,*) (np1(i),i=1,ite)
  875. C- WRITE(IOIMP,*) (np2(i),i=1,ite)
  876. C- WRITE(IOIMP,*) (np3(i),i=1,ite)
  877. DO i=1,NUMNP
  878. jD=np1(i)+2
  879. jF=np1(i+1)
  880. IF (jD.LE.jF) THEN
  881. 610 IPERM=0
  882. DO j=jD,jF
  883. k=j-1
  884. IF (np3(k).GT.np3(j)) THEN
  885. L1=np2(k)
  886. np2(k)=np2(j)
  887. np2(j)=L1
  888. L1=np3(k)
  889. np3(k)=np3(j)
  890. np3(j)=L1
  891. IPERM=1
  892. ENDIF
  893. ENDDO
  894. IF (IPERM.EQ.1) GOTO 610
  895. ENDIF
  896. ENDDO
  897. C- WRITE(IOIMP,*) ' dans zoneg2 apres permutation'
  898. C- WRITE(IOIMP,*) ' np1 ite apres perm ',ite
  899. C- WRITE(IOIMP,*) (np1(i),i=1,ite)
  900. C- WRITE(IOIMP,*) (np2(i),i=1,ite)
  901. C- WRITE(IOIMP,*) (np3(i),i=1,ite)
  902.  
  903. C A ce stade :
  904. C - np1(i)=j le debut de la super-zone i est en j-eme position + 1
  905. C - np2(j)=k le premier point dans la super zone i est le k-eme
  906. C - np3(j)=m sa vrai zone est la m-eme
  907. C Il ne reste plus qu'a faire le travail pour chaque element.
  908. C On va boucler sur les points qui sont dans les zones en dessous.
  909.  
  910. xpu(2)=0.D0
  911. xpu(3)=0.D0
  912. C SEGACT,ipt1
  913. DO isous=1,nbsous
  914. C- WRITE(IOIMP,*) ' boucle isous=',isous
  915. meleme=ipt1.lisous(isous)
  916. C SEGACT,meleme
  917. nbnn=num(/1)
  918.  
  919. SEGINI,wrk4
  920. SEGINI,inomil
  921.  
  922. kel=ITYPEL
  923. kd=KDEGRE(kel)
  924. C Elements quadratiques
  925. IF (kd.EQ.3) then
  926. jd=NSPOS(kel)
  927. jf=jd+NBSOM(kel)-1
  928. DO i=1,nbnn
  929. DO j=jd,jf
  930. IF (i.EQ.ibsom(j)) GOTO 620
  931. ENDDO
  932. inomil(**)=i
  933. 620 CONTINUE
  934. ENDDO
  935. C Elements lineaires
  936. ELSE IF (kd.EQ.2) THEN
  937. DO i=1,nbnn
  938. inomil(**)=i
  939. ENDDO
  940. ENDIF
  941.  
  942. DO iel=1,num(/2)
  943. CALL DOXE(XCOOR,IDIM,nbnn,num,iel,xel)
  944. j=num(1,iel)*IDIM1-IDIM
  945. xmi=XCOOR(j)
  946. xma=xmi
  947. DO i=2,nbnn
  948. j=num(i,iel)*IDIM1-IDIM
  949. xx=XCOOR(j)
  950. xmi=MIN(xmi,xx)
  951. xma=MAX(xma,xx)
  952. ENDDO
  953. C- IF (iel.EQ.1.AND.isous.EQ.1) THEN
  954. C- WRITE(IOIMP,*) ' element ',iel
  955. C- WRITE(IOIMP,*) (num(i,iel),i=1,nbnn)
  956. C- WRITE(IOIMP,*) ' boite reelle : xmi xma ',xmi,xma
  957. C- ENDIF
  958. xdec1=us3*(xma-xmi)
  959. xmi=xmi-xdec1
  960. xma=xma+xdec1
  961. xmi=MAX(xmi,xmin)
  962. xmi=MIN(xmax,xmi)
  963. xma=MAX(xma,xmin)
  964. xma=MIN(xma,xmax)
  965. C- IF (iel.EQ.1.AND.isous.EQ.1) THEN
  966. C- WRITE(IOIMP,*) ' element ',iel
  967. C- WRITE(IOIMP,*) ' boite apres modif : xmi xma ',xmi,xma
  968. C- ENDIF
  969. C- WRITE(IOIMP,*) ' xmi xma :',xmi,xma
  970.  
  971. C Recherche des points internes a l'element par balayage des super zones
  972. ixdep=INT((xmi-xmin)/xdec)+1
  973. ixpas=INT((xma-xmin)/xdec)+1
  974. ile=ixdep
  975. ilf=ile+ixpas
  976. izv=ile/nf+1
  977. C- IF (izv.GT.ite) WRITE(IOIMP,*) ' prob ',izv,ixpas,ixdep,ndec,xdec
  978. jd=np1(izv)+1
  979. DO 630 j=jd,numnp3
  980. izd=np3(j)
  981. IF (izd.LT.ile.OR.izd.GT.ilf) GOTO 630
  982. C Recherche si le point est interne
  983. ipo=np2(j)
  984. IF (idejvu(ipo).EQ.1) GOTO 630
  985. ip=ipt2.num(1,ipo)
  986. xpu(1)=XCOOR(ip*IDIM1-IDIM)
  987. IF (xpu(1).LT.xmi.OR.xpu(1).GT.xma) GOTO 630
  988. qsi(1) = 0.D0
  989. qsi(2) = 0.D0
  990. qsi(3) = 0.D0
  991. CALL QSIJS(xel,kel,nbnn,IDIM,xpu,shpp,qsi,IRET)
  992. C- WRITE(IOIMP,FMT='(8E12.5)') (shpp(1,i),i=1,nbnn)
  993. IF (IRET.NE.0) then
  994. irate=irate+1
  995. GOTO 630
  996. ENDIF
  997.  
  998. do i=1,inomil(/1)
  999. ilp= inomil(i)
  1000. if( shpp(1,ilp).lt.0.D0) goto 1110
  1001. enddo
  1002.  
  1003. goto 1100
  1004. 1110 xmap=1.D10
  1005. do i=1,inomil(/1)
  1006. ilp=inomil(i)
  1007. xmap = min ( xmap,shpp(1,ilp))
  1008. enddo
  1009. if( xmap . gt . ckrit(ipo) ) then
  1010. do 1180 jl =1,NUM(/1)
  1011. if(IP.NE.NUM(jl,iel)) go to 1180
  1012. if(IFL.EQ.0) go to 1190
  1013. idejvu(ipo)=1
  1014. itr=itr+1
  1015. go to 1181
  1016. 1180 continue
  1017. 1181 continue
  1018. ckrit(ipo)=xmap
  1019. ieint(1,ipo)=isous
  1020. ieint(2,ipo)=iel
  1021. QQQ(1,ipo) = qsi(1)
  1022. endif
  1023. 1190 continue
  1024. go to 630
  1025. 1100 CONTINUE
  1026.  
  1027. C on a trouve l element contenant le point ipo attention
  1028. C si le point appartient aux deux maillages on ne le stocke pas
  1029. C pour accro mais on le stocke pour proi
  1030. DO i=1,nbnn
  1031. IF ((ip.EQ.NUM(i,iel)).AND.(IFL.EQ.0)) GOTO 630
  1032. ENDDO
  1033. IF (idejvu(ipo).EQ.0) THEN
  1034. idejvu(ipo)=1
  1035. itr=itr+1
  1036. IEINT(1,ipo)=isous
  1037. IEINT(2,ipo)=iel
  1038. QQQ(1,ipo)=qsi(1)
  1039. C WRITE(IOIMP,*) 'Trouve point de numero global ',IP
  1040. C WRITE(IOIMP,*) 'dans le ',iel,'eme élément ',
  1041. C $ 'de la ',isous,'eme sous-zone'
  1042. C WRITE(IOIMP,*) 'Coordonnée globale = ',xpu(1)
  1043. C WRITE(IOIMP,*) 'Coordonnée locale = ',qsi(1)
  1044. CC- WRITE(IOIMP,*) 'Trouve point' ,IP, 'Element' ,J
  1045. CC- WRITE(IOIMP,603) ip,j,xpu(1),qsi(1)
  1046. C 603 FORMAT(2I4,2E12.5)
  1047. ENDIF
  1048. 630 CONTINUE
  1049. ENDDO
  1050. SEGSUP,inomil,wrk4
  1051. SEGDES,meleme
  1052. ENDDO
  1053. SEGDES IPT1
  1054. 800 CONTINUE
  1055. IACCRO=IACCRO+itr
  1056. C- WRITE(IOIMP,*) ' nb points non converges dans qsijs ',irate
  1057. SEGDES,IPT2
  1058. SEGSUP STRAV
  1059.  
  1060. RETURN
  1061. END
  1062.  
  1063.  
  1064.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales