Télécharger demait.eso

Retour à la liste

Numérotation des lignes :

  1. C DEMAIT SOURCE JC220346 16/11/29 21:15:09 9221
  2. C|-------------------------------------------------------------------|
  3. C| |
  4. C| PROGRAMME PRINCIPAL |
  5. C| MAIN DE DEMETER |
  6. C| |
  7. C|-------------------------------------------------------------------|
  8. C
  9. SUBROUTINE DEMAIT(IDCP,NPTBAS)
  10. C
  11. C
  12. IMPLICIT INTEGER(I-N)
  13. IMPLICIT REAL*8(a-h,o-z)
  14. -INC CCOPTIO
  15. -INC TDEMAIT
  16. SEGMENT IDCP(NPTCOM)
  17. REAL*8 XDEM
  18. LOGICAL REPONS,FACET,INTER,PINT,FERME,VAL,IN,DROIT,
  19. real*8 xval,xo(3),xa(3),xb(3),xc(3),epsi,epsj,xmesu,ymesu
  20. character*4 mcle(8)
  21. data mcle /'TCRI','CFAC','CDIS','TETR','EXPC','FINA','DIAC',
  22. > 'EXPF'/
  23. ymesu=0
  24. NPTORI=NPTMAX
  25. NFCORI=NFCMAX
  26. * flags type d'operation
  27. ipass=1
  28. tcrit=3.
  29. cfacei=6.0
  30. cfacet=cfacei
  31. cdist=0.125
  32. tetrl=2.75
  33. expcom=1.00
  34. diacri=0.93
  35. diacrd=diacri
  36. diacre=diacri
  37. expfac=sqrt(3.)/2.
  38. faccri=16
  39. volcri=0.01
  40. 10 continue
  41. call lirmot(mcle,8,imot,0)
  42. if (imot.eq.1) then
  43. call lirree(xval,1,iret)
  44. tcrit=xval
  45. endif
  46. if (imot.eq.2) then
  47. call lirree(xval,1,iret)
  48. cfacei=xval
  49. endif
  50. if (imot.eq.3) then
  51. call lirree(xval,1,iret)
  52. cdist=xval
  53. endif
  54. if (imot.eq.4) then
  55. call lirree(xval,1,iret)
  56. tetrl=xval
  57. endif
  58. if (imot.eq.5) then
  59. call lirree(xval,1,iret)
  60. expcom=xval
  61. endif
  62. if (imot.eq.7) then
  63. call lirree(xval,1,iret)
  64. diacre=xval
  65. endif
  66. if (imot.eq.8) then
  67. call lirree(xval,1,iret)
  68. expfac=xval
  69. endif
  70. if (imot.ne.0) goto 10
  71. C
  72. C INITIALISATION DU TABLEAU DES VOLUMES
  73. C -------------------------------------
  74. NFTOT=IFUT(/1)
  75. NVTOT=IVOL(/2)
  76. DO 130 J=1,NVTOT
  77. DO 120 I=1,9
  78. IVOL(I,J)=0
  79. 120 CONTINUE
  80. 130 CONTINUE
  81. C
  82. C
  83. C CONSTRUCTION DU TABLEAU NPF ( POINTS-FACETTES )
  84. C -----------------------------------------------
  85. DO 150 J=1,40
  86. DO 150 I=1,NPTMAX
  87. NPF(J,I)=0
  88. 150 CONTINUE
  89. DO 141 J=1,NFCMAX
  90. DO 140 K=1,4
  91. IP=NFC(K,J)
  92. IF (IP.EQ.0) GOTO 140
  93. L=-NPF(40,IP)+1
  94. IF (L.LE.0) CALL ERREUR(126)
  95. IF (IERR.NE.0) RETURN
  96. NPF(40,IP)=-L
  97. NPF(L,IP)=J
  98. 140 CONTINUE
  99. 141 CONTINUE
  100. DO 145 I=1,NPTMAX
  101. NPF(40,I)=MAX(0,NPF(40,I))
  102. 145 CONTINUE
  103. C
  104. C RECHERCHE DE LA TAILLE MOYENNE DE MAILLE ASSOCIEE A
  105. C CHAQUE POINT ( 4EME COMPOSANTE DU POINT )
  106. DO 190 I=1,NPTMAX
  107. DD=0.
  108. KK=0
  109. DO 170 J=1,40
  110. IF (NPF(J,I).EQ.0) GOTO 180
  111. if=npf(j,i)
  112. nc=4
  113. if (nfc(4,if).eq.0) nc=3
  114. jp=nfc(nc,if)
  115. do 175 ic=1,nc
  116. ip=nfc(ic,if)
  117. XX=(XYZ(1,IP)-XYZ(1,JP))**2
  118. YY=(XYZ(2,IP)-XYZ(2,JP))**2
  119. ZZ=(XYZ(3,IP)-XYZ(3,JP))**2
  120. DD=DD+SQRT(XX+YY+ZZ)
  121. KK=KK+1
  122. jp=ip
  123. 175 continue
  124. * IP=ISUCC(NPF(J,I),I)
  125. * XX=(XYZ(1,I)-XYZ(1,IP))**2
  126. * YY=(XYZ(2,I)-XYZ(2,IP))**2
  127. * ZZ=(XYZ(3,I)-XYZ(3,IP))**2
  128. * DD=DD+SQRT(XX+YY+ZZ)
  129. * KK=KK+1
  130. 170 CONTINUE
  131. 180 XYZ(4,I)=DD/KK
  132. 190 CONTINUE
  133. * regularisation locale
  134. * DO 182 I=1,NPTMAX
  135. * DD= XYZ(4,I)
  136. * KK=1
  137. * DO 184 J=1,40
  138. * IF (NPF(J,I).EQ.0) GOTO 186
  139. * IP=ISUCC(NPF(J,I),I)
  140. * DD=DD+XYZ(4,IP)
  141. * KK=KK+1
  142. *84 CONTINUE
  143. *86 XYZ(4,I)=DD/KK
  144. *82 CONTINUE
  145. * taille moyenne generale
  146. xmoy=0
  147. DO 181 I=1,NPTMAX
  148. XMOY=XMOY+LOG(XYZ(4,I))
  149. 181 CONTINUE
  150. XMOYG=EXP(XMOY/NPTMAX)
  151. IF (IVERB.EQ.1) WRITE (6,*) ' TAILLE MOYENNE VISEE ',XMOYG
  152. C
  153. C LE MAILLAGE DE LA SURFACE EST-IL FERME ?
  154. C ----------------------------------------
  155. REPONS=FERME(KKK)
  156. IF (REPONS) GOTO 210
  157. CALL ERREUR(127)
  158. RETURN
  159. 210 continue
  160. xmesu=0
  161. do 100 if=1,nfcmax
  162. xmesu=xmesu+vol(1,nfc(1,if),nfc(2,if),nfc(3,if))
  163. if (nfc(4,if).ne.0) then
  164. xmesu=xmesu+vol(1,nfc(1,if),nfc(3,if),nfc(4,if))
  165. ipass=2
  166. endif
  167. 100 continue
  168. xmesu=xmesu/(-6.)
  169. IF (IVERB.EQ.1) WRITE (6,*) ' volume de la piece ',xmesu
  170. C
  171. C NFACET : NOMBRE DE FACETTES DU MAILLAGE DE SURFACE
  172. C --------------------------------------------------
  173. NFACET=NFCMAX
  174. NPTCOM=NPTMAX
  175. NVOL=0
  176. NFCPRE=0
  177. NFCTRT=0
  178. NARET=0
  179. NPTOT=XYZ(/2)
  180. NTTRAV=NFCMAX
  181. SEGINI TRAV
  182. YMESU=0
  183. NVOLY=0
  184. NPTDEB=NPTMAX
  185. NPTDIS=1
  186. ICRTS=0
  187. 220 CONTINUE
  188. do 222 jvol=nvoly+1,nvol
  189. if (ivol(9,jvol).eq.25) then
  190. ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol),
  191. > ivol(3,jvol),ivol(4,jvol))/6.
  192. endif
  193. if (ivol(9,jvol).eq.35) then
  194. ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol),
  195. > ivol(3,jvol),ivol(5,jvol))/6.
  196. > +vol(ivol(1,jvol),ivol(3,jvol),
  197. > ivol(4,jvol),ivol(5,jvol))/6.
  198. endif
  199. if (ivol(9,jvol).eq.30) then
  200. ymesu=ymesu+vol(ivol(1,jvol),ivol(2,jvol),
  201. > ivol(3,jvol),ivol(4,jvol))/6.
  202. > +vol(ivol(2,jvol),ivol(3,jvol),
  203. > ivol(4,jvol),ivol(5,jvol))/6.
  204. > +vol(ivol(3,jvol),ivol(5,jvol),
  205. > ivol(6,jvol),ivol(4,jvol))/6.
  206. endif
  207. if (ivol(9,jvol).eq.20) then
  208. ymesu=ymesu-vol(ivol(1,jvol),ivol(3,jvol),
  209. > ivol(6,jvol),ivol(8,jvol))/6.
  210. > -vol(ivol(5,jvol),ivol(6,jvol),
  211. > ivol(8,jvol),ivol(1,jvol))/6.
  212. > -vol(ivol(2,jvol),ivol(6,jvol),
  213. > ivol(1,jvol),ivol(3,jvol))/6.
  214. ymesu=ymesu-vol(ivol(7,jvol),ivol(8,jvol),
  215. > ivol(6,jvol),ivol(3,jvol))/6.
  216. > -vol(ivol(4,jvol),ivol(1,jvol),
  217. > ivol(8,jvol),ivol(3,jvol))/6.
  218. endif
  219. 222 continue
  220. nvoly=nvol
  221. if (ymesu.gt.xmesu*1.01) goto 340
  222. if (ierr.ne.0) goto 340
  223. * AJUSTEMENT EVENTUEL DES DIMENSIONS DES TABLEAUX
  224. IF (NFCMAX+250.GE.NFTOT) THEN
  225. NFTOT=NFTOT+300
  226. SEGADJ NFC,NFV,IFUT,IFAT
  227. ENDIF
  228. IF (NVOL+10.GE.NVTOT) THEN
  229. NVTOT=NVTOT+50
  230. SEGADJ IVOL
  231. ENDIF
  232. IF (NPTMAX+50.GE.NPTOT) THEN
  233. NPTOT=NPTOT+100
  234. SEGADJ NPF,XYZ
  235. ENDIF
  236. * nouvelle methode de calcul de la taille locale
  237. DO 221 I=NPTDEB+1,NPTMAX
  238. call vcrit(i)
  239. 221 CONTINUE
  240. NPTDEB=NPTMAX
  241. IGAGNE=0
  242. IF (NFACET.EQ.0) GOTO 370
  243. NFCMA=NFCMAX
  244. C
  245. C
  246. C RECHERCHE DES DIEDRES
  247. C ---------------------
  248. * FAIRE ICI LE MENAGE DANS NARET(de temps en temps)
  249. IF (ICRTS.GE.100) THEN
  250. NPTDIS=NPTMAX
  251. ICRTS=0
  252. IVA=0
  253. DO 285 I=1,NARET
  254. if (IIGARD(I).LE.0) goto 285
  255. IF1=IF1GAR(I)
  256. IF (IFAT(IF1).EQ.0) GOTO 285
  257. IF2=IF2GAR(I)
  258. IF (IFAT(IF2).EQ.0) GOTO 285
  259. IVA=IVA+1
  260. II=IIGARD(I)
  261. NPTDIS=MIN(NPTDIS,II)
  262. IIGARD(IVA)=II
  263. JJ=JJGARD(I)
  264. NPTDIS=MIN(NPTDIS,JJ)
  265. JJGARD(IVA)=JJ
  266. IF1=IF1GAR(I)
  267. IF1GAR(IVA)=IF1
  268. IF2=IF2GAR(I)
  269. IF2GAR(IVA)=IF2
  270. ANGAR(IVA)=ANGAR(I)
  271. IF (IVA.NE.1) THEN
  272. ANGMA(IVA)=MAX(ANGAR(IVA),ANGMA(IVA-1))
  273. ELSE
  274. ANGMA(1)=ANGAR(1)
  275. ENDIF
  276. 285 CONTINUE
  277. * write (6,*) ' demait retassement effectue ',naret,iva
  278. NARET=IVA
  279. ENDIF
  280. DO 290 IF1=NFCTRT+1,NFCMAX
  281. IF (IFAT(IF1).EQ.0) GOTO 290
  282. NBD=4
  283. IF (NFC(4,IF1).EQ.0) NBD=3
  284. DO 292 IC=1,NBD
  285. IC1=IC-1
  286. IF (IC1.EQ.0) IC1=NBD
  287. II=NFC(IC1,IF1)
  288. JJ=NFC(IC,IF1)
  289. DO 294 I=1,40
  290. IF2=NPF(I,II)
  291. IF (IF2.EQ.0) GOTO 292
  292. IF (IF2.GE.IF1) GOTO 294
  293. IF (IPRED(IF2,II).NE.JJ) GOTO 294
  294. C
  295. C COMMENT SONT LES ELEMENTS DU DIEDRE ?
  296. C -------------------------------------
  297. C
  298. ANGLL=TETA(IF1,IF2,II,JJ)
  299. C
  300. * pour gagner du temps
  301. if (angll.le.-1d4) goto 294
  302. * if (if1.le.nfcori.or.if2.le.nfcori) angll=angll+1d6
  303. NARET=NARET+1
  304. IF (NARET.GT.NTTRAV) THEN
  305. NTTRAV=NARET+20
  306. SEGADJ TRAV
  307. ENDIF
  308. IIGARD(NARET)=II
  309. JJGARD(NARET)=JJ
  310. IF1GAR(NARET)=IF1
  311. IF2GAR(NARET)=IF2
  312. ANGAR(NARET)=ANGLL
  313. IF (NARET.NE.1) THEN
  314. ANGMA(NARET)=MAX(ANGAR(NARET),ANGMA(NARET-1))
  315. ELSE
  316. ANGMA(1)=ANGAR(1)
  317. ENDIF
  318. 294 CONTINUE
  319. 292 CONTINUE
  320. 290 CONTINUE
  321. NFCTRT=NFCMAX
  322. *
  323. * ON COMMENCE PAR L'ANGLE LE PLUS FERME
  324. *
  325. 315 CONTINUE
  326. ANGMAX=-1.E30
  327. IOK=0
  328. DO 310 I=NARET,1,-1
  329. IF(ANGMAX.GE.ANGMA(I)) GOTO 311
  330. IF(IIGARD(I).LE.0) GOTO 310
  331. IF(ANGMAX.GE.ANGAR(I)) GOTO 310
  332. IOK=I
  333. * if (angar(i).gt.0d0) ANGMAX=ANGAR(I)*0.99999D0
  334. * if (angar(i).lt.0d0) ANGMAX=ANGAR(I)*1.00001D0
  335. angmax=angar(i)-1d-6
  336.  
  337. 310 CONTINUE
  338. 311 CONTINUE
  339. IF (IOK.EQ.0) GOTO 320
  340. II=IIGARD(IOK)
  341. JJ=JJGARD(IOK)
  342. IF1=IF1GAR(IOK)
  343. IF2=IF2GAR(IOK)
  344. IIGARD(IOK)=-II
  345. ICRTS=ICRTS+1
  346. IF (IFAT(IF1).EQ.0) GOTO 315
  347. IF (IFAT(IF2).EQ.0) GOTO 315
  348. IGAGNE=0
  349. * write (6,*) 'traitement diedre ',ii,jj,if1,if2,angmax
  350. idiac=5
  351. 313 continue
  352. * on essaie d'abord de faire les hexaedres
  353. if (ipass.eq.2) then
  354. IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).NE.0)
  355. # CALL hexa(II,JJ,IF1,IF2,IGAGNE)
  356. GOTO 312
  357. else
  358. IF (NFC(4,IF1).EQ.0.AND.NFC(4,IF2).EQ.0)
  359. # CALL CONS33(II,JJ,IF1,IF2,IGAGNE,0)
  360. IF (IGAGNE.EQ.1) GOTO 312
  361. IF (NFC(4,IF1).EQ.0.AND.NFC(4,IF2).NE.0)
  362. # CALL CONS34(II,JJ,IF1,IF2,IGAGNE)
  363. IF (IGAGNE.EQ.1) GOTO 312
  364. IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).EQ.0)
  365. # CALL CONS34(JJ,II,IF2,IF1,IGAGNE)
  366. IF (IGAGNE.EQ.1) GOTO 312
  367. IF (NFC(4,IF1).NE.0.AND.NFC(4,IF2).NE.0)
  368. # CALL CONS44(II,JJ,IF1,IF2,IGAGNE)
  369. IF (IGAGNE.EQ.1) GOTO 312
  370. endif
  371. * write (6,*) ' demait relachement de diacrd'
  372. diacrd=1-0.3*(1-diacrd)
  373. diacri=diacrd
  374. cfacet=cfacet*1.5
  375. faccri=faccri*2
  376. cdist=0.085
  377. tetrl=tetrl*1.3
  378. idiac=idiac-1
  379. if (idiac.ne.0) goto 313
  380. diacrd=diacre
  381. diacri=diacrd
  382. cfacet=cfacei
  383. faccri=16
  384. cdist=0.125
  385. tetrl=2.75
  386. IF (IVERB.EQ.1) write (6,*) ' demait echec traitement diedre',
  387. > ii,jj,if1,if2,angmax
  388. goto 315
  389. 312 CONTINUE
  390. diacrd=diacre
  391. diacri=diacrd
  392. cfacet=cfacei
  393. faccri=16
  394. cdist=0.125
  395. tetrl=2.75
  396. if (ipass.eq.-2) then
  397. ipass=-1
  398. * cfacet=6
  399. * cdist=0.125
  400. * tetrl=2.75
  401. * faccri=16
  402. * volcri=0.2
  403. * diacri=0.90
  404. * diacrd=0.92
  405. endif
  406. GOTO 220
  407. 320 CONTINUE
  408. IF (ipass.eq.2) then
  409. ipass=1
  410. IF (IVERB.EQ.1) write (6,*) ' fin generation de cube '
  411. * on remet les compteurs a zero
  412. nfctrt=0
  413. naret =0
  414. goto 220
  415. endif
  416. IF (ipass.eq.1) then
  417. ipass=0
  418. IF (IVERB.EQ.1) write (6,*) ' strategie finale 1'
  419. cfacet=16
  420. cdist=0.085
  421. tetrl=9
  422. diacri=0.95
  423. diacrd=0.95
  424. faccri=64
  425. volcri=0.01
  426. * on remet les compteurs a zero
  427. nfctrt=0
  428. naret =0
  429. goto 220
  430. endif
  431. IF (ipass.eq.0) then
  432. ipass=-1
  433. IF (IVERB.EQ.1) write (6,*) ' strategie finale 2'
  434. cfacet=100
  435. cdist=0.050
  436. tetrl=9
  437. diacri=0.99
  438. diacrd=0.99
  439. faccri=64
  440. volcri=0.01
  441. * on remet les compteurs a zero
  442. nfctrt=0
  443. naret =0
  444. goto 220
  445. endif
  446. IF (ipass.eq.-1) then
  447. ipass=-2
  448. IF (IVERB.EQ.1) write (6,*) ' strategie finale 3'
  449. cfacet=1000
  450. cdist=0.005
  451. tetrl=9
  452. diacri=0.999
  453. diacrd=0.999
  454. faccri=64
  455. volcri=0.01
  456. * on remet les compteurs a zero
  457. nfctrt=0
  458. naret =0
  459. goto 220
  460. endif
  461. 340 continue
  462. DO 444 I=1,NFCMAX
  463. IF (IFAT(I).EQ.0) GOTO 444
  464. IF (IVERB.EQ.1) WRITE (6,*) ' FACETTE RESTANTE IF ',I,
  465. * NFC(1,I),NFC(2,I),NFC(3,I),NFC(4,I)
  466. 444 CONTINUE
  467. 4001 CONTINUE
  468. IF (NFACET.NE.0) CALL ERREUR(27)
  469. 370 CONTINUE
  470. IF (IVERB.EQ.1) WRITE (6,*) ' volume du maillage ',ymesu
  471. IF (IVERB.EQ.1) write (6,*) ' nb facette ',nfacet
  472. SEGSUP TRAV
  473. CALL OPTVOL
  474. RETURN
  475. C FIN DU PROGRAMME PRINCIPAL
  476. END
  477.  
  478.  
  479.  
  480.  
  481.  
  482.  
  483.  
  484.  
  485.  

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