Télécharger envvo2.eso

Retour à la liste

Numérotation des lignes :

  1. C ENVVO2 SOURCE CB215821 19/08/20 21:17:04 10287
  2. SUBROUTINE ENVVO2(icle)
  3. *
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. *
  7. * 20160713 : SG nouvelle programmation de envvol, simplifiée, en vue
  8. * de traiter les faces TRI7/QUA9 nouvellement rajoutées
  9. * dans le bdata. La simplification sera egalement utile pour
  10. * les subroutines inspirees de envvol : faced2, faced, envel1,
  11. * dicho3
  12. * icle=0 : renvoie l'enveloppe (operateur ENVE)
  13. * icle=1 : renvoie les face (operateur CHAN FACE)
  14.  
  15.  
  16. -INC PPARAM
  17. -INC CCOPTIO
  18. -INC SMELEME
  19. -INC CCGEOME
  20. -INC SMCOORD
  21. * Type de faces prises en compte: T3, Q4, T6, Q8, T7, Q9
  22. * Numero dans KDFAC 1 2 3 4 7 8
  23. * Pour ne pas se prendre la tête, on numerote pareil que KDFAC
  24. * Pour les 5 (non utilisé), 6 (polygone) et >8, ca restera à 0
  25. * NTYFAC=20, codé en dur dans CCGEOME pour KDFAC
  26. PARAMETER (NTYFAC=20)
  27. * Nb de faces de chaque type, sert également de compteur
  28. SEGMENT NBFAC(NTYFAC)
  29. * Tableau d'index de début fin dans les tableaux ???(NFAC)
  30. SEGMENT IDXFAC(NTYFAC+1)
  31. * Pointeurs sur des segments MELEME par type de face
  32. SEGMENT IPTFAC(NTYFAC)
  33. * Un segment pointant sur les IFACI
  34. SEGMENT IPOFAC(NTYFAC)
  35. * Segment IFACI contenant les noeuds et la couleurs des faces d'un
  36. * type donné
  37. SEGMENT IFACI(NNODE+1,NFACI)
  38. *
  39. SEGMENT IPPOL(NTPOL)
  40. SEGMENT NTFAC(NFAC)
  41. SEGMENT KFAK(NFAC)
  42. SEGMENT NAUX(max(2,NFAC))
  43. *SG
  44. * Logique loquaf : pour les faces TRI7 et QUA9, normalement, le
  45. * dernier noeud de la face est unique à la face : il peut donc
  46. * servir de clé de hachage et on peut éviter de vérifier l'égalité
  47. * de tous les autres noeuds lorsque l'on teste l'égalité des faces.
  48. * C'est ce qu'on fait si loquaf=vrai.
  49. *
  50. LOGICAL LOQUAF,LOPT
  51. PARAMETER (LOQUAF=.TRUE.)
  52. character*4 lnoid(1)
  53. * Pour chaque face dans KDFAC, le numéro d'élément associé
  54. * Ne se trouve pas dans CCGEOME, etonnant
  55. INTEGER ITYEL(NTYFAC)
  56. * T3, Q4, T6, Q8, ? , POLY, T7, Q9
  57. DATA ITYEL/4,8,6,10,0,0,7,11,12*0/
  58. data lnoid/'NOID'/
  59.  
  60. *dbg write(ioimp,*) 'coucou envvo2 IIMPI=',IIMPI
  61. c==== LECTURE ET OUVERTURE DU MELEME ==================================
  62. c et eventuelle boucle 10 sur les ojbets meleme elementaires
  63. CALL LIROBJ('MAILLAGE',MELEME,1,IRETOU)
  64. CALL ACTOBJ('MAILLAGE',MELEME,1)
  65. IF (IERR.NE.0) RETURN
  66. call lirmot(lnoid,1,inoid,0)
  67. c on compte le nombre d elements dont les faces sont de type 1 2 3 4
  68. c 7 8 dans NBFAC, on gère séparément les polygones
  69. SEGINI NBFAC
  70. NTPOL = 0
  71. IPT1=MELEME
  72. SEGACT MELEME
  73. DO 10 IOB=1,MAX(1,LISOUS(/1))
  74. IF (LISOUS(/1).NE.0) THEN
  75. IPT1=LISOUS(IOB)
  76. SEGACT IPT1
  77. ENDIF
  78. NBELEM=IPT1.NUM(/2)
  79. c LTEL,LDEL LFAC de CCGEOME remplis par bdata
  80. ILTEL=LTEL(1,IPT1.ITYPEL)
  81. IF (ILTEL.EQ.0) GOTO 10
  82. ILTAD=LTEL(2,IPT1.ITYPEL)
  83. c --- boucle sur les faces de chaque elements ---
  84. DO 13 IF=1,ILTEL
  85. IFT=LDEL(1,ILTAD+IF-1)
  86. IF (IFT.EQ.6) THEN
  87. NTPOL=NTPOL+1
  88. ELSE
  89. NBFAC(IFT)=NBFAC(IFT)+NBELEM
  90. ENDIF
  91. 13 CONTINUE
  92. c --- fin de boucle sur les faces de chaque elements ---
  93. 10 CONTINUE
  94. c write(6,*) 'dimension de NFAC3,4,6,8=',NFAC3,NFAC4,NFAC6,NFAC8
  95.  
  96. c==== CREATION DES FACES ==============================================
  97. SEGINI IPOFAC
  98. DO ITYFAC=1,NTYFAC
  99. NNODE=KDFAC(1,ITYFAC)
  100. IF (NNODE.GT.0) THEN
  101. NFACI=NBFAC(ITYFAC)
  102. SEGINI IFACI
  103. IPOFAC(ITYFAC)=IFACI
  104. ENDIF
  105. ENDDO
  106. SEGINI IPPOL
  107. c NBFAC sert maintenant de compteur
  108. DO ITYFAC=1,NTYFAC
  109. NBFAC(ITYFAC)=0
  110. ENDDO
  111. NTPOL=0
  112. c eventuelle boucle 50 sur les objets meleme elementaires
  113. DO 50 IOB=1,MAX(1,LISOUS(/1))
  114. IF (LISOUS(/1).NE.0) IPT1=LISOUS(IOB)
  115. segact ipt1
  116. NBELEM=IPT1.NUM(/2)
  117. ILTEL=LTEL(1,IPT1.ITYPEL)
  118. IF (ILTEL.EQ.0) GOTO 52
  119. ILTAD=LTEL(2,IPT1.ITYPEL)
  120. c --- boucle 60 sur faces d'1 type d'element ------------------
  121. DO 60 IF=1,ILTEL
  122. ITYFAC=LDEL(1,ILTAD+IF-1)
  123. IAD=LDEL(2,ILTAD+IF-1)
  124. NNODE=KDFAC(1,ITYFAC)
  125. IF (NNODE.GT.0) THEN
  126. IFACI=IPOFAC(ITYFAC)
  127. c --- boucle 66 sur elements ---------------------------------
  128. DO 66 IEL=1,NBELEM
  129. NBFAC(ITYFAC)=NBFAC(ITYFAC)+1
  130. j=NBFAC(ITYFAC)
  131. DO i=1,NNODE
  132. IFACI(i,j)=IPT1.NUM(LFAC(IAD+i-1),IEL)
  133. ENDDO
  134. IFACI(NNODE+1,j)=IPT1.ICOLOR(IEL)
  135. 66 CONTINUE
  136. c --- fin de boucle 66 sur elements ---------------------------------
  137. ENDIF
  138. *SG Avant ce if était après le 52 CONTINUE mais alors ITYFAC etait
  139. * potentiellement non initialise
  140. IF (ITYFAC.EQ.6) THEN
  141. C Polygone
  142. NTPOL = NTPOL+1
  143. IPPOL(NTPOL)= IPT1
  144. ENDIF
  145. 60 CONTINUE
  146. c --- fin de boucle 60 sur faces d'1 type d'element ------------------
  147. 52 CONTINUE
  148. 50 CONTINUE
  149.  
  150. C ======================================================================
  151. C IL FAUT MAINTENANT RETASSER ET CLASSER LES TABLEAUX DES FACES
  152. C PROBLEME ELLES NE SONT PAS FORCEMENT DECRITES DE LA MEME FACON
  153. C SG 20160712 NTFAC sert de cle de hachage, elle est égale à la
  154. C somme des numeros de noeuds de la face
  155. NFAC=0
  156. SEGINI IDXFAC
  157. IDXFAC(1)=NFAC+1
  158. DO ITYFAC=1,NTYFAC
  159. NFAC=NFAC+NBFAC(ITYFAC)
  160. IDXFAC(ITYFAC+1)=NFAC+1
  161. * write(ioimp,*) 'ityfac=',ityfac,' nbfac=',NBFAC(ITYFAC)
  162. ENDDO
  163. SEGINI NTFAC,KFAK
  164. IFAC=0
  165. DO ITYFAC=1,NTYFAC
  166. NNODE=KDFAC(1,ITYFAC)
  167. IF (NNODE.GT.0) THEN
  168. LOPT=LOQUAF.AND.(ITYFAC.EQ.7.OR.ITYFAC.EQ.8)
  169. IFACI=IPOFAC(ITYFAC)
  170. DO I=1,NBFAC(ITYFAC)
  171. IFAC=IFAC+1
  172. IF (LOPT) THEN
  173. NTFAC(IFAC)=IFACI(NNODE,I)
  174. ELSE
  175. DO J=1,NNODE
  176. NTFAC(IFAC)=NTFAC(IFAC)+IFACI(J,I)
  177. ENDDO
  178. ENDIF
  179. KFAK(IFAC)=I
  180. ENDDO
  181. ENDIF
  182. ENDDO
  183. C IL N'Y A PLUS QU'A TRIER ET RETASSER KFAK SUIVANT NTFAC
  184. SEGINI NAUX
  185. DO 300 ITYFAC=1,NTYFAC
  186. IDEB=IDXFAC(ITYFAC)
  187. IFIN=IDXFAC(ITYFAC+1)-1
  188. IF (IFIN.LE.IDEB) GOTO 300
  189. NAUX(1)=IDEB
  190. NAUX(2)=IFIN
  191. *
  192. IZ=2
  193. 208 IZ=IZ-1
  194. IF (IZ.LE.0) GOTO 209
  195. IPB=NAUX(IZ*2-1)
  196. IPH=NAUX(IZ*2)
  197. IF(IPB.GE.IPH) GOTO 208
  198. JPB=IPB-1
  199. JPH=IPH+1
  200. C CALCUL DU PIVOT
  201. NPV=0
  202. * DO 207 J=IPB,IPH
  203. * NPV=NPV+NTFAC(J)
  204. * 207 CONTINUE
  205. * NPV=NPV/(IPH-IPB+1)
  206. NPV=(NTFAC(IPB)+NTFAC(IPH))/2
  207. 242 JPB=JPB+1
  208. IF (JPH.EQ.JPB) GOTO 245
  209. IF (NTFAC(JPB).LE.NPV) GOTO 243
  210. GOTO 242
  211. 243 JPH=JPH-1
  212. IF (JPH.EQ.JPB) GOTO 245
  213. IF (NTFAC(JPH).GE.NPV) GOTO 244
  214. GOTO 243
  215. 244 IAUX=KFAK(JPB)
  216. KFAK(JPB)=KFAK(JPH)
  217. KFAK(JPH)=IAUX
  218. NTAUX=NTFAC(JPB)
  219. NTFAC(JPB)=NTFAC(JPH)
  220. NTFAC(JPH)=NTAUX
  221. GOTO 242
  222. 245 IF (JPB.GE.IPB) GOTO 247
  223. JPB=JPB+1
  224. JPH=JPH+2
  225. GOTO 248
  226. 247 IF (JPH.LE.IPH) GOTO 249
  227. JPB=JPB-2
  228. JPH=JPH-1
  229. GOTO 248
  230. 249 IF (NTFAC(JPB).LE.NPV) GOTO 250
  231. IF (JPH.EQ.IPH) GOTO 251
  232. 252 JPH=JPH+1
  233. GOTO 248
  234. 250 IF (JPB.EQ.IPB) GOTO 252
  235. 251 JPB=JPB-1
  236. 248 IF (JPB.EQ.IPB) GOTO 253
  237. NAUX(2*IZ)=JPB
  238. IZ=IZ+1
  239. 253 IF (JPH.EQ.IPH) GOTO 208
  240. NAUX(2*IZ)=IPH
  241. NAUX(2*IZ-1)=JPH
  242. IZ=IZ+1
  243. GOTO 208
  244. 209 CONTINUE
  245. 300 CONTINUE
  246.  
  247. C ======================================================================
  248. C LES FACES SONT CLASSEES DANS KFAK IL FAUT ELIMINER LES FACES EN DOUBL
  249. C ELLES SONT PAR TYPE LES UNES DERRIERES LES AUTRES LES PLUS HAUTES
  250. C DEVANT
  251. IF (IIMPI.NE.0) WRITE (IOIMP,9111) (KFAK(I),NTFAC(I),I=1,NFAC)
  252. 9111 FORMAT(5(2X,2I6))
  253. DO 400 ITYFAC=1,NTYFAC
  254. IDEB=IDXFAC(ITYFAC)
  255. IFIN=IDXFAC(ITYFAC+1)-1
  256. IF (IFIN.LE.IDEB) GOTO 400
  257. NNODE=KDFAC(1,ITYFAC)
  258. * A cette etape on doit avoir nnode.gt.0
  259. IF (NNODE.LE.0) THEN
  260. CALL ERREUR(5)
  261. RETURN
  262. ENDIF
  263. LOPT=LOQUAF.AND.(ITYFAC.EQ.7.OR.ITYFAC.EQ.8)
  264. IFACI=IPOFAC(ITYFAC)
  265. *
  266. IFINM=IFIN-1
  267. DO 450 I1=IDEB,IFINM
  268. NTI1=NTFAC(I1)
  269. IF (NTI1.EQ.0) GOTO 450
  270. IDEB1=I1+1
  271. DO 460 I2=IDEB1,IFIN
  272. NTI2=NTFAC(I2)
  273. IF (NTI2.EQ.0) GOTO 460
  274. IF (NTI2.NE.NTI1) GOTO 450
  275. IF (.NOT.LOPT) THEN
  276. IR1=KFAK(I1)
  277. IR2=KFAK(I2)
  278. DO 471 J1=1,NNODE
  279. INU=IFACI(J1,IR1)
  280. DO 472 J2=1,NNODE
  281. IF (INU.EQ.IFACI(J2,IR2)) GOTO 471
  282. 472 CONTINUE
  283. GOTO 460
  284. 471 CONTINUE
  285. ENDIF
  286. C DEUX FACES EGALES ON LES SUPPRIMENT (ENVE)
  287. if (icle.eq.0) then
  288. NTFAC(I1)=0
  289. NTFAC(I2)=0
  290. C DEUX FACES EGALES ON EN SUPPRIME UNE (CHAN FACE)
  291. elseif (icle.eq.1) then
  292. * il semble plus astucieux de supprimer la premiere pour le cas qui
  293. * ne doit pas arriver ici ou il y ait plus que deux faces egales
  294. NTFAC(I1)=0
  295. else
  296. call erreur(5)
  297. return
  298. endif
  299. GOTO 450
  300. 460 CONTINUE
  301. 450 CONTINUE
  302. 400 CONTINUE
  303. *
  304. IF (IIMPI.NE.0) WRITE (IOIMP,9111) (KFAK(I),NTFAC(I),I=1,NFAC)
  305.  
  306. SEGINI IPTFAC
  307. NBSOUS=0
  308. NBREF=0
  309. NBSOU2=0
  310. DO 600 ITYFAC=1,NTYFAC
  311. IDEB=IDXFAC(ITYFAC)
  312. IFIN=IDXFAC(ITYFAC+1)-1
  313. * write(ioimp,*) 'ityfac=',ityfac,' ideb=',ideb,' ifin=',ifin
  314. IF (IFIN.LT.IDEB) GOTO 600
  315. NNODE=KDFAC(1,ITYFAC)
  316. * A cette etape on doit avoir nnode.gt.0
  317. IF (NNODE.LE.0) THEN
  318. CALL ERREUR(5)
  319. RETURN
  320. ENDIF
  321. IFACI=IPOFAC(ITYFAC)
  322. NBELEM=0
  323. DO 611 I=IDEB,IFIN
  324. IF (NTFAC(I).NE.0) NBELEM=NBELEM+1
  325. 611 CONTINUE
  326. * write(ioimp,*) 'nbelem=',nbelem,' nnode=',nnode
  327. IF (NBELEM.EQ.0) GOTO 600
  328. NBSOU2=NBSOU2+1
  329. NBNN=NNODE
  330. SEGINI IPT1
  331. IPT1.ITYPEL=ITYEL(ITYFAC)
  332. JAUX=0
  333. DO 612 J=IDEB,IFIN
  334. IF (NTFAC(J).EQ.0) GOTO 612
  335. JAUX=JAUX+1
  336. IPT1.ICOLOR(JAUX)=IFACI(NNODE+1,KFAK(J))
  337. DO 613 I=1,NBNN
  338. IPT1.NUM(I,JAUX)=IFACI(I,KFAK(J))
  339. 613 CONTINUE
  340. 612 CONTINUE
  341. IPTFAC(ITYFAC)=IPT1
  342. * write(ioimp,*) 'ipt1=',ipt1
  343. 600 CONTINUE
  344. * on rajoute les points et les segments qui pouvaient etre dans le
  345. * maillage initial
  346. ipt5=0
  347. segact meleme
  348. ipt6=meleme
  349. do 710 io=1,max(1,lisous(/1))
  350. if (lisous(/1).ne.0) ipt6=lisous(io)
  351. segact ipt6
  352. if (ipt6.itypel.le.3) then
  353. nbsou2=nbsou2+1
  354. ipt5=ipt6
  355. endif
  356. 710 continue
  357. * write(ioimp,*) 'nbsou2=',nbsou2
  358. IF (NBSOU2.EQ.0.AND.NTPOL.EQ.0.AND.inoid.eq.0) CALL ERREUR(26)
  359. IF (NBSOU2.NE.1.OR.NTPOL.GT.0) THEN
  360. NBREF=0
  361. NBSOUS=NBSOU2+NTPOL
  362. NBNN=0
  363. NBELEM=0
  364. SEGINI IPT5
  365. I=0
  366. DO ITYFAC=1,NTYFAC
  367. IPT1=IPTFAC(ITYFAC)
  368. IF (IPT1.NE.0) THEN
  369. I=I+1
  370. IPT5.LISOUS(I)=IPT1
  371. ENDIF
  372. ENDDO
  373. segact meleme
  374. ipt1=meleme
  375. do 711 io=1,max(1,lisous(/1))
  376. if (lisous(/1).ne.0) ipt1=lisous(io)
  377. segact ipt1
  378. if (ipt1.itypel.le.3) then
  379. I=I+1
  380. IPT5.LISOUS(I)=IPT1
  381. endif
  382. 711 continue
  383. DO 720, IO = 1, NTPOL
  384. I = I+1
  385. IPT5.LISOUS(I)=IPPOL(IO)
  386. 720 CONTINUE
  387. if (ipt5.lisous(/1).eq.1) then
  388. ipt6=ipt5
  389. ipt5=ipt6.lisous(1)
  390. segsup ipt6
  391. endif
  392. ELSE
  393. * Cas ou il n'y a qu'un meleme
  394. IPTF=0
  395. DO ITYFAC=1,NTYFAC
  396. IPTF=IPTF+IPTFAC(ITYFAC)
  397. * write(ioimp,*) 'iptf=',iptf
  398. ENDDO
  399. IPTF=IPTF+IPT5
  400. IPT5=IPTF
  401. ENDIF
  402. CALL ACTOBJ('MAILLAGE',IPT5,1)
  403. CALL ECROBJ('MAILLAGE',IPT5)
  404.  
  405. SEGSUP IPTFAC,NAUX,NTFAC,KFAK,IDXFAC,IPPOL
  406. DO ITYFAC=1,NTYFAC
  407. IFACI=IPOFAC(ITYFAC)
  408. IF (IFACI.NE.0) THEN
  409. SEGSUP IFACI
  410. ENDIF
  411. ENDDO
  412. SEGSUP IPOFAC,NBFAC
  413. END
  414.  
  415.  
  416.  

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