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

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