Télécharger inclus.eso

Retour à la liste

Numérotation des lignes :

  1. C INCLUS SOURCE BP208322 16/11/18 21:17:43 9177
  2. C
  3. C CE SOUS PROGRAMME EXTRAIT DU PREMIER OBJET L'ENSEMBLE DES ELEMENTS
  4. C STRICTEMENT INCLUS DANS LE DEUXIEME OBJET
  5. C
  6. SUBROUTINE INCLUS
  7. IMPLICIT INTEGER(I-N)
  8. IMPLICIT REAL*8 (A-H,O-Z)
  9. -INC CCREEL
  10. *-
  11. -INC SMELEME
  12. -INC CCOPTIO
  13. -INC SMCOORD
  14. -INC CCGEOME
  15. CHARACTER*4 LEMOT(5)
  16. DATA LEMOT/'BARY','STRI','LARG','VOLU','NOID'/
  17. SEGMENT /FER/(NFI(ITT),MAI(IPP),ITOUR)
  18. SEGMENT ICPR(XCOOR(/1)/(IDIM+1))
  19.  
  20. CALL LIROBJ('MAILLAGE',IPT1,1,IRETOU)
  21. CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU)
  22. IF (IERR.NE.0) RETURN
  23.  
  24. * On fixe les options
  25. * Par défaut : non barycentrique, appartenance large, et erreur si
  26. * maillage vide
  27. IBARY=0
  28. IFRON=0
  29. IVERI=0
  30. DO I=1,3
  31. CALL LIRMOT(LEMOT,5,IRET,0)
  32. IF (IRET.EQ.1) IBARY=1
  33. IF (IRET.EQ.2) IFRON=1
  34. IF (IRET.EQ.3) IFRON=0
  35. IF (IRET.EQ.4) THEN
  36. * Cas de maillage volumique
  37. CALL INCLU3(IPT1,IPT2)
  38. RETURN
  39. ENDIF
  40. IF (IRET.EQ.5) IVERI=1
  41. IF (IRET.EQ.0) GOTO 20
  42. ENDDO
  43.  
  44. 20 CONTINUE
  45.  
  46. C CRITERE D'INCLUSION PAS UTILISE EN 2D...
  47. CALL LIRREE(XCRITT,0,IRET)
  48. XCRITT=0.D0
  49.  
  50. * Cas d'une surface plane.
  51. SEGACT IPT2
  52.  
  53. C Appel éventuel à CONTOUR
  54. IF ((KSURF(IPT2.ITYPEL).EQ.0).AND.(IPT2.LISOUS(/1).EQ.0)) GOTO 25
  55. CALL ECROBJ('MAILLAGE',IPT2)
  56. CALL PRCONT
  57. CALL LIROBJ('MAILLAGE',IPT2,1,IRETOU)
  58. IF (IERR.NE.0) RETURN
  59. SEGACT IPT2
  60.  
  61. 25 CONTINUE
  62. * Transformation éventuelle en SEG2 (SEG3 ou maillage complexe)
  63. IF(IPT2.ITYPEL.EQ.3.OR.IPT2.LISOUS(/1).NE.0)
  64. $ CALL CHANGE (IPT2,2)
  65.  
  66. * Remplissage segment FER
  67. ICLE=0
  68. CALL AVTRSF(IPT2,FER,ICLE)
  69. IF (IERR.NE.0) RETURN
  70.  
  71. * Détermination du plan par son centre de gravité et son vecteur
  72. * normal unitaire (pour changement de variable)
  73. SEGACT MCOORD
  74. SEGINI ICPR
  75. DO 1 I=1,ICPR(/1)
  76. ICPR(I)=0
  77. 1 CONTINUE
  78. XGRAV=0.
  79. YGRAV=0.
  80. ZGRAV=0.
  81. DO 30 I=MAI(1)+1,MAI(ITOUR+1)
  82. IREF=NFI(I)*(IDIM+1)-IDIM
  83. XGRAV=XGRAV+XCOOR(IREF)
  84. YGRAV=YGRAV+XCOOR(IREF+1)
  85. ZGRAV=ZGRAV+XCOOR(IREF+2)
  86. 30 CONTINUE
  87. XGRAV=XGRAV/FLOAT(MAI(ITOUR+1)-MAI(1))
  88. YGRAV=YGRAV/FLOAT(MAI(ITOUR+1)-MAI(1))
  89. ZGRAV=ZGRAV/FLOAT(MAI(ITOUR+1)-MAI(1))
  90. IF (IDIM.EQ.2) ZGRAV=0.
  91. XNORM=0.D0
  92. YNORM=0.D0
  93. ZNORM=0.D0
  94. DO 42 IT=1,ITOUR
  95. IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM
  96. XV1=XCOOR(IREF)-XGRAV
  97. YV1=XCOOR(IREF+1)-YGRAV
  98. ZV1=XCOOR(IREF+2)-ZGRAV
  99. IF (IDIM.EQ.2) ZV1=0.
  100. DO 421 I=MAI(IT-1+1)+1,MAI(IT+1)
  101. IREF=(NFI(I))*(IDIM+1)-IDIM
  102. XV2=XCOOR(IREF)-XGRAV
  103. YV2=XCOOR(IREF+1)-YGRAV
  104. ZV2=0.
  105. IF (IDIM.GE.3) ZV2=XCOOR(IREF+2)-ZGRAV
  106. XNORM=XNORM+YV1*ZV2-ZV1*YV2
  107. YNORM=YNORM+ZV1*XV2-XV1*ZV2
  108. ZNORM=ZNORM+XV1*YV2-YV1*XV2
  109. XV1=XV2
  110. YV1=YV2
  111. ZV1=ZV2
  112. 421 CONTINUE
  113. 42 CONTINUE
  114. DNORM=SQRT(XNORM**2+YNORM**2+ZNORM**2)
  115. XNORM=XNORM/DNORM
  116. YNORM=YNORM/DNORM
  117. ZNORM=ZNORM/DNORM
  118.  
  119. MELEME=IPT1
  120. SEGACT MELEME
  121. NBS =LISOUS(/1)
  122. NBSOUS=NBS
  123. IPT8=0
  124. IF (NBSOUS.NE.0) THEN
  125. NBREF=0
  126. NBNN=0
  127. NBELEM=0
  128. SEGINI IPT8
  129. IOB1=0
  130. ENDIF
  131.  
  132. IF(IBARY.EQ.0) THEN
  133. * Eligibilité d'après chaque point de l'élément
  134. DO 7001 IOB=1,MAX(1,NBS)
  135. IF (NBS.NE.0) THEN
  136. IPT1=LISOUS(IOB)
  137. SEGACT IPT1
  138. ENDIF
  139. NBELEM=IPT1.NUM(/2)
  140. NBNN=IPT1.NUM(/1)
  141. NBREF=0
  142. NBSOUS=0
  143. SEGINI IPT4
  144. IPT4.ITYPEL=IPT1.ITYPEL
  145. DO 2 I=1,NBNN
  146. DO 2222 J=1,NBELEM
  147. IPT4.NUM(I,J)=0
  148. 2222 CONTINUE
  149. 2 CONTINUE
  150.  
  151. ICO=0
  152. DO 10 J=1,NBELEM
  153. DO 11 I=1,NBNN
  154. IP=IPT1.NUM(I,J)
  155. IF (ICPR(IP).EQ.-1) GOTO 10
  156. IF (ICPR(IP).EQ.1) GOTO 11
  157. IREF=IP*(IDIM+1)-IDIM
  158. XP=XCOOR(IREF)
  159. YP=XCOOR(IREF+1)
  160. ZP=XCOOR(IREF+2)
  161. IF (IDIM.EQ.2) ZP=0.
  162. ANG=0.D0
  163. DO 12 IT=1,ITOUR
  164. IF (IP.EQ.NFI(MAI(IT+1)))THEN
  165. IF(IFRON.EQ.0) GO TO 13
  166. GO TO 18
  167. ENDIF
  168. IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM
  169. XV1=XCOOR(IREF)-XP
  170. YV1=XCOOR(IREF+1)-YP
  171. ZV1=XCOOR(IREF+2)-ZP
  172. IF (IDIM.EQ.2) ZV1=0.
  173. DO 120 M=MAI(IT-1+1)+1,MAI(IT+1)
  174. IF (IP.EQ.NFI(M))THEN
  175. IF(IFRON.EQ.0) GO TO 13
  176. GO TO 18
  177. ENDIF
  178. IREF=NFI(M)*(IDIM+1)-IDIM
  179. XV2=XCOOR(IREF)-XP
  180. YV2=XCOOR(IREF+1)-YP
  181. ZV2=XCOOR(IREF+2)-ZP
  182. IF (IDIM.EQ.2) ZV2=0.
  183. V1=XNORM*(YV1*ZV2-ZV1*YV2)+YNORM*(ZV1*XV2-XV1
  184. $ *ZV2)+ZNORM*(XV1*YV2-YV1*XV2)
  185. V2=XV1*XV2+YV1*YV2+ZV1*ZV2
  186. IF( IFRON.EQ.1) THEN
  187. V6 = ( XV1 -XV2)*( XV1 -XV2) + ( YV1 -YV2)* (
  188. $ YV1 -YV2)+ ( ZV1-ZV2)* ( ZV1-ZV2)
  189. IF( V2.LE.0. . OR. ( V2 . LT . V6* 1.D-4))
  190. $ THEN
  191. V4 =SQRT (XV1* XV1 + YV1 * YV1 + ZV1 * ZV1
  192. $ )
  193. V5 =SQRT (XV2* XV2 + YV2 * YV2 + ZV2
  194. $ * ZV2 )
  195. V3 = ((YV1*ZV2-ZV1*YV2)*(YV1*ZV2-ZV1*YV2)
  196. $ +(ZV1*XV2-XV1*ZV2)*(ZV1*XV2-XV1*ZV2)
  197. $ +(XV1*YV2-YV1*XV2) *(XV1*YV2-YV1*XV2)
  198. $ ) **0.5D0
  199. IF( V4* 1.D4. LT. V5 . OR . V5*1.D4 .LT.V4
  200. $ ) GO TO 18
  201. IF( V3*1.D4 . LT . ( ABS ( V2) ) )
  202. # go to 18
  203. ENDIF
  204. ENDIF
  205. IF (V1.EQ.0.D0.AND.V2.EQ.0.D0) GOTO 121
  206. ANG=ANG+ATAN2(V1,V2)
  207. 121 CONTINUE
  208. XV1=XV2
  209. YV1=YV2
  210. ZV1=ZV2
  211. 120 CONTINUE
  212. 12 CONTINUE
  213. IF (ABS(ANG).GT.XPI) GOTO 13
  214. 18 ICPR(IP)=-1
  215. GOTO 10
  216. 13 ICPR(IP)=1
  217. 11 CONTINUE
  218. * Cet élément est ok, on le met de côté.
  219. DO 14 L=1,NBNN
  220. IPT4.NUM(L,J)=IPT1.NUM(L,J)
  221. 14 CONTINUE
  222. IPT4.ICOLOR(J)=IPT1.ICOLOR(J)
  223. ICO=ICO+1
  224. 10 CONTINUE
  225.  
  226. NBELEM=ICO
  227. IF (NBELEM.EQ.0) GOTO 7001
  228.  
  229. * sauvegarde partition ok
  230. SEGINI IPT3
  231. IPT3.ITYPEL=IPT1.ITYPEL
  232. JJ=0
  233. DO 21 J=1,IPT4.NUM(/2)
  234. IF (IPT4.NUM(1,J).EQ.0) GOTO 21
  235. JJ=JJ+1
  236. IPT3.ICOLOR(JJ)=IPT4.ICOLOR(J)
  237. DO 22 I=1,NBNN
  238. IPT3.NUM(I,JJ)=IPT4.NUM(I,J)
  239. 22 CONTINUE
  240. 21 CONTINUE
  241. SEGSUP IPT4
  242.  
  243. IF (NBS.NE.0) THEN
  244. SEGDES IPT1,IPT3
  245. IOB1=IOB1+1
  246. IPT8.LISOUS(IOB1)=IPT3
  247. ENDIF
  248. 7001 CONTINUE
  249. ELSE
  250. * Eligibilité d'après la position du barycentre
  251. DO 7051 IOB=1,MAX(1,NBS)
  252. IF (NBS.NE.0) THEN
  253. IPT1=LISOUS(IOB)
  254. SEGACT IPT1
  255. ENDIF
  256. NBELEM=IPT1.NUM(/2)
  257. NBNN=IPT1.NUM(/1)
  258. NBREF=0
  259. NBSOUS=0
  260. SEGINI IPT4
  261. IPT4.ITYPEL=IPT1.ITYPEL
  262. DO 52 I=1,NBNN
  263. DO 521 J=1,NBELEM
  264. IPT4.NUM(I,J)=0
  265. 521 CONTINUE
  266. 52 CONTINUE
  267. ICO=0
  268. DO 60 J=1,NBELEM
  269. XP=0.D0
  270. YP=0.D0
  271. ZP=0.D0
  272. DO 61 I=1,NBNN
  273. IP=IPT1.NUM(I,J)
  274. IREF=IP*(IDIM+1)-IDIM
  275. XP=XCOOR(IREF) +XP
  276. YP=XCOOR(IREF+1)+YP
  277. ZP=XCOOR(IREF+2) +ZP
  278. 61 CONTINUE
  279. IF (IDIM.EQ.2) ZP=0.D0
  280. XP = XP / FLOAT(NBNN)
  281. YP = YP / FLOAT(NBNN)
  282. ZP = ZP / FLOAT(NBNN)
  283. ANG=0.D0
  284. DO 62 IT=1,ITOUR
  285. IREF=NFI(MAI(IT+1))*(IDIM+1)-IDIM
  286. XV1=XCOOR(IREF)-XP
  287. YV1=XCOOR(IREF+1)-YP
  288. ZV1=XCOOR(IREF+2)-ZP
  289. IF (IDIM.EQ.2) ZV1=0.
  290. DO 621 M=MAI(IT-1+1)+1,MAI(IT+1)
  291. IREF=NFI(M)*(IDIM+1)-IDIM
  292. XV2=XCOOR(IREF)-XP
  293. YV2=XCOOR(IREF+1)-YP
  294. ZV2=XCOOR(IREF+2)-ZP
  295. IF (IDIM.EQ.2) ZV2=0.
  296. V1=XNORM*(YV1*ZV2-ZV1*YV2)+YNORM*(ZV1*XV2-XV1*ZV2)
  297. $ +ZNORM*(XV1*YV2-YV1*XV2)
  298. V2=XV1*XV2+YV1*YV2+ZV1*ZV2
  299. IF (V1.EQ.0.D0.AND.V2.EQ.0.D0) GOTO 171
  300. ANG=ANG+ATAN2(V1,V2)
  301. 171 CONTINUE
  302. XV1=XV2
  303. YV1=YV2
  304. ZV1=ZV2
  305. 621 CONTINUE
  306. 62 CONTINUE
  307. IF (ABS(ANG).GT.XPI) GOTO 63
  308. GOTO 60
  309. 63 CONTINUE
  310. * Cet élément est ok, on le met de côté.
  311. DO 64 L=1,NBNN
  312. IPT4.NUM(L,J)=IPT1.NUM(L,J)
  313. 64 CONTINUE
  314. IPT4.ICOLOR(J)=IPT1.ICOLOR(J)
  315. ICO=ICO+1
  316. 60 CONTINUE
  317.  
  318. NBELEM=ICO
  319. IF (NBELEM.EQ.0) GOTO 7051
  320.  
  321. * sauvegarde partition ok
  322. SEGINI IPT3
  323. IPT3.ITYPEL=IPT1.ITYPEL
  324. JJ=0
  325. DO 71 J=1,IPT4.NUM(/2)
  326. IF (IPT4.NUM(1,J).EQ.0) GOTO 71
  327. JJ=JJ+1
  328. IPT3.ICOLOR(JJ)=IPT4.ICOLOR(J)
  329. DO 72 I=1,NBNN
  330. IPT3.NUM(I,JJ)=IPT4.NUM(I,J)
  331. 72 CONTINUE
  332. 71 CONTINUE
  333. SEGSUP IPT4
  334. IF (NBS.NE.0) THEN
  335. SEGDES IPT1,IPT3
  336. IOB1=IOB1+1
  337. IPT8.LISOUS(IOB1)=IPT3
  338. ENDIF
  339. 7051 CONTINUE
  340. ENDIF
  341.  
  342. * Ecriture du maillage résultat et sortie
  343. * ---------------------------------------
  344. IF (((NBS.EQ.0).AND.(NBELEM.EQ.0)).OR.
  345. & ((NBS.NE.0).AND.(IOB1.EQ.0))) THEN
  346. * on n'a rien trouvé
  347. IF (IVERI.EQ.0) THEN
  348. * Tache impossible. Probablement données erronées
  349. CALL ERREUR(26)
  350. RETURN
  351. ELSE
  352. * on renvoie un maillage vide
  353. GOTO 7006
  354. ENDIF
  355. ENDIF
  356. IF (NBS.EQ.0) THEN
  357. * Si maillage simple en entrée
  358. * On écrit directement IPT3
  359. GOTO 7002
  360. ELSE
  361. * Si maillage complexe en entrée
  362. * IOB1 : nombre de partitions retenues
  363. IF (IOB1.EQ.NBS) THEN
  364. GOTO 7003
  365. ELSEIF (IOB1.EQ.1) THEN
  366. GOTO 7005
  367. ELSE
  368. GOTO 7004
  369. ENDIF
  370. ENDIF
  371.  
  372. * Maillage complexe, avec moins de partitions que l'objet entré
  373. 7004 CONTINUE
  374. NBSOUS=IOB1
  375. NBREF=0
  376. NBNN=0
  377. NBELEM=0
  378. SEGINI IPT3
  379. DO IOB=1,NBSOUS
  380. IPT3.LISOUS(IOB)=IPT8.LISOUS(IOB)
  381. ENDDO
  382. SEGSUP IPT8
  383. GOTO 7002
  384.  
  385. * Maillage simple, avec moins de partitions que l'objet entré
  386. 7005 CONTINUE
  387. IPT3=IPT8.LISOUS(1)
  388. SEGSUP IPT8
  389. GOTO 7002
  390.  
  391. * Maillage avec autant de partitions que l'objet entré
  392. * (inutile de recréer un nouveau maillage)
  393. 7003 CONTINUE
  394. IPT3=IPT8
  395. GOTO 7002
  396.  
  397. * Maillage vide
  398. 7006 CONTINUE
  399. NBSOUS=0
  400. NBREF=0
  401. NBNN=0
  402. NBELEM=0
  403. SEGINI IPT3
  404. IF (NBS.NE.0) SEGSUP IPT8
  405.  
  406. 7002 CONTINUE
  407. SEGDES MELEME,IPT3
  408. SEGSUP ICPR,FER
  409. CALL ECROBJ('MAILLAGE',IPT3)
  410. SEGDES IPT3
  411.  
  412. RETURN
  413. END
  414.  
  415.  
  416.  
  417.  
  418.  
  419.  
  420.  
  421.  
  422.  
  423.  
  424.  
  425.  
  426.  
  427.  

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