Télécharger zoneg2.eso

Retour à la liste

Numérotation des lignes :

zoneg2
  1. C ZONEG2 SOURCE CB215821 22/07/20 15:39:47 11411
  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 CCREEL
  22.  
  23. -INC PPARAM
  24. -INC CCOPTIO
  25. -INC SMCOORD
  26. -INC SMELEME
  27. -INC CCGEOME
  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.  
  43. SEGMENT ICOUNT
  44. INTEGER IEINT(2,NODES),IDEJVU(NODES),IACCRO
  45. REAL*8 QQQ(3,NODES),CKRIT(NODES)
  46. ENDSEGMENT
  47. EXTERNAL SHAPE
  48. REAL*8 ZDIST
  49.  
  50. ZDIST=0.D0
  51. ip=1
  52. ipo=1
  53. * WRITE(IOIMP,*) ' entree dans zoneg2'
  54. itr=0
  55. irate=0
  56. IDIM1 = IDIM + 1
  57. IPT2 = IP2
  58. ipt1=ip1
  59.  
  60. * on fait le zonage sur les points
  61. SEGACT IPT2
  62. Nodes = IPT2.NUM(/2)
  63. * WRITE(IOIMP,*) ' nodes ' , nodes
  64.  
  65. numnp=nodes
  66. ITE=numnp+3+1
  67. numnp3=ite-1
  68. * write(6,*) ' numnp ite numnp3 ' , numnp , ite , numnp3
  69. NF=100
  70. NMULT=100
  71. segini strav
  72. C
  73. C Cas 1D : idem qu'en 2D ou 3D mais pour eviter multiplications
  74. C -------- des tests (IF) qui alourdissent, on duplique ZONEG2
  75. C pour le 1D en 600
  76. IF (IDIM.EQ.1) GOTO 600
  77.  
  78. C recherche des min et max
  79. C
  80. iref = ipt2.NUM(1,1)*idim1-idim
  81. XMIN=xcoor(iref)
  82. YMIN=xcoor(iref+1)
  83. ZMIN=xcoor(iref+2)
  84. YMAx=ymin
  85. XMAx=xmin
  86. ZMAx=zmin
  87. DO 201 Iyr = 2,nodes
  88. IB=ipt2.NUM(1,Iyr)
  89. IA=IB*IDIM1-idim
  90. IF(XCOOR(IA).LT.XMIN) XMIn=XCOOR(IA)
  91. IF(XCOOR(IA).GT.XMAx) XMAx=XCOOR(IA)
  92. IF(XCOOR(IA+1).LT.YMIn) YMIn=XCOOR(IA+1)
  93. IF(XCOOR(IA+1).GT.YMAx) YMAx=XCOOR(IA+1)
  94. IF( IDIM.EQ.3 ) THEN
  95. IF(XCOOR(IA+2).LT.ZMIn) ZMIn=XCOOR(IA+2)
  96. IF(XCOOR(IA+2).GT.ZMAx) ZMAx=XCOOR(IA+2)
  97. ENDIF
  98. 201 CONTINUE
  99. if(idim.eq.2) then
  100. ZMIN=0.D0
  101. ZMAX=0.D0
  102. ZDIST=1.D0
  103. ZDISTV=1.D0
  104. endif
  105. XDIST=XMAx-XMIn
  106. YDIST=YMAx-YMIn
  107. xdistv=xdist
  108. ydistv=ydist
  109. IF(XDIST.EQ.0.D0) THEN
  110. NDEX=1
  111. XDIST=1.D0
  112. ELSE
  113. NDEX=0
  114. ENDIF
  115. IF(YDIST.EQ.0.D0) THEN
  116. NDEY=1
  117. YDIST=1.D0
  118. ELSE
  119. NDEY=0
  120. ENDIF
  121. * SMO=1E-5*XDIST*YDIST/NUMNP
  122. IF( IDIM.EQ.3 ) then
  123. zdist = zmax-zmin
  124. zdistv=zdist
  125. IF(ZDIST.EQ.0.D0) THEN
  126. NDEZ=1
  127. ZDIST=1.D0
  128. * SMO=1E-5*XDIST*YDIST*ZDIST/NUMNP
  129. ELSE
  130. NDEZ=0
  131. ENDIF
  132. endif
  133. * write(ioimp,*) 'Maillage POI1'
  134. * write(ioimp,*) 'xmin,xmax,ymin,ymax,zmin,zmax=',xmin,xmax,ymin
  135. * $ ,ymax,zmin,zmax
  136. * write(ioimp,*) 'xdistv,ydistv,zdistv=',xdistv,ydistv,zdistv
  137. * write(ioimp,*) 'xdist,ydist,zdist=',xdist,ydist,zdist
  138.  
  139. C
  140. C*** recherche sur les elements de la plus petite taille en x y et z
  141. C
  142. ipt1=ip1
  143. segact ipt1
  144. dminx=1.D10
  145. dminy=1.D10
  146. dminz=1.D10
  147. do 453 isous=1,ipt1.lisous(/1)
  148. meleme=ipt1.lisous(isous)
  149. segact meleme
  150. do 454 iel=1,num(/2)
  151. dmix=1.D10
  152. dmax=-1.D10
  153. dmiy=1.D10
  154. dmay=-1.D10
  155. dmiz=1.D10
  156. dmaz=-1.D10
  157. do 455 ip=1,num(/1)
  158. ire=num(ip,iel)*idim1-idim
  159. if(xcoor(ire).le.dmix) dmix=xcoor(ire)
  160. if(xcoor(ire).gt.dmax) dmax=xcoor(ire)
  161. if(xcoor(ire+1).le.dmiy) dmiy=xcoor(ire+1)
  162. if(xcoor(ire+1).gt.dmay) dmay=xcoor(ire+1)
  163. if(xcoor(ire+2).le.dmiz) dmiz=xcoor(ire+2)
  164. if(xcoor(ire+2).gt.dmaz) dmaz=xcoor(ire+2)
  165. 455 continue
  166. if(dminx.gt.(dmax-dmix) ) dminx= dmax-dmix
  167. if(dminy.gt.(dmay-dmiy) ) dminy= dmay-dmiy
  168. if(dminz.gt.(dmaz-dmiz) ) dminz= dmaz-dmiz
  169. 454 continue
  170. 453 continue
  171. * write(ioimp,*) 'Maillage Massif'
  172. * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz
  173. * WRITE(IOIMP,*) ' xdist,ydist,zdist',xdist,ydist,zdist
  174.  
  175. if( dminx.eq.0.D0) then
  176. if(xdistv.eq.0.D0) then
  177. dminx=1.D0
  178. else
  179. dminx = xdistv/1000.D0
  180. endif
  181. endif
  182. if( dminy.eq.0.D0) then
  183. if(ydistv.eq.0.D0) then
  184. dminy=1.D0
  185. else
  186. dminy = ydistv/1000.D0
  187. endif
  188. endif
  189. if (idim.eq.3) then
  190. if( dminz.eq.0.D0) then
  191. if(zdistv.eq.0.D0) then
  192. dminz=1.D0
  193. else
  194. dminz = zdistv/1000.D0
  195. endif
  196. endif
  197. endif
  198. *
  199. * write(ioimp,*) 'Apres correction 1'
  200. * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz
  201.  
  202. dminx = max( dminx , xdist/10000.D0)
  203. dminy = max(dminy,ydist/10000.D0)
  204. if( idim.eq.3) dminz = max(dminz,zdist/10000.D0)
  205.  
  206. * write(ioimp,*) 'Apres correction 2'
  207. * WRITE(IOIMP,*) 'dminx,dminy,dminz=',dminx,dminy,dminz
  208.  
  209. if (ndex.eq.0) ndex=xdist/dminx
  210. ndex=max(1,ndex)
  211. if (ndey.eq.0) ndey=ydist/dminy
  212. ndey=max(1,ndey)
  213. if(idim.eq.3) then
  214. if (ndez.eq.0) ndez=zdist/dminz
  215. ndez=max(1,ndez)
  216. else
  217. ndez=1
  218. zdist=1.D0
  219. dminz=1.01D0
  220. endif
  221. * WRITE(IOIMP,*) ' ndex ndey ndez',ndex,ndey,ndez
  222.  
  223. C CALCUL DU DECOUPAGE EN ZONE
  224. NBZONE=NUMNP*NF
  225. * if( ndex*ndey*ndez.lt.nbzone) then
  226. * WRITE(IOIMP,*) ' ndex ndey ndez nbzone=',ndex,ndey,ndez,nbzone
  227. xnprop = real(ndex)*real(ndey)*real(ndez)
  228. * dimension effective : là où les découpages sont plus grands que 1
  229. idimf=0
  230. if (ndex.gt.1) idimf=idimf+1
  231. if (ndey.gt.1) idimf=idimf+1
  232. if (ndez.gt.1) idimf=idimf+1
  233. * write(ioimp,*) ' idimf=',idimf
  234. * if(idim.eq.3) then
  235. * rap = (real(nbzone)/xnprop)**0.333333333D0
  236. * else
  237. * rap= (real(nbzone)/xnprop)** 0.5D0
  238. * endif
  239. if (idimf.gt.0) then
  240. rap = (real(nbzone)/xnprop)**(1.D0/real(idimf))
  241. * write(ioimp,*) ' rap=',rap
  242. if(rap.gt.1.D0) rap=1.d0
  243. else
  244. rap=1.d0
  245. endif
  246. * write(ioimp,*) ' rap=',rap
  247. if (ndex.gt.1) then
  248. ndec=ndex*rap
  249. ndec=max(1,ndec)
  250. else
  251. ndec=1
  252. endif
  253. *
  254. if (ndey.gt.1) then
  255. mdec=ndey*rap
  256. mdec = max(1,mdec)
  257. else
  258. mdec=1
  259. endif
  260. *
  261. if(idim.eq.3) then
  262. if (ndez.gt.1) then
  263. ldec=ndez*rap
  264. ldec = max(1,ldec)
  265. else
  266. ldec=1
  267. endif
  268. else
  269. ldec=1
  270. endif
  271. *
  272. nbzone=ndec*mdec*ldec
  273. nf= nbzone / numnp3 +1
  274. nmult=nf
  275. * write(ioimp,*) 'Apres correction 3'
  276. * WRITE(IOIMP,*) ' ndec mdec ldec nbzone=',ndec,mdec,ldec,nbzone
  277. * else
  278. * RAPY=YDIST/XDIST
  279. * NDEC=MAX(INT(SQRT(REAL(NBZONE)/RAPY)),1)
  280. * MDEC=NBZONE/NDEC
  281. * LDEC=1
  282. * IF(IDIM.EQ.3) THEN
  283. * NDEC=max(int((real(nbzone)*xdist*xdist/(ydist*zdist))
  284. * $ **0.3333333D0),1)
  285. * MDEC= max(int((real(nbzone)*ydist*ydist/(xdist*zdist))
  286. * $ **0.3333333D0),1)
  287. * LDEC=NBZONE/(MDEC*NDEC)
  288. * ENDIF
  289. NBZONE=NDEC*MDEC*LDEC
  290. * endif
  291. * WRITE (IOIMP,8933) NDEC,MDEC,ldec,NBZONE,nmult
  292. * 8933 FORMAT(' DECOUPAGE EN X ',I6,' EN Y ',I6,' EN Z ',i6,
  293. * $ ' SOIT ',I8,' ZONES ')
  294. * write(6,*) 'diviseur des zones pour aller aux super zones',nmult
  295. * write(6,*) ' nombre max de super zones ', nbzone/nmult+1
  296. XDEC=XDIST/NDEC*1.0000001D0
  297. YDEC=YDIST/MDEC*1.0000001D0
  298. ZDEC=ZDIST/LDEC*1.0000001D0
  299. * WRITE (IOIMP,8935) XMIN,XMAX,YMIN,YMAX,zmin, zmax,XDEC,YDEC,zdec
  300. * 8935 FORMAT (' XMIN,XMAX,YMIN,YMAX,XDEC,YDEC ',2(6G12.5))
  301. NMDEC=NDEC*MDEC
  302. if(idim.eq.2) THEN
  303. DO 2 I=1,NUMNP
  304. IB=IPT2.num(1,i)
  305. IA=IB*IDIM1-idim
  306. IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+
  307. & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC
  308. IPOINT=IZONE/NMULT+1
  309. * WRITE(IOIMP,*) ' point zone ipoint' ,ib,izone,ipoint
  310. np1(ipoint)=np1(ipoint)+1
  311. 2 CONTINUE
  312. ELSE
  313. DO 3 I=1,NUMNP
  314. IB=IPT2.num(1,i)
  315. IA=IB*IDIM1 -idim
  316. IZONE=INT(real((XCOOR(IA)-XMIn)/XDEC))+1+
  317. & INT(real((XCOOR(IA+1)-YMIn)/YDEC))*NDEC +
  318. & INT(real((XCOOR(IA+2)-ZMIn)/ZDEC))*NMDEC
  319. IPOINT=IZONE/NMULT+1
  320. * WRITE(IOIMP,*) ' point zone ipoint' ,ib,izone,ipoint
  321. np1(ipoint)=np1(ipoint)+1
  322. 3 CONTINUE
  323. endif
  324. C np1= NB DE POINTS PAR SUPER-ZONE DE CLASSEMENT
  325. DO 4 I=2,ite
  326. np1(i)=np1(i-1)+np1(i)
  327. 4 CONTINUE
  328. * WRITE(IOIMP,*) ' np1 avant descente'
  329. * WRITE(IOIMP,*) ( np1( iyou),iyou=1,ite)
  330. if(idim.eq.2) then
  331. do 35 i=1,numnp
  332. ia = ipt2.num(1,i)*idim1-idim
  333. IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+
  334. & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC
  335. IPOINT=IZONE/NMULT+1
  336. J=np1(ipoint)
  337. np3(j)=izone
  338. np2(j)=i
  339. np1(ipoint)=np1(ipoint)-1
  340. 35 continue
  341. else
  342. do 36 i=1,numnp
  343. ia = ipt2.num(1,i)*idim1-idim
  344. IZONE=INT(real((XCOOR(IA)-XMIN)/XDEC))+1+
  345. & INT(real((XCOOR(IA+1)-YMIN)/YDEC))*NDEC +
  346. & INT(real((XCOOR(IA+2)-ZMIn)/ZDEC))*NMDEC
  347. IPOINT=IZONE/NMULT+1
  348. J=np1(ipoint)
  349. np3(j)=izone
  350. np2(j)=i
  351. np1(ipoint)=np1(ipoint)-1
  352. 36 continue
  353. endif
  354. C
  355. * WRITE(IOIMP,*) ' np1 ite avant permutation' ,ite
  356. * WRITE(IOIMP,*) ( np1( iyuo),iyuo=1,ite)
  357. * WRITE(IOIMP,*) ' np2 ite avant permutation' ,ite
  358. * WRITE(IOIMP,*)( np2( iyuo),iyuo=1,ite)
  359. * WRITE(IOIMP,*) ' np3 ite avant permutation' ,ite
  360. * WRITE(IOIMP,*)( np3( iyuo),iyuo=1,ite)
  361. DO 40 I=1,NUMNP3
  362. JD=NP1(I)+2
  363. JF=NP1(I+1)
  364. * write(6,*) ' i jd jf ' ,i,jd,jf
  365. IF (JD.LE.JF) THEN
  366. 45 IPERM=0
  367. DO 50 J=JD,JF
  368. IF(NP3(J-1).GT.NP3(J)) THEN
  369. L1=NP2(J-1)
  370. NP2(J-1)=NP2(J)
  371. NP2(J)=L1
  372. L1=NP3(J-1)
  373. NP3(J-1)=NP3(J)
  374. NP3(J)=L1
  375. IPERM=1
  376. ENDIF
  377. 50 CONTINUE
  378. IF (IPERM.EQ.1) GO TO 45
  379. ENDIF
  380. 40 CONTINUE
  381.  
  382. * WRITE(IOIMP,*) ' dans zoneg2 apres boucle 40'
  383. C
  384. C a ce stade np1(i)=j le debut de la super zone i est en j-ieme position
  385. C +1
  386. C np2(j)=k le premier point dans la super zone i est le k-ieme
  387. C np3(j)=m sa vrai zone est la m-ieme
  388. C
  389. *
  390. * IL NE RESTE PLUS QU'A FAIRE LE TRAVAIL PROPREMENT DIT POUR CHAQUE element
  391. * on va boucler sur les points qui sont dans les zones en dessous
  392. *
  393. * WRITE(IOIMP,*) ' np1 ite apres peerm ',ite
  394. * WRITE(IOIMP,*) ( np1( iyuo),iyuo=1,ite)
  395. * WRITE(IOIMP,*) ' np2 ite apres peerm ',ite
  396. * WRITE(IOIMP,*)( np2( iyuo),iyuo=1,ite)
  397. * WRITE(IOIMP,*) ' np3 ite apres peerm ',ite
  398. * WRITE(IOIMP,*)( np3( iyuo),iyuo=1,ite)
  399.  
  400. segact ipt1
  401. xpu(3)=0.d0
  402. do 200 isous=1,ipt1.lisous(/1)
  403. * WRITE(IOIMP,*) ' boucle 200 isous' , isous
  404. meleme = ipt1.lisous(isous)
  405. segact meleme
  406. KEL =itypel
  407. nbnn=num(/1)
  408. C
  409. kd = kdegre(kel)
  410. nso = nbsom(kel)
  411. iad=nspos(kel)-1
  412. segini wrk4
  413. segini inomil
  414. C element quaf ? si oui kell=numero lineaire correspondant
  415. C si non kell=0
  416. CALL NQF2NL(kel,kell)
  417. IF (IERR.NE.0) RETURN
  418. IF (KELL.NE.0) THEN
  419. nbnnl=NBNNE(KELL)
  420. segini wrkl
  421. do i=1,nbnnl
  422. inomil(**)=i
  423. enddo
  424. ELSE
  425. if(kd.eq.3) then
  426. C element quadratiques
  427. do 762 i=1,nbnn
  428. do 763 j=1,nso
  429. iso = ibsom(iad+j)
  430. if(i.eq.iso) goto 762
  431. 763 continue
  432. inomil(**)= i
  433. 762 continue
  434. elseif(kd.eq.2) then
  435. C elements lineaires
  436. do i=1,nbnn
  437. inomil(**)=i
  438. enddo
  439. endif
  440. ENDIF
  441. itest=1
  442. do 101 iel=1,num(/2)
  443. call doxe(xcoor,idim,nbnn,num,iel,xel)
  444. * cas quaf
  445. if (kell.ne.0) then
  446. nsom=nbsom(kel)
  447. if(nsom.ne.nbnnl) then
  448. call erreur(5)
  449. return
  450. endif
  451. idx=nspos(kel)-1
  452. do innl=1,nbnnl
  453. do iid=1,idim
  454. xell(iid,innl)=xel(iid,ibsom(idx+innl))
  455. enddo
  456. enddo
  457. endif
  458.  
  459. iref = num(1,iel)*idim1 - idim
  460. xmi=xcoor(iref)
  461. ymi=xcoor(iref+1)
  462. zmi=xcoor(iref+2)
  463. xma=xmi
  464. yma=ymi
  465. zma=zmi
  466. do 102 ipo=2,num(/1)
  467. iref=num(ipo,iel)*idim1-idim
  468. if(xcoor(iref).lt.xmi) xmi = xcoor(iref)
  469. if(xcoor(iref).gt.xma) xma = xcoor(iref)
  470. if(xcoor(iref+1).lt.ymi) ymi = xcoor(iref+1)
  471. if(xcoor(iref+1).gt.yma) yma = xcoor(iref+1)
  472. if(xcoor(iref+2).lt.zmi) zmi = xcoor(iref+2)
  473. if(xcoor(iref+2).gt.zma) zma = xcoor(iref+2)
  474. 102 continue
  475. * if(iel.lt.2.and.isous.eq.1) then
  476. * WRITE(IOIMP,*) ' element ',iel
  477. * WRITE(IOIMP,*) ( num(iuo,iel),iuo=1,num(/1))
  478. * WRITE(IOIMP,*) ' xmi ymi zmi xma yma zma'
  479. * WRITE(IOIMP,*) ' boite reelle'
  480. * WRITE(IOIMP,*) xmi,ymi,zmi,xma,yma,zma
  481. * endif
  482. xdec1=(xma-xmi)/3.D0
  483. ydec1=(yma-ymi)/3.D0
  484. zdec1=(zma-zmi)/3.D0
  485. xmi=xmi-xdec1
  486. xma=xma+xdec1
  487. ymi=ymi-ydec1
  488. yma=yma+ydec1
  489. zmi=zmi-zdec1
  490. zma=zma+zdec1
  491. if(xmi.lt.xmin) xmi=xmin
  492. if(xmi.gt.xmax) xmi=xmax
  493. if(xma.lt.xmin) xma=xmin
  494. if(xma.gt.xmax) xma=xmax
  495. if(ymi.lt.ymin) ymi=ymin
  496. if(ymi.gt.ymax) ymi=ymax
  497. if(yma.lt.ymin) yma=ymin
  498. if(yma.gt.ymax) yma=ymax
  499. if(zmi.lt.zmin) zmi=zmin
  500. if(zmi.gt.zmax) zmi=zmax
  501. if(zma.lt.zmin) zma=zmin
  502. if(zma.gt.zmax) zma=zmax
  503. * if(iel.lt.2.and.isous.eq.1) then
  504. *
  505. * WRITE(IOIMP,*) ' boite apres modif'
  506. * WRITE(IOIMP,*) xmi,ymi,zmi,xma,yma,zma
  507. * endif
  508. * WRITE(IOIMP,*) ' xmi xma ymi yma zmi zma',xmi,xma,ymi,yma,zmi,zma
  509. * if(idim.eq.2) then
  510. * zmi=0.d0
  511. * zma=0.d0
  512. * endif
  513. C
  514. C
  515. C recherche des points interne à l'element
  516. C pour cela on balaye les super zones
  517. C
  518. if(idim.eq.2) then
  519.  
  520. IXDEP=INT(real((Xmi-XMIN)/XDEC))+1+
  521. & INT(real((ymi-YMIN)/YDEC))*NDEC
  522.  
  523.  
  524. ixpas=INT(real((Xma-XMIN)/XDEC))+1+
  525. & INT(real((ymi-YMIN)/YDEC))*NDEC - ixdep
  526.  
  527. iny=INT(real((yma-YMIn)/YDEC))-INT(real((ymi-YMIn)
  528. $ /YDEC))+1
  529.  
  530. do 120 ily=1,iny
  531. ile=ixdep +( ily -1 )* ndec
  532. ilf=ile + ixpas
  533. izv= ile/nf+1
  534. * if (izv.gt.ite) WRITE(IOIMP,*) ' prob ',izv,ixpas,ixdep,
  535. * > ndec,mdec,xdec,ydec,ymi,ymin
  536. jd=np1(izv)+1
  537. do 122 j=jd,numnp3
  538. izd=np3(j)
  539. if(izd.lt.ile) go to 122
  540. if(izd.gt.ilf) go to 122
  541. ipo= np2(j)
  542. C recherche si le point est interne
  543. if(idejvu(ipo).eq.1) go to 122
  544. ip=ipt2.num(1,ipo)
  545. iref=ip*idim1-idim
  546. xpu(1)=xcoor(iref)
  547. xpu(2)=xcoor(iref+1)
  548. if(xpu(1).lt.xmi.or.xpu(1).gt.xma) go to 122
  549. if(xpu(2).lt.ymi.or.xpu(2).gt.yma) go to 122
  550. * cas quaf
  551. if (kell.ne.0) then
  552. call qsijs(xell,kell,nbnnl,idim,xpu,
  553. $ SHPL,qsi,iret)
  554. call qsijs2(xel,kel,nbnn,idim,xpu,
  555. $ SHPP,qsi,iret)
  556. else
  557. call qsijs(xel,kel,nbnn,idim,xpu,SHPP,qsi,iret)
  558. endif
  559. IF (IRET.NE.0) then
  560. irate=irate+1
  561. GOTO 122
  562. endif
  563. IF (KELL.NE.0) THEN
  564. CALL SHAPE(qsi(1),qsi(2),qsi(3),kell,shpl
  565. $ ,iret2)
  566. if (iret2.eq.0) then
  567. MOTERR(1:4)=NOMS(KELL)
  568. CALL ERREUR(68)
  569. RETURN
  570. endif
  571. ENDIF
  572.  
  573. * WRITE(IOIMP,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn)
  574.  
  575. do i=1,inomil(/1)
  576. ilp= inomil(i)
  577. if (kell.ne.0) then
  578. if( shpl(1,ilp).lt.0.D0) goto 110
  579. else
  580. if( shpp(1,ilp).lt.0.D0) goto 110
  581. endif
  582. enddo
  583. C
  584. goto 100
  585. 110 xmap=1.D10
  586. do i=1,inomil(/1)
  587. ilp=inomil(i)
  588. if (kell.ne.0) then
  589. xmap = min ( xmap,shpl(1,ilp))
  590. else
  591. xmap = min ( xmap,shpp(1,ilp))
  592. endif
  593. enddo
  594. if( xmap . gt . ckrit(ipo) ) then
  595. do 180 jl =1,NUM(/1)
  596. if(IP.NE.NUM(jl,iel)) go to 180
  597. if(IFL.EQ.0) go to 190
  598. idejvu(ipo)=1
  599. itr=itr+1
  600. go to 181
  601. 180 continue
  602. 181 continue
  603. ckrit(ipo)=xmap
  604. ieint(1,ipo)=isous
  605. ieint(2,ipo)=iel
  606. QQQ(1,ipo) = qsi(1)
  607. QQQ(2,ipo) = qsi(2)
  608. endif
  609. 190 continue
  610. go to 122
  611. 100 CONTINUE
  612.  
  613. C on a trouve l element contenant le point ipo attention
  614. C si le point appartient aux deux maillages on ne le stocke pas
  615. C pour accro mais on le stocke pour proi
  616. do 18 jl =1,NUM(/1)
  617. if((IP.EQ.NUM(jl,iel)).AND.(IFL.EQ.0)) GOTO
  618. $ 19
  619. 18 continue
  620. if(idejvu(ipo).EQ.0) then
  621. idejvu(ipo) = 1
  622. itr = itr + 1
  623. IEINT(1,ipo) = isous
  624. IEINT(2,ipo) = iel
  625. QQQ(1,ipo) = qsi(1)
  626. QQQ(2,ipo) = qsi(2)
  627. * WRITE(IOIMP,*) 'trouve point' ,IP, 'Element' ,J
  628. * WRITE(IOIMP,2375)ip,j,(xpu(i1),i1=1,idim),(qsi(i2),i2=1,idim)
  629. endif
  630. 19 continue
  631. * 2375 format(2i4,6e12.5)
  632. 122 continue
  633. 120 continue
  634. endif
  635. if(idim.eq.3) then
  636.  
  637.  
  638.  
  639. IXDEP=INT(real((Xmi-XMIN)/XDEC))+1+
  640. & INT(real((ymi-YMIN)/YDEC))*NDEC +
  641. & INT(real((zmi-zmin)/ZDEC))*NMDEC
  642.  
  643.  
  644. ixpas=INT(real((Xma-XMIN)/XDEC))+1+
  645. & INT(real((ymi-YMIN)/YDEC))*NDEC +
  646. & INT(real((zmi-zmin)/ZDEC))*NMDEC - ixdep
  647.  
  648.  
  649. iny=INT(real((yma-YMIn)/YDEC))-INT(real((ymi-YMIn)
  650. $ /YDEC))+1
  651. inz= INT(real((zma-zmin)/ZDEC))-INT(real((zmi-zmin)
  652. $ /ZDEC))+1
  653.  
  654. * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) ' inz iny',inz,iny
  655. do 240 itz=1,inz
  656. idexx= ixdep + (itz-1)*nmdec
  657. do 220 ily=1,iny
  658. ile=idexx +( ily -1 )* ndec
  659. ilf=ile + ixpas
  660. * if(iel.lt.2.and.isous.eq.1)WRITE(IOIMP,*)' ile,ilf' , ile,ilf
  661. izv= ile/nf + 1
  662. jd=np1(izv)+1
  663. * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) ' jd,jf ',jd,jf
  664. do 222 j=jd,numnp3
  665. izd=np3(j)
  666. * WRITE(IOIMP,*) ' izd ', izd
  667. if(izd.lt.ile) go to 222
  668. if(izd.gt.ilf) go to 223
  669. ipo= np2(j)
  670. C recherche si le point est interne
  671. if(idejvu(ipo).eq.1) go to 222
  672. ip=ipt2.num(1,ipo)
  673. * if(iel.lt.2.and.isous.eq.1)WRITE(IOIMP,*) ' ipo ip ',ipo,ip
  674. iref=ip*idim1-idim
  675. xpu(1)=xcoor(iref)
  676. xpu(2)=xcoor(iref+1)
  677. xpu(3)=xcoor(iref+2)
  678.  
  679. * if(iel.lt.2.and.isous.eq.1) WRITE(IOIMP,*) 'cooor '
  680. * $ ,xpu(1),xpu(2),xpu(3)
  681. if(xpu(1).lt.xmi.or.xpu(1).gt.xma) go to 222
  682. if(xpu(2).lt.ymi.or.xpu(2).gt.yma) go to 222
  683. if(xpu(3).lt.zmi.or.xpu(3).gt.zma) go to 222
  684. call qsijs(xel,kel,nbnn,idim,xpu,SHPP,qsi
  685. $ ,iret)
  686. IF (IRET.NE.0) THEN
  687. irate=irate+1
  688. if(irate.eq.itest) then
  689. * WRITE(IOIMP,*) 'solution pas trouvee', irate
  690. itest=itest *2
  691. endif
  692. GOTO 222
  693. endif
  694. IF (KELL.NE.0) THEN
  695. CALL SHAPE(qsi(1),qsi(2),qsi(3),kell,shpl
  696. $ ,iret2)
  697. ENDIF
  698. C
  699. * WRITE(IOIMP,fmt='(8E12.5)')(shpp(1,ih),ih=1,nbnn)
  700. do i=1,inomil(/1)
  701. ilp= inomil(i)
  702. if (kell.ne.0) then
  703. if( shpl(1,ilp).lt.0.D0) goto 210
  704. else
  705. if( shpp(1,ilp).lt.0.D0) goto 210
  706. endif
  707. enddo
  708. C
  709. goto 300
  710. 210 xmap=1.D10
  711. do i=1,inomil(/1)
  712. ilp=inomil(i)
  713. if (kell.ne.0) then
  714. xmap = min ( xmap,shpl(1,ilp))
  715. else
  716. xmap = min ( xmap,shpp(1,ilp))
  717. endif
  718. enddo
  719. if( xmap . gt . ckrit(ipo) ) then
  720. do 280 jl =1,NUM(/1)
  721. if(IP.NE.NUM(jl,iel)) go to 280
  722. if(IFL.EQ.0) go to 290
  723. idejvu(ipo)=1
  724. itr=itr+1
  725. * if(iel.lt.3.and.isous.eq.1) WRITE(IOIMP,*)'egalite de point'
  726. C
  727. C
  728. C
  729. go to 281
  730. 280 continue
  731. 281 continue
  732. ckrit(ipo)=xmap
  733. ieint(1,ipo)=isous
  734. ieint(2,ipo)=iel
  735. QQQ(1,ipo) = qsi(1)
  736. QQQ(2,ipo) = qsi(2)
  737. QQQ(3,ipo) = qsi(3)
  738. endif
  739. 290 continue
  740. go to 222
  741. 300 CONTINUE
  742.  
  743. C on a trouve l element contenant le point ipo attention
  744. C si le point appartient aux deux maillages on ne le stocke pas
  745. C pour accro mais on le stocke pour proi
  746. do 28 jl =1,NUM(/1)
  747. if((IP.EQ.NUM(jl,iel)).AND.(IFL.EQ.0))
  748. $ GOTO 29
  749. 28 continue
  750. if(idejvu(ipo).EQ.0) then
  751.  
  752. idejvu(ipo) = 1
  753. * if(iel.lt.3.and.isous.eq.1) WRITE(IOIMP,*)'point interne'
  754. itr = itr + 1
  755. IEINT(1,ipo) = isous
  756. IEINT(2,ipo) = iel
  757. QQQ(1,ipo) = qsi(1)
  758. QQQ(2,ipo) = qsi(2)
  759. QQQ(3,ipo) = qsi(3)
  760. * if(iel.lt.3.and.isous.eq.1) then
  761. * WRITE(IOIMP,*) 'trouve point' ,IP, 'Element' ,iel
  762. * WRITE(IOIMP,2375)ip,iel,(xpu(i1),i1=1,idim),(qsi(i2),i2=1,idim
  763. C )
  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. CALL QSIJS(xel,kel,nbnn,IDIM,xpu,shpp,qsi,IRET)
  989. C- WRITE(IOIMP,FMT='(8E12.5)') (shpp(1,i),i=1,nbnn)
  990. IF (IRET.NE.0) then
  991. irate=irate+1
  992. GOTO 630
  993. ENDIF
  994.  
  995. do i=1,inomil(/1)
  996. ilp= inomil(i)
  997. if( shpp(1,ilp).lt.0.D0) goto 1110
  998. enddo
  999. C
  1000. goto 1100
  1001. 1110 xmap=1.D10
  1002. do i=1,inomil(/1)
  1003. ilp=inomil(i)
  1004. xmap = min ( xmap,shpp(1,ilp))
  1005. enddo
  1006. if( xmap . gt . ckrit(ipo) ) then
  1007. do 1180 jl =1,NUM(/1)
  1008. if(IP.NE.NUM(jl,iel)) go to 1180
  1009. if(IFL.EQ.0) go to 1190
  1010. idejvu(ipo)=1
  1011. itr=itr+1
  1012. go to 1181
  1013. 1180 continue
  1014. 1181 continue
  1015. ckrit(ipo)=xmap
  1016. ieint(1,ipo)=isous
  1017. ieint(2,ipo)=iel
  1018. QQQ(1,ipo) = qsi(1)
  1019. endif
  1020. 1190 continue
  1021. go to 630
  1022. 1100 CONTINUE
  1023.  
  1024.  
  1025.  
  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. RETURN
  1060. END
  1061.  
  1062.  
  1063.  
  1064.  
  1065.  
  1066.  
  1067.  
  1068.  
  1069.  
  1070.  
  1071.  
  1072.  
  1073.  
  1074.  
  1075.  
  1076.  
  1077.  
  1078.  
  1079.  

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