Télécharger crcoup.eso

Retour à la liste

Numérotation des lignes :

  1. C CRCOUP SOURCE BP208322 16/11/18 21:16:04 9177
  2. C CALCUL D'UN OBJET COUPE VIRTUEL
  3. C
  4. C CRCOUP REMPLACE MELEME PAR L'OBJET COUPE
  5. C MCOUP ET LE DERNIER COMPOSANT DE MELEME REMPLACE DECRIVENT
  6. C L'INTERSECTION
  7. C IL FAUDRA METTRE UNE INDICATION AUX TRIANGLES SE TROUVANT DANS LE PLA
  8. C DE COUPE POUR REUSSIR LES PARTIES VUES PARTIES CACHEES
  9. C
  10. C SG 2016/07/18 : traitement du cas des faces tri7/qua9
  11. C
  12. SUBROUTINE CRCOUP(JOEIL,ICOUP1,ICOUP2,ICOUP3,MELEME,MCOUP,VCHCA,
  13. # MELEM2,MCOU2,mchamo,isect)
  14. IMPLICIT INTEGER(I-N)
  15. -INC CCGEOME
  16. -INC CCREEL
  17. -INC CCOPTIO
  18. -INC SMCOORD
  19. -INC SMELEME
  20. -INC SMCHAML
  21. REAL*8 XN,YN,ZN,SIGA,SIGB,SIGC,SIGRAV
  22. * sg
  23. REAL*8 xo,yo,zo,xct,yct,zct,xc2,yc2,zc2,xc3,yc3,zc3
  24. REAL*8 xv,yv,zv,xw,yw,zw,xp,yp,zp,xgrav,ygrav,zgrav
  25. REAL*8 xa,ya,za,xb,yb,zb,xc,yc,zc
  26. REAL*8 rn,sig,chgrav
  27. * sg
  28. SEGMENT VCHCA(0)
  29. SEGMENT ITR(0)
  30. SEGMENT ITRT(0)
  31. SEGMENT XTR(0)
  32. SEGMENT XTRT(0)
  33. SEGMENT IQUAT(0)
  34. SEGMENT XQUAT(0)
  35. SEGMENT ISEGM(0)
  36. SEGMENT XSEGM(NBPOIN)
  37. SEGMENT MCOUP(0)
  38. SEGMENT IJGRAV(2,100)
  39. DIMENSION JTT(9),ISEG(9)
  40. LOGICAL VOLUM
  41. logical lquaf
  42. *
  43. *dbg write(ioimp,*) 'coucou crcoup isect=',isect
  44. n2ptel=0
  45. mcham=mchamo
  46. MEDSAU=0
  47. SEGINI ITR,ISEGM,MCOUP,IJGRAV,ITRT,IQUAT
  48. SEGACT MCOORD
  49. NBPOIN=XCOOR(/1)/4
  50. if (mcham.ne.0) segini xtr,xtrt,xquat,xsegm
  51. IREF=(JOEIL-1)*4
  52. XO=XCOOR(IREF+1)
  53. YO=XCOOR(IREF+2)
  54. ZO=XCOOR(IREF+3)
  55. IREF=(ICOUP1-1)*4
  56. XCT=XCOOR(IREF+1)
  57. YCT=XCOOR(IREF+2)
  58. ZCT=XCOOR(IREF+3)
  59. IREF=(ICOUP2-1)*4
  60. XC2=XCOOR(IREF+1)
  61. YC2=XCOOR(IREF+2)
  62. ZC2=XCOOR(IREF+3)
  63. IREF=(ICOUP3-1)*4
  64. XC3=XCOOR(IREF+1)
  65. YC3=XCOOR(IREF+2)
  66. ZC3=XCOOR(IREF+3)
  67. C NORMALE AU PLAN :
  68. XV=XC2-XCT
  69. YV=YC2-YCT
  70. ZV=ZC2-ZCT
  71. XW=XC3-XCT
  72. YW=YC3-YCT
  73. ZW=ZC3-ZCT
  74. XN=YV*ZW-ZV*YW
  75. YN=ZV*XW-XV*ZW
  76. ZN=XV*YW-YV*XW
  77. RN=XN**2+YN**2+ZN**2
  78. * IF (RN.EQ.0.) CALL ERREUR(21)
  79. IF (ABS(RN).LT.(XPETIT**0.25D0)) CALL ERREUR(21)
  80. IF (IERR.NE.0) RETURN
  81. RN=SQRT(RN)
  82. XN=XN/RN
  83. YN=YN/RN
  84. ZN=ZN/RN
  85. SIG=(XO-XCT)*XN+(YO-YCT)*YN+(ZO-ZCT)*ZN
  86. IF (SIG.LE.0.) THEN
  87. XN=-XN
  88. YN=-YN
  89. ZN=-ZN
  90. ENDIF
  91. 5000 CONTINUE
  92. SEGACT MELEME
  93. NBSOU=LISOUS(/1)
  94. NBSOUS=MAX(2,NBSOU+1)
  95. NBREF=0
  96. if (mcham.ne.0) nbref=nbsous
  97. NBNN=0
  98. NBELEM=0
  99. SEGINI IPT8
  100. IPT1=MELEME
  101. ISU=0
  102. * prob optimiseur
  103. C melva1=meleme
  104. C melva2=meleme
  105. DO 10 ISOUS=1,MAX(1,NBSOU)
  106. IF (NBSOU.NE.0) THEN
  107. IPT1=LISOUS(ISOUS)
  108. if (mcham.ne.0) then
  109. melva1=lisref(isous)
  110. segact melva1
  111. lva1=melva1.velche(/1)
  112. endif
  113. SEGACT IPT1
  114. ENDIF
  115. NBREF=0
  116. NBSOUS=0
  117. NBNN=IPT1.NUM(/1)
  118. NBELEM=IPT1.NUM(/2)
  119. SEGINI IPT2
  120. IPT2.ITYPEL=IPT1.ITYPEL
  121. if (mcham.ne.0) then
  122. n1ptel=nbnn
  123. n1el=nbelem
  124. n2ptel=0
  125. n2el=0
  126. segini melva2
  127. segact melva1
  128. idi2= melva1.velche(/2)
  129. endif
  130. JELN=0
  131. DO 20 JEL=1,NBELEM
  132. DO 25 INOEU=1,NBNN
  133. IPT=IPT1.NUM(INOEU,JEL)
  134. IREF=4*(IPT-1)
  135. XP=XCOOR(IREF+1)
  136. YP=XCOOR(IREF+2)
  137. ZP=XCOOR(IREF+3)
  138. SIG=(XP-XCT)*XN+(YP-YCT)*YN+(ZP-ZCT)*ZN
  139. IF (SIG.GT.0.) GOTO 30
  140. 25 CONTINUE
  141. * tous les points sont du bon cote du plan de coupe
  142. if (isect.ne.0) goto 20
  143. JELN=JELN+1
  144. DO 27 INOEU=1,NBNN
  145. IPT2.NUM(INOEU,JELN)=IPT1.NUM(INOEU,JEL)
  146. if (mcham.ne.0) melva2.velche(INOEU,JELN)=
  147. $ melva1.velche(min(INOEU,lva1),JEL)
  148. 27 CONTINUE
  149. IPT2.ICOLOR(JELN)=IPT1.ICOLOR(JEL)
  150. GOTO 20
  151. 30 CONTINUE
  152. VOLUM=KSURF(IPT1.ITYPEL).NE.IPT1.ITYPEL
  153. * IF (VOLUM) THEN
  154. IPT=IPT1.NUM(1,JEL)
  155. IREF=(IPT-1)*4
  156. XGRAV=XCOOR(IREF+1)
  157. YGRAV=XCOOR(IREF+2)
  158. ZGRAV=XCOOR(IREF+3)
  159. IF (VCHCA.NE.0) VGRAV=VCHCA(IPT)
  160. if (mcham.ne.0) chgrav=melva1.velche(1,min(JEL,idi2))
  161. DO 21 INOEU=2,NBNN
  162. IPT=IPT1.NUM(INOEU,JEL)
  163. IREF=(IPT-1)*4
  164. XGRAV=XGRAV+XCOOR(IREF+1)
  165. YGRAV=YGRAV+XCOOR(IREF+2)
  166. ZGRAV=ZGRAV+XCOOR(IREF+3)
  167. IF (VCHCA.NE.0) VGRAV=VGRAV+VCHCA(IPT)
  168. if (mcham.ne.0) chgrav=chgrav+
  169. $ melva1.velche(min(INOEU,lva1),min(idi2,JEL))
  170. 21 CONTINUE
  171. XGRAV=XGRAV/NBNN
  172. YGRAV=YGRAV/NBNN
  173. ZGRAV=ZGRAV/NBNN
  174. C ON PREND COMME POINT DU PLAN LA PROJECTION DU CENTRE GRAVITE
  175. SIGRAV=(XGRAV-XCT)*XN+(YGRAV-YCT)*YN+(ZGRAV-ZCT)*ZN
  176. XCT=XGRAV-SIGRAV*XN
  177. YCT=YGRAV-SIGRAV*YN
  178. ZCT=ZGRAV-SIGRAV*ZN
  179. SIGRAV=(XGRAV-XCT)*XN+(YGRAV-YCT)*YN+(ZGRAV-ZCT)*ZN
  180. IF (ABS(SIGRAV).LE.1D-30) SIGRAV=1D-30
  181. NJGRAV=0
  182. IF (VCHCA.NE.0) VGRAV=VGRAV/NBNN
  183. if (mcham.ne.0) CHGRAV=CHGRAV/NBNN
  184. * ENDIF
  185. C DECOMPOSITION EN FACE PUIS EN TRIANGLE
  186. NBFAC=LTEL(1,IPT1.ITYPEL)
  187. IAD=LTEL(2,IPT1.ITYPEL)-1
  188. IF (NBFAC.EQ.0) GOTO 20
  189. DO 161 IFAC=1,NBFAC
  190. ITYP=LDEL(1,IAD+IFAC)
  191. * SG Pour les TRI7/QUA9, il y a un traitement particulier pour ne pas
  192. * tracer les segments relies au noeud central. Il prend en compte le
  193. * fait que le decoupage est barycentrique
  194. lquaf=(ityp.eq.7.OR.ityp.eq.8)
  195. JAD=LDEL(2,IAD+IFAC)-1
  196. NPFAC=KDFAC(1,ITYP)
  197. IDEP=KDFAC(2,ITYP)
  198. IFEP=IDEP+3*(KDFAC(3,ITYP)-1)
  199. * write(ioimp,*) 'ifac,npfac=',ifac,npfac
  200. DO 160 ITRIAN=IDEP,IFEP,3
  201. kfa=kfac(itrian)
  202. kfb=kfac(itrian+1)
  203. kfc=kfac(itrian+2)
  204. JFA=JAD+KFAC(ITRIAN)
  205. JFB=JAD+KFAC(ITRIAN+1)
  206. JFC=JAD+KFAC(ITRIAN+2)
  207. ISEG(1)=0
  208. ISEG(2)=0
  209. ISEG(3)=0
  210. ISEG(4)=0
  211. ISEG1=0
  212. ISEG2=0
  213. ISEG3=0
  214. if (lquaf) then
  215. if (KFA.NE.NPFAC.AND.KFB.NE.NPFAC) ISEG1=1
  216. if (KFB.NE.NPFAC.AND.KFC.NE.NPFAC) ISEG2=1
  217. if (KFA.NE.NPFAC.AND.KFC.NE.NPFAC) ISEG3=1
  218. else
  219. * imod1 = MOD(JFB-JFA+NPFAC,NPFAC)
  220. * imod2 = MOD(JFC-JFB+NPFAC,NPFAC)
  221. * imod3 = MOD(JFA-JFC+NPFAC,NPFAC)
  222. IF (MOD(JFB-JFA+NPFAC,NPFAC).EQ.1) ISEG1=1
  223. IF (MOD(JFB-JFA+NPFAC,NPFAC).EQ.NPFAC-1) ISEG1=1
  224. IF (MOD(JFC-JFB+NPFAC,NPFAC).EQ.1) ISEG2=1
  225. IF (MOD(JFC-JFB+NPFAC,NPFAC).EQ.NPFAC-1) ISEG2=1
  226. IF (MOD(JFA-JFC+NPFAC,NPFAC).EQ.1) ISEG3=1
  227. IF (MOD(JFA-JFC+NPFAC,NPFAC).EQ.NPFAC-1) ISEG3=1
  228. * write(ioimp,*) ' itrian,kfa,kfb,kfc=',itrian
  229. * $ ,kfac(itrian),kfac(itrian+1),kfac(itrian+2)
  230. * write(ioimp,*) ' iseg1,iseg2,iseg3=',iseg1,iseg2,iseg3
  231. * write(ioimp,*) ' imod1,imod2,imod3=',imod1,imod2,imod3
  232. * write(ioimp,*) ' imod1p,imod2p,imod3p=',imod1p,imod2p
  233. * $ ,imod3p
  234. endif
  235. IAFA=LFAC(JFA)
  236. IBFA=LFAC(JFB)
  237. ICFA=LFAC(JFC)
  238. IA=IPT1.NUM(IAFA,JEL)
  239. IB=IPT1.NUM(IBFA,JEL)
  240. IC=IPT1.NUM(ICFA,JEL)
  241. IREF=4*(IA-1)
  242. XA=XCOOR(IREF+1)
  243. YA=XCOOR(IREF+2)
  244. ZA=XCOOR(IREF+3)
  245. IREF=4*(IB-1)
  246. XB=XCOOR(IREF+1)
  247. YB=XCOOR(IREF+2)
  248. ZB=XCOOR(IREF+3)
  249. IREF=4*(IC-1)
  250. XC=XCOOR(IREF+1)
  251. YC=XCOOR(IREF+2)
  252. ZC=XCOOR(IREF+3)
  253. C TESTONS SI LES POINTS SONT DU BON COTE
  254. SIGA=(XA-XCT)*XN+(YA-YCT)*YN+(ZA-ZCT)*ZN
  255. IF (ABS(SIGA).LE.1D-30) SIGA=1D-30
  256. SIGB=(XB-XCT)*XN+(YB-YCT)*YN+(ZB-ZCT)*ZN
  257. IF (ABS(SIGB).LE.1D-30) SIGB=1D-30
  258. SIGC=(XC-XCT)*XN+(YC-YCT)*YN+(ZC-ZCT)*ZN
  259. IF (ABS(SIGC).LE.1D-30) SIGC=1D-30
  260. * if (isect.eq.1) goto 350
  261. IF (SIGA.GT.0..OR.SIGB.GT.0..OR.SIGC.GT.0.) GOTO 103
  262. IF (VOLUM) THEN
  263. IA1=MIN(IA,IB,IC)
  264. IB1=MAX(IA,IB,IC)
  265. IC1=IA+IB+IC
  266. DO 102 IDJF=ITRT(/1)-4,1,-5
  267. ITR1=MIN(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
  268. ITR2=MAX(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
  269. ITR3=ITRT(IDJF)+ITRT(IDJF+1)+ITRT(IDJF+2)
  270. IF (IA1.NE.ITR1) GOTO 102
  271. IF (IB1.NE.ITR2) GOTO 102
  272. IF (IC1.NE.ITR3) GOTO 102
  273. * WRITE (6,*) ' TRIANGLE EN DOUBLE CAS 1 ',IA,IB,IC
  274. C LE TRIANGLE EST DECLARE NON TRACABLE
  275. ITRT(IDJF+4)=MOD(ITRT(IDJF+4),16)+16
  276. GOTO 301
  277. 102 CONTINUE
  278. ENDIF
  279. if (isect.eq.0) then
  280. * WRITE (6,*) ' NOUVEAU TRIANGLE ',IA,IB,IC
  281. ITRT(**)=IA
  282. ITRT(**)=IB
  283. ITRT(**)=IC
  284. ITRT(**)=IPT1.ICOLOR(JEL)
  285. ITRT(**)=ISEG1+2*ISEG2+4*ISEG3
  286. if (mcham.ne.0) then
  287. xtrt(**)=melva1.velche(min(lva1,iafa),min(idi2,jel))
  288. xtrt(**)=melva1.velche(min(lva1,ibfa),min(idi2,jel))
  289. xtrt(**)=melva1.velche(min(lva1,icfa),min(idi2,jel))
  290. endif
  291. endif
  292. 301 CONTINUE
  293. GOTO 101
  294. 103 CONTINUE
  295. IF (SIGA.GT.0..AND.SIGB.GT.0..AND.SIGC.GT.0.) GOTO 101
  296. IPOS=0
  297. IF (SIGA.LE.0.) THEN
  298. IPOS=IPOS+1
  299. JTT(IPOS)=IA
  300. ISEG(IPOS)=ISEG1
  301. if(mcham.ne.0) xsegm(ia)=
  302. $ melva1.velche(min(lva1,iafa),min(idi2,jel))
  303. ENDIF
  304. IF (SIGA*SIGB.LT.0..AND.(.NOT.(SIGA.LE.0..AND.SIGB.LE.0.)))
  305. $ THEN
  306. IA1=MIN(IA,IB)
  307. IB1=MAX(IA,IB)
  308. DO 104 NSEGM=ISEGM(/1)-2,1,-3
  309. IF (IA1.NE.ISEGM(NSEGM)) GOTO 104
  310. IF (IB1.NE.ISEGM(NSEGM+1)) GOTO 104
  311. ISUP=ISEGM(NSEGM+2)
  312. if (mcham.ne.0) xsegm(isup)=
  313. $ (melva1.velche(min(lva1,iafa),min(idi2,jel))*sigb
  314. $ -melva1.velche(min(lva1,ibfa),min(idi2,jel))*siga)
  315. $ /(sigb-siga)
  316. * WRITE (6,*) ' POINT EN DOUBLE ',IA1,IB1,ISUP
  317. GOTO 105
  318. 104 CONTINUE
  319. NBPTS=XCOOR(/1)/4+1
  320. SEGADJ MCOORD
  321. XCOOR((NBPTS-1)*4+1)=(XA*SIGB-XB*SIGA)/(SIGB-SIGA)
  322. XCOOR((NBPTS-1)*4+2)=(YA*SIGB-YB*SIGA)/(SIGB-SIGA)
  323. XCOOR((NBPTS-1)*4+3)=(ZA*SIGB-ZB*SIGA)/(SIGB-SIGA)
  324. XCOOR(NBPTS*(IDIM+1))=0.D0
  325. IF (VCHCA.NE.0)
  326. $ VCHCA(**)=(VCHCA(IA)*SIGB-VCHCA(IB)*SIGA)
  327. $ /(SIGB-SIGA)
  328. ISUP=NBPTS
  329. ISEGM(**)=IA1
  330. ISEGM(**)=IB1
  331. ISEGM(**)=ISUP
  332. if (mcham.ne.0)
  333. > xsegm(**)=(melva1.velche(min(lva1,iafa),min(idi2,jel))*sigb-
  334. > melva1.velche(min(lva1,ibfa),min(idi2,jel))*siga)/(sigb-siga)
  335. 105 CONTINUE
  336. IPOS=IPOS+1
  337. JTT(IPOS)=ISUP
  338. IF (SIGB.LE.0.) THEN
  339. ISEG(IPOS)=ISEG1
  340. ELSE
  341. ISEG(IPOS)=1
  342. ENDIF
  343. ENDIF
  344. IF (SIGB.LE.0.) THEN
  345. IPOS=IPOS+1
  346. JTT(IPOS)=IB
  347. ISEG(IPOS)=ISEG2
  348. if (mcham.ne.0)
  349. $ xsegm(ib)=melva1.velche(min(lva1,ibfa),min(idi2,jel))
  350. ENDIF
  351. IF (SIGB*SIGC.LT.0..AND.(.NOT.(SIGB.LE.0..AND.SIGC.LE.0.)))
  352. $ THEN
  353. IB1=MIN(IB,IC)
  354. IC1=MAX(IB,IC)
  355. DO 106 NSEGM=ISEGM(/1)-2,1,-3
  356. IF (IB1.NE.ISEGM(NSEGM)) GOTO 106
  357. IF (IC1.NE.ISEGM(NSEGM+1)) GOTO 106
  358. ISUP=ISEGM(NSEGM+2)
  359. if (mcham.ne.0)
  360. $ xsegm(isup)=(
  361. $ melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigc
  362. $ -melva1.velche(min(lva1,icfa),min(idi2,jel))*sigb)
  363. $ /(sigc-sigb)
  364. * WRITE (6,*) ' POINT EN DOUBLE ',IB1,IC1,ISUP
  365. GOTO 107
  366. 106 CONTINUE
  367. NBPTS=XCOOR(/1)/4+1
  368. SEGADJ MCOORD
  369. XCOOR((NBPTS-1)*4+1)=(XB*SIGC-XC*SIGB)/(SIGC-SIGB)
  370. XCOOR((NBPTS-1)*4+2)=(YB*SIGC-YC*SIGB)/(SIGC-SIGB)
  371. XCOOR((NBPTS-1)*4+3)=(ZB*SIGC-ZC*SIGB)/(SIGC-SIGB)
  372. XCOOR(NBPTS*4)=0.D0
  373. IF (VCHCA.NE.0)
  374. $ VCHCA(**)=(VCHCA(IB)*SIGC-VCHCA(IC)*SIGB)/(SIGC
  375. $ -SIGB)
  376. ISUP=NBPTS
  377. ISEGM(**)=IB1
  378. ISEGM(**)=IC1
  379. ISEGM(**)=ISUP
  380. if (mcham.ne.0)
  381. $ xsegm(**)=(
  382. $ melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigc
  383. $ -melva1.velche(min(lva1,icfa),min(idi2,jel))*sigb)
  384. $ /(sigc-sigb)
  385. 107 CONTINUE
  386. IPOS=IPOS+1
  387. JTT(IPOS)=ISUP
  388. IF (SIGC.LE.0.) THEN
  389. ISEG(IPOS)=ISEG2
  390. ELSE
  391. ISEG(IPOS)=1
  392. ENDIF
  393. ENDIF
  394. IF (SIGC.LE.0.) THEN
  395. IPOS=IPOS+1
  396. JTT(IPOS)=IC
  397. ISEG(IPOS)=ISEG3
  398. if (mcham.ne.0)
  399. $ xsegm(ic)=melva1.velche(min(lva1,icfa),min(idi2,jel))
  400. ENDIF
  401. IF (SIGC*SIGA.LT.0..AND.(.NOT.(SIGC.LE.0..AND.SIGA.LE.0.)))
  402. $ THEN
  403. IC1=MIN(IC,IA)
  404. IA1=MAX(IC,IA)
  405. DO 108 NSEGM=ISEGM(/1)-2,1,-3
  406. IF (IC1.NE.ISEGM(NSEGM)) GOTO 108
  407. IF (IA1.NE.ISEGM(NSEGM+1)) GOTO 108
  408. ISUP=ISEGM(NSEGM+2)
  409. if (mcham.ne.0)
  410. $ xsegm(isup)=(
  411. $ melva1.velche(min(lva1,icfa),min(idi2,jel))*siga
  412. $ -melva1.velche(min(lva1,iafa),min(idi2,jel))*sigc)
  413. $ /(siga-sigc)
  414. * WRITE (6,*) ' POINT EN DOUBLE ',IC1,IA1,ISUP
  415. GOTO 109
  416. 108 CONTINUE
  417. NBPTS=XCOOR(/1)/4+1
  418. SEGADJ MCOORD
  419. XCOOR((NBPTS-1)*4+1)=(XC*SIGA-XA*SIGC)/(SIGA-SIGC)
  420. XCOOR((NBPTS-1)*4+2)=(YC*SIGA-YA*SIGC)/(SIGA-SIGC)
  421. XCOOR((NBPTS-1)*4+3)=(ZC*SIGA-ZA*SIGC)/(SIGA-SIGC)
  422. XCOOR(NBPTS*4)=0.D0
  423. IF (VCHCA.NE.0)
  424. $ VCHCA(**)=(VCHCA(IC)*SIGA-VCHCA(IA)*SIGC)/(SIGA
  425. $ -SIGC)
  426. ISUP=NBPTS
  427. ISEGM(**)=IC1
  428. ISEGM(**)=IA1
  429. ISEGM(**)=ISUP
  430. if (mcham.ne.0)
  431. $ xsegm(**)=(
  432. $ melva1.velche(min(lva1,icfa),min(idi2,jel))*siga
  433. $ -melva1.velche(min(lva1,iafa),min(idi2,jel))*sigc)
  434. $ /(siga-sigc)
  435. 109 CONTINUE
  436. IPOS=IPOS+1
  437. JTT(IPOS)=ISUP
  438. IF (SIGA.LE.0.) THEN
  439. ISEG(IPOS)=ISEG3
  440. ELSE
  441. ISEG(IPOS)=1
  442. ENDIF
  443. ENDIF
  444. IF (IPOS.GT.4) WRITE (6,*) ' ANOMALIE IPOS SURFACE ',IPOS
  445. IF (IPOS.LT.3) WRITE (6,*) ' ANOMALIE IPOS SURFACE ',IPOS
  446. IF (IPOS.LT.3) GOTO 101
  447. IF (IPOS.EQ.3) THEN
  448. IF (VOLUM) THEN
  449. IA1=MIN(JTT(1),JTT(2),JTT(3))
  450. IB1=MAX(JTT(1),JTT(2),JTT(3))
  451. IC1=JTT(1)+JTT(2)+JTT(3)
  452. DO 110 IDJF=ITRT(/1)-4,1,-5
  453. ITR1=MIN(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
  454. ITR2=MAX(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
  455. ITR3=ITRT(IDJF)+ITRT(IDJF+1)+ITRT(IDJF+2)
  456. IF (IA1.NE.ITR1) GOTO 110
  457. IF (IB1.NE.ITR2) GOTO 110
  458. IF (IC1.NE.ITR3) GOTO 110
  459. * WRITE (6,*) ' TRIANGLE EN DOUBLE CAS 2 ',JTT(1),JTT(2),JTT(3)
  460. C LE TRIANGLE EST DECLARE NON TRACABLE
  461. ITRT(IDJF+4)=MOD(ITRT(IDJF+4),16)+16
  462. GOTO 303
  463. 110 CONTINUE
  464. ENDIF
  465. * WRITE (6,*) ' NOUVEAU TRIANGLE ',JTT(1),JTT(2),JTT(3)
  466. if (isect.eq.1.and.jtt(1).le.nbpoin) jtt(1)=isup
  467. ITRT(**)=JTT(1)
  468. if (isect.eq.1.and.jtt(2).le.nbpoin) jtt(2)=isup
  469. ITRT(**)=JTT(2)
  470. if (isect.eq.1.and.jtt(3).le.nbpoin) jtt(3)=isup
  471. ITRT(**)=JTT(3)
  472. ITRT(**)=IPT1.ICOLOR(JEL)
  473. ITRT(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)
  474. if (mcham.ne.0) then
  475. xtrt(**)=xsegm(jtt(1))
  476. xtrt(**)=xsegm(jtt(2))
  477. xtrt(**)=xsegm(jtt(3))
  478. endif
  479. ENDIF
  480. 303 CONTINUE
  481. IF (IPOS.EQ.4) THEN
  482. IF (VOLUM) THEN
  483. IA1=MIN(JTT(1),JTT(2),JTT(3),JTT(4))
  484. IB1=MAX(JTT(1),JTT(2),JTT(3),JTT(4))
  485. IC1=JTT(1)+JTT(2)+JTT(3)+JTT(4)
  486. ID1=JTT(1)*JTT(2)*JTT(3)*JTT(4)
  487. DO 111 IDJF=IQUAT(/1)-5,1,-6
  488. ITR1=MIN(IQUAT(IDJF),IQUAT(IDJF+1),IQUAT(IDJF+2),IQUAT(IDJF+3))
  489. ITR2=MAX(IQUAT(IDJF),IQUAT(IDJF+1),IQUAT(IDJF+2),IQUAT(IDJF+3))
  490. ITR3=IQUAT(IDJF)+IQUAT(IDJF+1)+IQUAT(IDJF+2)+IQUAT(IDJF+3)
  491. ITR4=IQUAT(IDJF)*IQUAT(IDJF+1)*IQUAT(IDJF+2)*IQUAT(IDJF+3)
  492. IF (IA1.NE.ITR1) GOTO 111
  493. IF (IB1.NE.ITR2) GOTO 111
  494. IF (IC1.NE.ITR3) GOTO 111
  495. IF (ID1.NE.ITR4) GOTO 111
  496. * WRITE (6,*) ' QUADRILA EN DOUBLE ',JTT(1),JTT(2),JTT(3),JTT(4)
  497. C LE QUADRILATERE EST DECLARE NON TRACABLE
  498. IQUAT(IDJF+5)=MOD(IQUAT(IDJF+5),16)+16
  499. GOTO 302
  500. 111 CONTINUE
  501. ENDIF
  502. * WRITE (6,*) ' NOUVEAU QUADRILATERE ',JTT(1),JTT(2),JTT(3),JTT(4)
  503. if (isect.eq.1.and.jtt(1).le.nbpoin) jtt(1)=isup
  504. IQUAT(**)=JTT(1)
  505. if (isect.eq.1.and.jtt(2).le.nbpoin) jtt(2)=isup
  506. IQUAT(**)=JTT(2)
  507. if (isect.eq.1.and.jtt(3).le.nbpoin) jtt(3)=isup
  508. IQUAT(**)=JTT(3)
  509. if (isect.eq.1.and.jtt(4).le.nbpoin) jtt(4)=isup
  510. IQUAT(**)=JTT(4)
  511. IQUAT(**)=IPT1.ICOLOR(JEL)
  512. IQUAT(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)+8*ISEG(4)
  513. if (mcham.ne.0) then
  514. xquat(**)=xsegm(jtt(1))
  515. xquat(**)=xsegm(jtt(2))
  516. xquat(**)=xsegm(jtt(3))
  517. xquat(**)=xsegm(jtt(4))
  518. endif
  519. ENDIF
  520. 302 CONTINUE
  521. GOTO 101
  522. 101 CONTINUE
  523. 350 CONTINUE
  524. IF (.NOT.VOLUM) GOTO 160
  525. C ETUDE POUR LE CAS DES VOLUMES DU TETRAEDRE S'APPUYANT SUR LE
  526. C TRIANGLE ET LE CENTRE DE GRAVITE
  527. C TETRAEDRE
  528. IF (SIGA.LT.0..AND.SIGB.LT.0..AND.SIGC.LT.0..AND.SIGRAV
  529. # .LT.0.) GOTO 160
  530. IF (SIGA.GT.0..AND.SIGB.GT.0..AND.SIGC.GT.0..AND.SIGRAV
  531. # .GT.0.) GOTO 160
  532. IPOS=0
  533. ISEG(1)=0
  534. ISEG(2)=0
  535. ISEG(3)=0
  536. ISEG(4)=0
  537. IF (SIGRAV*SIGA.LE.0.) THEN
  538. IF (SIGA.EQ.0.) THEN
  539. ISUP=IA
  540. if (mcham.ne.0)
  541. > xsegm(isup)=melva1.velche(min(lva1,iafa),min(idi2,jel))
  542. GOTO 359
  543. ENDIF
  544. DO 358 NSEGM=1,NJGRAV
  545. IF (IA.NE.IJGRAV(1,NSEGM)) GOTO 358
  546. ISUP=IJGRAV(2,NSEGM)
  547. if (mcham.ne.0)
  548. > xsegm(isup)=(melva1.velche(min(lva1,iafa),min(idi2,jel))*sigrav-
  549. > chgrav*siga)/(sigrav-siga)
  550. GOTO 359
  551. 358 CONTINUE
  552. NBPTS=XCOOR(/1)/4+1
  553. SEGADJ MCOORD
  554. XCOOR((NBPTS-1)*4+1)=(XA*SIGRAV-XGRAV*SIGA)/(SIGRAV-SIGA)
  555. XCOOR((NBPTS-1)*4+2)=(YA*SIGRAV-YGRAV*SIGA)/(SIGRAV-SIGA)
  556. XCOOR((NBPTS-1)*4+3)=(ZA*SIGRAV-ZGRAV*SIGA)/(SIGRAV-SIGA)
  557. XCOOR(NBPTS*4)=0.D0
  558. IF (VCHCA.NE.0)
  559. # VCHCA(**)=(VCHCA(IA)*SIGRAV-VGRAV*SIGA)/(SIGRAV-SIGA)
  560. if (mcham.ne.0)
  561. > xsegm(**)=(melva1.velche(min(lva1,iafa),min(idi2,jel))*sigrav-
  562. > chgrav*siga)/(sigrav-siga)
  563. ISUP=NBPTS
  564. NJGRAV=NJGRAV+1
  565. IJGRAV(1,NJGRAV)=IA
  566. IJGRAV(2,NJGRAV)=ISUP
  567. 359 CONTINUE
  568. IPOS=IPOS+1
  569. JTT(IPOS)=ISUP
  570. ENDIF
  571. IF (SIGA*SIGB.LT.0..AND.
  572. # (.NOT.(SIGRAV*SIGA.LE.0..AND.SIGRAV*SIGB.LE.0.))) THEN
  573. IA1=MIN(IA,IB)
  574. IB1=MAX(IA,IB)
  575. DO 404 NSEGM=ISEGM(/1)-2,1,-3
  576. IF (IA1.NE.ISEGM(NSEGM)) GOTO 404
  577. IF (IB1.NE.ISEGM(NSEGM+1)) GOTO 404
  578. ISUP=ISEGM(NSEGM+2)
  579. * WRITE (6,*) ' POINT EN DOUBLE ',IA1,IB1,ISUP
  580. if (mcham.ne.0)
  581. > xsegm(isup)=(melva1.velche(min(lva1,iafa),min(idi2,jel))*sigb-
  582. > melva1.velche(min(lva1,ibfa),min(idi2,jel))*siga)/(sigb-siga)
  583. GOTO 405
  584. 404 CONTINUE
  585. NBPTS=XCOOR(/1)/4+1
  586. SEGADJ MCOORD
  587. XCOOR((NBPTS-1)*4+1)=(XA*SIGB-XB*SIGA)/(SIGB-SIGA)
  588. XCOOR((NBPTS-1)*4+2)=(YA*SIGB-YB*SIGA)/(SIGB-SIGA)
  589. XCOOR((NBPTS-1)*4+3)=(ZA*SIGB-ZB*SIGA)/(SIGB-SIGA)
  590. XCOOR(NBPTS*4)=0.D0
  591. IF (VCHCA.NE.0)
  592. # VCHCA(**)=(VCHCA(IA)*SIGB-VCHCA(IB)*SIGA)/(SIGB-SIGA)
  593. ISUP=NBPTS
  594. ISEGM(**)=IA1
  595. ISEGM(**)=IB1
  596. ISEGM(**)=ISUP
  597. if (mcham.ne.0)
  598. > xsegm(**)=(melva1.velche(min(lva1,iafa),jel)*sigb-
  599. > melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga)
  600. 405 CONTINUE
  601. IPOS=IPOS+1
  602. JTT(IPOS)=ISUP
  603. IF (SIGRAV*SIGB.GT.0.) THEN
  604. ISEG(IPOS)=1
  605. ENDIF
  606. ENDIF
  607. IF (SIGRAV*SIGB.LE.0.) THEN
  608. IF (SIGB.EQ.0.) THEN
  609. ISUP=IB
  610. if (mcham.ne.0)
  611. > xsegm(isup)=melva1.velche(min(lva1,ibfa),min(idi2,jel))
  612. GOTO 369
  613. ENDIF
  614. DO 368 NSEGM=1,NJGRAV
  615. IF (IB.NE.IJGRAV(1,NSEGM)) GOTO 368
  616. ISUP=IJGRAV(2,NSEGM)
  617. GOTO 369
  618. 368 CONTINUE
  619. NBPTS=XCOOR(/1)/4+1
  620. SEGADJ MCOORD
  621. XCOOR((NBPTS-1)*4+1)=(XB*SIGRAV-XGRAV*SIGB)/(SIGRAV-SIGB)
  622. XCOOR((NBPTS-1)*4+2)=(YB*SIGRAV-YGRAV*SIGB)/(SIGRAV-SIGB)
  623. XCOOR((NBPTS-1)*4+3)=(ZB*SIGRAV-ZGRAV*SIGB)/(SIGRAV-SIGB)
  624. XCOOR(NBPTS*4)=0.D0
  625. IF (VCHCA.NE.0)
  626. # VCHCA(**)=(VCHCA(IB)*SIGRAV-VGRAV*SIGB)/(SIGRAV-SIGB)
  627. if (mcham.ne.0)
  628. > xsegm(**)=(melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigrav-
  629. > chgrav*sigb)/(sigrav-sigb)
  630. ISUP=NBPTS
  631. NJGRAV=NJGRAV+1
  632. IJGRAV(1,NJGRAV)=IB
  633. IJGRAV(2,NJGRAV)=ISUP
  634. 369 CONTINUE
  635. IPOS=IPOS+1
  636. JTT(IPOS)=ISUP
  637. if (mcham.ne.0)
  638. > xsegm(isup)=(melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigrav-
  639. > chgrav*sigb)/(sigrav-sigb)
  640. ENDIF
  641. IF (SIGB*SIGC.LT.0..AND.
  642. # (.NOT.(SIGRAV*SIGB.LE.0..AND.SIGRAV*SIGC.LE.0.))) THEN
  643. IB1=MIN(IB,IC)
  644. IC1=MAX(IB,IC)
  645. DO 406 NSEGM=ISEGM(/1)-2,1,-3
  646. IF (IB1.NE.ISEGM(NSEGM)) GOTO 406
  647. IF (IC1.NE.ISEGM(NSEGM+1)) GOTO 406
  648. ISUP=ISEGM(NSEGM+2)
  649. * WRITE (6,*) ' POINT EN DOUBLE ',IB1,IC1,ISUP
  650. if (mcham.ne.0)
  651. > xsegm(isup)=(melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigc-
  652. > melva1.velche(min(lva1,icfa),min(idi2,jel))*sigb)/(sigc-sigb)
  653. GOTO 407
  654. 406 CONTINUE
  655. NBPTS=XCOOR(/1)/4+1
  656. SEGADJ MCOORD
  657. XCOOR((NBPTS-1)*4+1)=(XB*SIGC-XC*SIGB)/(SIGC-SIGB)
  658. XCOOR((NBPTS-1)*4+2)=(YB*SIGC-YC*SIGB)/(SIGC-SIGB)
  659. XCOOR((NBPTS-1)*4+3)=(ZB*SIGC-ZC*SIGB)/(SIGC-SIGB)
  660. XCOOR(NBPTS*4)=0.D0
  661. IF (VCHCA.NE.0)
  662. # VCHCA(**)=(VCHCA(IB)*SIGC-VCHCA(IC)*SIGB)/(SIGC-SIGB)
  663. ISUP=NBPTS
  664. ISEGM(**)=IB1
  665. ISEGM(**)=IC1
  666. ISEGM(**)=ISUP
  667. if (mcham.ne.0)
  668. > xsegm(**)=(melva1.velche(min(lva1,ibfa),min(idi2,jel))*sigc-
  669. > melva1.velche(min(lva1,icfa),min(idi2,jel))*sigb)/(sigc-sigb)
  670. 407 CONTINUE
  671. IPOS=IPOS+1
  672. JTT(IPOS)=ISUP
  673. IF (SIGRAV*SIGC.GT.0.) THEN
  674. ISEG(IPOS)=1
  675. ENDIF
  676. ENDIF
  677. IF (SIGRAV*SIGC.LE.0.) THEN
  678. IF (SIGC.EQ.0.) THEN
  679. ISUP=IC
  680. if (mcham.ne.0)
  681. > xsegm(isup)=melva1.velche(min(lva1,icfa),min(idi2,jel))
  682. GOTO 379
  683. ENDIF
  684. DO 378 NSEGM=1,NJGRAV
  685. IF (IC.NE.IJGRAV(1,NSEGM)) GOTO 378
  686. ISUP=IJGRAV(2,NSEGM)
  687. if (mcham.ne.0)
  688. > xsegm(isup)=(melva1.velche(min(lva1,icfa),min(idi2,jel))*sigrav-
  689. > chgrav*sigc)/(sigrav-sigc)
  690. GOTO 379
  691. 378 CONTINUE
  692. NBPTS=XCOOR(/1)/4+1
  693. SEGADJ MCOORD
  694. XCOOR((NBPTS-1)*4+1)=(XC*SIGRAV-XGRAV*SIGC)/(SIGRAV-SIGC)
  695. XCOOR((NBPTS-1)*4+2)=(YC*SIGRAV-YGRAV*SIGC)/(SIGRAV-SIGC)
  696. XCOOR((NBPTS-1)*4+3)=(ZC*SIGRAV-ZGRAV*SIGC)/(SIGRAV-SIGC)
  697. XCOOR(NBPTS*4)=0.D0
  698. IF (VCHCA.NE.0)
  699. # VCHCA(**)=(VCHCA(IC)*SIGRAV-VGRAV*SIGC)/(SIGRAV-SIGC)
  700. if (mcham.ne.0)
  701. > xsegm(**)=(melva1.velche(min(lva1,icfa),min(idi2,jel))*sigrav-
  702. > chgrav*sigc)/(sigrav-sigc)
  703. ISUP=NBPTS
  704. NJGRAV=NJGRAV+1
  705. IJGRAV(1,NJGRAV)=IC
  706. IJGRAV(2,NJGRAV)=ISUP
  707. 379 CONTINUE
  708. IPOS=IPOS+1
  709. JTT(IPOS)=ISUP
  710. ENDIF
  711. IF (SIGC*SIGA.LT.0..AND.
  712. # (.NOT.(SIGRAV*SIGC.LE.0..AND.SIGRAV*SIGA.LE.0.))) THEN
  713. IC1=MIN(IC,IA)
  714. IA1=MAX(IC,IA)
  715. DO 408 NSEGM=ISEGM(/1)-2,1,-3
  716. IF (IC1.NE.ISEGM(NSEGM)) GOTO 408
  717. IF (IA1.NE.ISEGM(NSEGM+1)) GOTO 408
  718. ISUP=ISEGM(NSEGM+2)
  719. * WRITE (6,*) ' POINT EN DOUBLE ',IC1,IA1,ISUP
  720. if (mcham.ne.0)
  721. > xsegm(isup)=(melva1.velche(min(lva1,icfa),min(idi2,jel))*siga-
  722. > melva1.velche(min(lva1,iafa),min(idi2,jel))*sigc)/(siga-sigc)
  723. GOTO 409
  724. 408 CONTINUE
  725. NBPTS=XCOOR(/1)/4+1
  726. SEGADJ MCOORD
  727. XCOOR((NBPTS-1)*4+1)=(XC*SIGA-XA*SIGC)/(SIGA-SIGC)
  728. XCOOR((NBPTS-1)*4+2)=(YC*SIGA-YA*SIGC)/(SIGA-SIGC)
  729. XCOOR((NBPTS-1)*4+3)=(ZC*SIGA-ZA*SIGC)/(SIGA-SIGC)
  730. XCOOR(NBPTS*4)=0.D0
  731. IF (VCHCA.NE.0)
  732. # VCHCA(**)=(VCHCA(IC)*SIGA-VCHCA(IA)*SIGC)/(SIGA-SIGC)
  733. ISUP=NBPTS
  734. ISEGM(**)=IC1
  735. ISEGM(**)=IA1
  736. ISEGM(**)=ISUP
  737. if (mcham.ne.0)
  738. > xsegm(**)=(melva1.velche(min(lva1,icfa),min(idi2,jel))*siga-
  739. > melva1.velche(min(lva1,iafa),min(idi2,jel))*sigc)/(siga-sigc)
  740. 409 CONTINUE
  741. IPOS=IPOS+1
  742. JTT(IPOS)=ISUP
  743. IF (SIGRAV*SIGA.GT.0.) THEN
  744. ISEG(IPOS)=1
  745. ENDIF
  746. ENDIF
  747. IF (IPOS.GT.4) WRITE (6,*) ' ANOMALIE IPOS VOLUME ',IPOS
  748. IF (IPOS.LT.3) WRITE (6,*) ' ANOMALIE IPOS VOLUME ',IPOS
  749. IF (IPOS.LT.3) GOTO 160
  750. ITR(**)=JTT(1)
  751. ITR(**)=JTT(2)
  752. ITR(**)=JTT(3)
  753. ITR(**)=IPT1.ICOLOR(JEL)
  754. if (mcham.ne.0) then
  755. xtr(**)=xsegm(jtt(1))
  756. xtr(**)=xsegm(jtt(2))
  757. xtr(**)=xsegm(jtt(3))
  758. endif
  759. * WRITE (6,*) ' TRIANGLE (CAS 1) DU PLAN DE COUPE '
  760. * WRITE (6,*) JTT(1),JTT(2),JTT(3)
  761. C TRIANGLE DU PLAN DE COUPE
  762. IF (IPOS.EQ.3) THEN
  763. MCOUP(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)+8
  764. ELSE
  765. MCOUP(**)=ISEG(1)+2*ISEG(2)+4*0 +8
  766. ENDIF
  767. 403 CONTINUE
  768. IF (IPOS.GE.4) THEN
  769. ITR(**)=JTT(1)
  770. ITR(**)=JTT(3)
  771. ITR(**)=JTT(4)
  772. ITR(**)=IPT1.ICOLOR(JEL)
  773. if (mcham.ne.0) then
  774. xtr(**)=xsegm(jtt(1))
  775. xtr(**)=xsegm(jtt(3))
  776. xtr(**)=xsegm(jtt(4))
  777. endif
  778. * WRITE (6,*) ' TRIANGLE (CAS 2) DU PLAN DE COUPE '
  779. * WRITE (6,*) JTT(1),JTT(3),JTT(4)
  780. C TRIANGLE DU PLAN DE COUPE
  781. MCOUP(**)=0 +2*ISEG(3)+4*ISEG(4)+8
  782. ENDIF
  783. 402 CONTINUE
  784. 160 CONTINUE
  785. 161 CONTINUE
  786. 20 CONTINUE
  787. NBELEM=JELN
  788. SEGADJ IPT2
  789. ISU=ISU+1
  790. IPT8.LISOUS(ISU)=IPT2
  791. if (mcham.ne.0) then
  792. n1el=nbelem
  793. segadj melva2
  794. ipt8.lisref(isu)=melva2
  795. endif
  796. SEGDES IPT1
  797. 10 CONTINUE
  798. IF (NBSOU.NE.0) SEGDES MELEME
  799. NBSOUS=ISU+1
  800. * IF (ITR(/1).EQ.0) NBSOUS=ISU
  801. NBELEM=0
  802. NBNN=0
  803. NBREF=0
  804. if (mcham.ne.0) nbref=nbsous
  805. SEGADJ IPT8
  806. MELEME=IPT8
  807. * IF (ITR(/1).EQ.0) RETURN
  808. * WRITE(6,*) ' NOMBRE DE TRIANGLE ENGENDREES DANS LE PLAN DE COUPE',
  809. * $ ITR(/1)/4
  810. * WRITE(6,*) ' NOMBRE DE TRIANGLE ENGENDREES HORS LE PLAN DE COUPE',
  811. * $ ITRT(/1)/5
  812. * WRITE(6,*) ' NOMBRE DE QUADRILA ENGENDREES HORS LE PLAN DE COUPE',
  813. * $ IQUAT(/1)/6
  814. NBELEM=ITR(/1)/4+ITRT(/1)/5+(IQUAT(/1)/6)*2
  815. NBNN=3
  816. NBREF=0
  817. NBSOUS=0
  818. SEGINI IPT4
  819. IPT4.ITYPEL=4
  820. if (mcham.ne.0) then
  821. n1el=nbelem
  822. n1ptel=nbnn
  823. segini melva4
  824. endif
  825. ITRL1=ITR(/1)/4
  826. DO 200 IEL=1,ITRL1
  827. DO 201 NBNN=1,3
  828. IPT4.NUM(NBNN,IEL)=ITR((IEL-1)*4+NBNN)
  829. if (mcham.ne.0) melva4.velche(nbnn,iel)=xtr((IEL-1)*3+nbnn)
  830. 201 CONTINUE
  831. IPT4.ICOLOR(IEL)=ITR(IEL*4)
  832. 200 CONTINUE
  833. ITRL2=ITRT(/1)/5
  834. DO 202 IEL=1,ITRL2
  835. DO 203 NBNN=1,3
  836. IPT4.NUM(NBNN,IEL+ITRL1)=ITRT((IEL-1)*5+NBNN)
  837. if (mcham.ne.0)
  838. > melva4.velche(nbnn,iel+itrl1)=xtrt((IEL-1)*3+nbnn)
  839. 203 CONTINUE
  840. IPT4.ICOLOR(IEL+ITRL1)=ITRT(IEL*5-1)
  841. MCOUP(**)=ITRT(IEL*5)
  842. 202 CONTINUE
  843. ITRL3=(IQUAT(/1)/6)
  844. DO 204 IEL=1,ITRL3
  845. IEL1=2*IEL-1
  846. DO 205 NBNN=1,3
  847. IPT4.NUM(NBNN,IEL1+ITRL1+ITRL2)=IQUAT((IEL-1)*6+NBNN)
  848. if (mcham.ne.0)
  849. > melva4.velche(nbnn,iel1+itrl1+itrl2)=xquat((IEL-1)*4+nbnn)
  850. 205 CONTINUE
  851. IPT4.ICOLOR(IEL1+ITRL1+ITRL2)=IQUAT(IEL*6-1)
  852. MCOUP(**)=MOD(IQUAT(IEL*6),4)+16*(IQUAT(IEL*6)/16)
  853. IEL2=2*IEL
  854. IPT4.NUM(1,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+1)
  855. IPT4.NUM(2,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+3)
  856. IPT4.NUM(3,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+4)
  857. IPT4.ICOLOR(IEL2+ITRL1+ITRL2)=IQUAT(IEL*6-1)
  858. MCOUP(**)=0+MOD(IQUAT(IEL*6)/4,2)*2+MOD(IQUAT(IEL*6)/8,2)*4
  859. # +16*(IQUAT(IEL*6)/16)
  860. if (mcham.ne.0) then
  861. melva4.velche(1,iel2+itrl1+itrl2)=xquat((IEL-1)*4+1)
  862. melva4.velche(2,iel2+itrl1+itrl2)=xquat((IEL-1)*4+3)
  863. melva4.velche(3,iel2+itrl1+itrl2)=xquat((IEL-1)*4+4)
  864. endif
  865. 204 CONTINUE
  866. * WRITE (6,*) ' LISTE DES ELEMENTS RAJOUTES '
  867. * DO 9834 J=1,NBELEM
  868. * WRITE (6,*) (IPT4.NUM(I,J),I=1,3)
  869. *9834 CONTINUE
  870. LISOUS(ISU+1)=IPT4
  871. if (mcham.ne.0) lisref(ISU+1)=melva4
  872. IF (MELEM2.NE.0) THEN
  873. IF (MEDSAU.EQ.0) THEN
  874. C ON COUPE MELEM2
  875. mcham=0
  876. MEDSAU=MELEME
  877. MEDCOU=MCOUP
  878. MELEME=MELEM2
  879. SEGSUP ITR,ITRT,IQUAT
  880. SEGINI MCOUP,ITR,ITRT,IQUAT
  881. GOTO 5000
  882. ELSE
  883. MELEM2=MELEME
  884. MELEME=MEDSAU
  885. MCOU2=MCOUP
  886. MCOUP=MEDCOU
  887. ENDIF
  888. ENDIF
  889. SEGSUP ITR,ISEGM,IJGRAV,ITRT,IQUAT
  890. if (mchamo.ne.0) segsup xtr,xtrt,xquat,xsegm
  891. RETURN
  892. END
  893.  
  894.  
  895.  

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