Télécharger poiext.eso

Retour à la liste

Numérotation des lignes :

  1. C POIEXT SOURCE BP208322 16/11/30 21:15:11 9223
  2.  
  3. c=======================================================================
  4. C Sert a extraire un point d'un objet MAILLAGE
  5. C
  6. C Novembre 1985 :
  7. C Extrait d'un objet tous les points appartenant a une DROIte ou a
  8. C un PLAN (possibilite de critere de distance)
  9. C
  10. C Fevrier 1986 :
  11. C Extrait d'un objet tous les points appartenant a un CYLIndre, un
  12. C CONE, une SPHEre ou un TORE
  13. C
  14. C 09/2003 :
  15. C Introduction de la DIMEnsion 1. Ajout de la syntaxe pour creer un
  16. C point a partir de ses coordonnees. Inhibation d'options non compa-
  17. C tibles avec la dimension 1. Option PROC adapte au cas IDIM=1.
  18. C Ajout de tests sur le type d'element pour POINT Entier,
  19. C POINT INITIAL et POINT FINAL (conforme a la notice).
  20. C
  21. C 01/2004 :
  22. C Modifications des methodes de recherche des points appartenant a
  23. C une DROIte, un PLAN, un CYLIndre, un CONE. Methodes moins sensibles
  24. C aux imprecisions numeriques.
  25. C
  26. C 07/2013 GOUNAND
  27. C Les options INIT FINA et PROC détectent les maillages vides
  28. C en entrée et provoquent une erreur
  29. C
  30. C 04/2015 CB215821
  31. C Ajout de la possibilite de generer un MAILLAGE de POI1 en fournissant
  32. C des LISTREELS de coordonnees et optionnellement un LISTREEL de densite
  33. C
  34. C 11/2016 JCARDO
  35. C Ajout de l'option 'JONC' (en liaison avec les developpements dans
  36. C l'operateur PART)
  37. C
  38. c 11/2016 BP
  39. c ajout de la syntaxe : ptij = POIN ieme_noeud jeme_element maillage;
  40. c
  41. c=======================================================================
  42.  
  43. SUBROUTINE POIEXT
  44.  
  45. IMPLICIT INTEGER(I-N)
  46. IMPLICIT REAL*8 (A-H,O-Z)
  47.  
  48. -INC CCOPTIO
  49. -INC SMELEME
  50. -INC SMLREEL
  51. -INC SMCOORD
  52. -INC CCGEOME
  53.  
  54. SEGMENT SLREEL
  55. INTEGER LISLRE(NBLREE)
  56. ENDSEGMENT
  57.  
  58. C- DIMENSION XNORM(3)
  59. CHARACTER*4 MOT,IPRO
  60. CHARACTER*4 MCLE(10),MOTM(9),MOABS(1),MOTAV(2)
  61. SEGMENT INDIC(NBELEM)
  62.  
  63. SEGMENT NETN(XCOOR(/1)/(IDIM+1)+1)
  64.  
  65. INTEGER IERREU
  66. INTEGER ICOUR
  67.  
  68. C COMMON de sauvegarde lors d'un appel par SYMTRI
  69. COMMON / CSYMTR / XNORM(3)
  70.  
  71. DATA MCLE / 'INIT','FINA','JONC','PROC',
  72. . 'DROI','PLAN','CYLI','CONE','SPHE','TORE' /
  73. DATA MOTM / 'MAXI','MINI','SUPE','EGSU','EGAL','EGIN','INFE',
  74. . 'DIFF','COMP' /
  75. DATA MOTAV / 'AVEC','SANS' /
  76. DATA MOABS / 'ABS ' /
  77.  
  78.  
  79. ICOUR = 0
  80. IERREU= 0
  81.  
  82. c Lecture des mots cles correspondant a la 3eme syntaxe (CHPOINT)
  83. CALL LIRMOT(MOTM,9,IMM,0)
  84. IF (IMM.NE.0) GOTO 300
  85. c syntaxe 3 ==> GOTO 300
  86.  
  87.  
  88. idimp1=IDIM+1
  89. c Lecture des mots cles geometrique
  90. CALL LIRMOT(MCLE(5),6,IMC,0)
  91. IF (IMC.NE.0) GOTO 250
  92. c syntaxe 2.b ==> GOTO 250
  93.  
  94.  
  95. c Lecture d'un maillage en entree
  96. CALL LIROBJ('MAILLAGE',MELEME,0,IRetou)
  97. IF (IERR.NE.0) RETURN
  98. IF (IRetou.NE.0) GOTO 200
  99. c syntaxe 2.a ==> GOTO 200
  100.  
  101.  
  102. 100 CONTINUE
  103. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) 'POINT syntaxe 1:'
  104. C -------------------------------
  105. C Operateur POINT - syntaxe 1 :
  106. C -------------------------------
  107. C CB215821 Ajout de la possibilité de générer un MAILLAGE de POI1 en
  108. C fournissant directement un LISTREEL par COORDONNEES plus
  109. C optionnellement un LISTREEL pour la densité
  110.  
  111. CALL LIROBJ('LISTREEL',MLREEL,0,IRETOU)
  112. IF (IRETOU.NE.0) THEN
  113. IF(IDIM .EQ. 0) THEN
  114. C Dimension non définie
  115. CALL ERREUR(14)
  116. RETURN
  117. ENDIF
  118. NBLREE = IDIM+1
  119. SEGINI,SLREEL
  120.  
  121. LISLRE(1)=MLREEL
  122. SEGACT,MLREEL
  123. NBVAL = PROG(/1)
  124. JG = NBVAL
  125. C Lecture OBLIGATOIRE d'un nombre de LISTREEL egal a IDIM
  126. C Lecture OPTIONNELLE de la densité comme dernier LISTREEL
  127. DO I=2,IDIM+1
  128. IF (I .LE. IDIM) THEN
  129. CALL LIROBJ('LISTREEL',MLREEL,1,IRETOU)
  130. IF (IRETOU .EQ. 0) RETURN
  131. SEGACT,MLREEL
  132. ELSE
  133. C Lecture de la densite OPTIONNELLE
  134. CALL LIROBJ('LISTREEL',MLREEL,0,IRETOU)
  135. IF (IRETOU .EQ. 0) THEN
  136. SEGINI,MLREEL
  137. C Remplissage du LISTREEL avec la densite globale
  138. DO J=1,NBVAL
  139. PROG(J)=DENSIT
  140. ENDDO
  141. ELSE
  142. SEGACT,MLREEL
  143. ENDIF
  144. ENDIF
  145.  
  146. LISLRE(I)=MLREEL
  147.  
  148. C Verification de la longeur des LISTREELS lus
  149. IF (PROG(/1) .NE. NBVAL) THEN
  150. ICOUR = I
  151. IERREU = 1
  152. GOTO 101
  153. ENDIF
  154. ENDDO
  155.  
  156. C On remplit le MELEME de POI1
  157. IF (NBVAL .EQ. 0) THEN
  158. C Gestion des MAILLAGES VIDES
  159. NBNN = NBNNE(ILCOUR)
  160. NBELEM = 0
  161. NBSOUS = 0
  162. NBREF = 0
  163. SEGINI,MELEME
  164. ITYPEL = ILCOUR
  165. ELSE
  166. C Ajout des NOEUDS au SEGMENT MCOORD et ajustement
  167. SEGACT,MCOORD*MOD
  168. NBANC=XCOOR(/1)/(IDIM+1)
  169.  
  170. NBPTS=NBANC+NBVAL
  171. SEGADJ,MCOORD
  172.  
  173. NBNN = 1
  174. NBELEM = NBVAL
  175. NBSOUS = 0
  176. NBREF = 0
  177. SEGINI,MELEME
  178. ITYPEL = 1
  179.  
  180. NBAD = 0
  181. IDINI = NBANC*(IDIM+1)
  182. DO I=1,NBVAL
  183. DO J=1,IDIM+1
  184. C Le cas (IDIM+1) sert à remplir la densité
  185. NBAD = NBAD+1
  186. MLREEL=LISLRE(J)
  187. XCOOR(IDINI+NBAD)=PROG(I)
  188. ENDDO
  189. NUM(1,I) =NBANC+I
  190. ICOLOR(I)=IDCOUL
  191. ENDDO
  192. ENDIF
  193.  
  194. CALL ECROBJ('MAILLAGE',MELEME)
  195. SEGDES,MELEME
  196. ICOUR = IDIM
  197.  
  198. C Desactivation / Suppression des Segments de travail
  199. 101 CONTINUE
  200. DO J=1,ICOUR
  201. MLREEL=LISLRE(J)
  202. SEGDES,MLREEL
  203. ENDDO
  204. SEGSUP SLREEL
  205. IF (IERREU .EQ. 1) CALL ERREUR(217)
  206. RETURN
  207.  
  208.  
  209. ELSE
  210. C Creation d'un seul POINT par trois flottants
  211. C Pour IDIM = 1, seul moyen de creer un point.
  212. CALL MESLIR(-149)
  213. CALL LIRREE(Val1,1,IRetou)
  214. CALL MESLIR(-150)
  215. IF (IDIM.EQ.2) THEN
  216. CALL LIRREE(Val2,1,IRetou)
  217. ELSE IF (IDIM.EQ.3) THEN
  218. CALL LIRREE(Val2,1,IRetou)
  219. IF (IERR.NE.0) RETURN
  220. CALL MESLIR(-151)
  221. CALL LIRREE(Val3,1,IRetou)
  222. ENDIF
  223. IF (IERR.NE.0) RETURN
  224. SEGACT MCOORD*MOD
  225. NbPts=1+(XCOOR(/1)/idimp1)
  226. SEGADJ MCOORD
  227. IRef=(NbPts-1)*idimp1
  228. XCOOR(IRef+1)=Val1
  229. IF (IDIM.EQ.2) THEN
  230. XCOOR(IRef+2)=Val2
  231. ELSE IF (IDIM.EQ.3) THEN
  232. XCOOR(IRef+2)=Val2
  233. XCOOR(IRef+3)=Val3
  234. ENDIF
  235. ENDIF
  236. XCOOR(IRef+idimp1)=DENSIT
  237. CALL ECROBJ('POINT ',NbPts)
  238. RETURN
  239.  
  240.  
  241. 200 CONTINUE
  242. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) 'POINT syntaxe 2:'
  243. C -------------------------------
  244. C Operateur POINT - syntaxe 2 :
  245. C -------------------------------
  246. C 1) Recuperation du POINT1 du maillage GEO1 par son numero
  247. C Point1 = POINT Geo1 Entier1 ;
  248. C Note : le maillage est constitue uniquement de POI1 ou SEG2.
  249. C 1') Recuperation de POINT1 = Entier1^eme noeud du Entier2^eme element de GEO1
  250. C Point1 = POINT Geo1 Entier1 Entier2;
  251. C Note : le maillage est constitue uniquement de POI1 ou SEG2.
  252. C 2) Recuperation du point INITIAL ou FINAL du maillage GEO1
  253. C Point1 = POINT Geo1 'INITIAL' | 'FINAL' ;
  254. C Note : le maillage est constitue uniquement de POI1, SEG2 ou SEG3.
  255. C 3) Recuperation du POINT1 du maillage GEO1 le plus proche du POINT2.
  256. C Point1 = POINT Geo1 'PROC' Point2 ;
  257.  
  258. C Lecture facultative d'un entier (numero du point a trouver)
  259. CALL LIRENT(JPTR,0,IRetou)
  260.  
  261. c cas ou on a bien lu un entier
  262. IF (IRetou.NE.0) THEN
  263.  
  264. IF (JPTR.LE.0) THEN
  265. CALL ERREUR(262)
  266. RETURN
  267. ENDIF
  268. MOT=' '
  269. IPRO=' '
  270. C Lecture d'un 2eme entier
  271. CALL LIRENT(KPTR,0,IRetou)
  272. IF (IRetou.NE.0) THEN
  273. IF (KPTR.LE.0) THEN
  274. CALL ERREUR(262)
  275. RETURN
  276. ENDIF
  277. MOT='MELE'
  278. ENDIF
  279.  
  280. ELSE
  281. C Lecture d'un des mots cles ('INIT','FINA','JONC','PROC')
  282.  
  283. CALL LIRCHA(MOT,0,IRetou)
  284. IF (IRetou.EQ.0) THEN
  285. C Le type de l'operande est incorrect
  286. CALL QUETYP(MOTERR(1:8),0,IRetou)
  287. IF (IRetou.NE.0) THEN
  288. CALL ERREUR (39)
  289. ELSE
  290. CALL ERREUR(533)
  291. ENDIF
  292. RETURN
  293. ENDIF
  294. C point INIT ou FINA => IPRO=' '
  295. IF ((MOT.EQ.MCLE(1)).OR.(MOT.EQ.MCLE(2))) THEN
  296. IPRO=' '
  297. C point JONC ou PROCH
  298. ELSEIF (MOT.EQ.MCLE(3)) THEN
  299. IPRO=MOT
  300. GOTO 203
  301. ELSEIF (MOT.EQ.MCLE(4)) THEN
  302. IPRO=MOT
  303. GOTO 204
  304. ELSE
  305. CALL ERREUR(7)
  306. RETURN
  307. ENDIF
  308.  
  309. ENDIF
  310.  
  311. C Traitement des syntaxes POINT Geo1 Entier1 | 'INIT' | 'FINA'
  312. C Erreur si objet maillage complexe
  313. SEGACT MELEME
  314. IF (LISOUS(/1).NE.0) THEN
  315. CALL ERREUR(25)
  316. SEGDES MELEME
  317. RETURN
  318. ENDIF
  319. IF (NUM(/1).LE.0.OR.NUM(/2).LE.0) THEN
  320. MOTERR(1:8)='MAILLAGE'
  321. C Une donnée de type %m1:8 est vide
  322. CALL ERREUR(1027)
  323. SEGDES MELEME
  324. RETURN
  325. ENDIF
  326.  
  327. c --- syntaxe POINT Geo1 Entier1 Entier2 ---
  328. IF (MOT.EQ.'MELE') THEN
  329. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option Entier1 Entier2'
  330. IPTR=JPTR
  331. JPTR=KPTR
  332. IF (IPTR.GT.NUM(/1).OR.JPTR.GT.NUM(/2)) THEN
  333. CALL ERREUR(262)
  334. SEGDES MELEME
  335. RETURN
  336. ENDIF
  337. c --- syntaxe POINT Geo1 Entier1 ---
  338. ELSEIF (MOT.EQ.' ') THEN
  339. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option Entier1'
  340. IF ((ITYPEL.NE.1).AND.(JPTR.EQ.(NUM(/2)+1))) THEN
  341. IPTR=NUM(/1)
  342. JPTR=NUM(/2)
  343. ELSE
  344. IPTR=1
  345. ENDIF
  346. IF (JPTR.GT.NUM(/2)) THEN
  347. CALL ERREUR(262)
  348. SEGDES MELEME
  349. RETURN
  350. ENDIF
  351. c --- syntaxe POINT Geo1 'INIT' ---
  352. ELSE IF (MOT.EQ.MCLE(1)) THEN
  353. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option INIT'
  354. IPTR=1
  355. JPTR=1
  356. c --- syntaxe POINT Geo1 'FINA' ---
  357. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option FINA'
  358. ELSE IF (MOT.EQ.MCLE(2)) THEN
  359. IPTR=NUM(/1)
  360. JPTR=NUM(/2)
  361. ENDIF
  362. c ---partie commune aux syntaxes Entier(s) INIT FINA ---
  363. IPoin=NUM(IPTR,JPTR)
  364. CALL ECROBJ('POINT ',IPoin)
  365. SEGDES MELEME
  366. RETURN
  367.  
  368.  
  369. C --- Traitement de l'option POIN 'JONC' ---
  370. 203 CONTINUE
  371. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option JONC'
  372. SEGACT MELEME
  373. NBSOUS=LISOUS(/1)
  374. NBS=MAX(1,NBSOUS)
  375. IPT1=MELEME
  376.  
  377. NBELEM=0
  378. SEGINI,NETN
  379. DO IO=1,NBS
  380. IF (NBSOUS.GT.0) THEN
  381. IPT1=LISOUS(IO)
  382. SEGACT,IPT1
  383. ENDIF
  384. SEGACT,IPT1
  385. DO 222 J=1,IPT1.NUM(/2)
  386. DO 222 I=1,IPT1.NUM(/1)
  387. IF (NETN(IPT1.NUM(I,J)).EQ.2) NBELEM=NBELEM+1
  388. NETN(IPT1.NUM(I,J))=NETN(IPT1.NUM(I,J))+1
  389. 222 CONTINUE
  390. IF (NBSOUS.NE.0) SEGDES,IPT1
  391. ENDDO
  392. SEGDES,MELEME
  393.  
  394. NBNN=1
  395. NBREF=0
  396. NBSOUS=0
  397. SEGINI,IPT1
  398. IPT1.ITYPEL=1
  399. K=0
  400. DO I=1,NETN(/1)
  401. IF (NETN(I).GT.2) THEN
  402. K=K+1
  403. IPT1.NUM(1,K)=I
  404. ENDIF
  405. ENDDO
  406.  
  407. SEGSUP,NETN
  408.  
  409. CALL ECROBJ('MAILLAGE',IPT1)
  410.  
  411. RETURN
  412.  
  413. C --- Traitement de l'option POIN 'PROC' ---
  414. 204 CONTINUE
  415. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option PROCH'
  416. CALL LIROBJ('POINT ',IPoin1,1,IRetou)
  417. IF (IERR.NE.0) THEN
  418. RETURN
  419. ENDIF
  420. Val1=1.E30
  421. IRef=(IPoin1-1)*idimp1
  422. XP1=XCOOR(IRef+1)
  423. YP1=XCOOR(IRef+2)
  424. ZP1=0.D0
  425. IF (IDIM.GE.3) ZP1=XCOOR(IRef+3)
  426. SEGACT MELEME
  427. IPT1=MELEME
  428. Ipoin=0
  429. DO IOB=1,MAX(1,LISOUS(/1))
  430. IF (LISOUS(/1).NE.0) THEN
  431. IPT1=LISOUS(IOB)
  432. SEGACT IPT1
  433. ENDIF
  434. DO i=1,IPT1.NUM(/2)
  435. DO j=1,IPT1.NUM(/1)
  436. IRef=(IPT1.NUM(j,i)-1)*idimp1
  437. IF (IDIM.EQ.1) THEN
  438. DVN=ABS(XP1-XCOOR(IRef+1))
  439. ELSE
  440. XVN=XP1-XCOOR(IRef+1)
  441. YVN=YP1-XCOOR(IRef+2)
  442. ZVN=0.D0
  443. IF (IDIM.GE.3) ZVN=ZP1-XCOOR(IRef+3)
  444. DVN=XVN*XVN+YVN*YVN+ZVN*ZVN
  445. ENDIF
  446. IF (DVN.LE.Val1) THEN
  447. Val1=DVN
  448. IPoin=IPT1.NUM(j,i)
  449. ENDIF
  450. ENDDO
  451. ENDDO
  452. IF (LISOUS(/1).NE.0) SEGDES,IPT1
  453. ENDDO
  454. IF (IPoin.EQ.0) THEN
  455. MOTERR(1:8)='MAILLAGE'
  456. C Une donnée de type %m1:8 est vide
  457. CALL ERREUR(1027)
  458. ELSE
  459. IF (IERR.EQ.0) CALL ECROBJ('POINT ',IPoin)
  460. ENDIF
  461. SEGDES MELEME
  462. RETURN
  463.  
  464.  
  465. 250 CONTINUE
  466. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) 'POINT syntaxe 2.b :'
  467. C -------------------------------
  468. C Operateur POINT - syntaxe 2.b :
  469. C -------------------------------
  470. C Extraction du(des) point(s) du maillage GEO1 situe(s) sur une
  471. C DROIte, un PLAN, un CYLIndre, un CONE, une SPHEre ou un TORE.
  472. C L'appartenance du point a l'objet est definie suivant un critere
  473. C de distance non nul (appartenance stricte : critere tres faible).
  474. C Geof = POINT Geo1 MOT1 Points ... (Critere)
  475.  
  476. C En dimension 1, ces criteres n'ont aucun sens.
  477. IF (IDIM.EQ.1) THEN
  478. MOTERR(1:4)=MCLE(3+IMC)
  479. INTERR(1)=IDIM
  480. CALL ERREUR(971)
  481. RETURN
  482. ENDIF
  483. C Lectures obligatoire du maillage et facultative du critere
  484. CALL LIROBJ('MAILLAGE',MELEME,1,IRetou)
  485. IF (IERR.NE.0) RETURN
  486. CRIT=DENSIT/10.D0
  487. CALL LIRREE(CRIT,0,IRetou)
  488. IF (CRIT.LE.0.D0) CALL ERREUR(17)
  489. IF (IERR.NE.0) RETURN
  490. C Transformation du maillage en 'POI1'
  491. SEGACT MELEME
  492. CALL CHANGE(MELEME,1)
  493. SEGACT MELEME
  494. NBELEM=NUM(/2)
  495. C Intersection du maillage avec une DROIte
  496. C --------------------------
  497. IF (IMC.EQ.1) THEN
  498. CALL LIROBJ('POINT ',IPoin1,1,IRetou)
  499. CALL LIROBJ('POINT ',IPoin2,1,IRetou)
  500. IF (IERR.NE.0) THEN
  501. SEGDES MELEME
  502. RETURN
  503. ENDIF
  504. C Recherche du vecteur directeur de la DROIte
  505. IRef1=(IPoin1-1)*idimp1
  506. IRef2=(IPoin2-1)*idimp1
  507. XP1=XCOOR(IRef1+1)
  508. XVN=XCOOR(IRef2+1)-XP1
  509. YP1=XCOOR(IRef1+2)
  510. YVN=XCOOR(IRef2+2)-YP1
  511. ZP1=0.D0
  512. ZVN=0.D0
  513. IF (IDIM.GE.3) THEN
  514. ZP1=XCOOR(IRef1+3)
  515. ZVN=XCOOR(IRef2+3)-ZP1
  516. ENDIF
  517. DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
  518. IF (DVN.EQ.0.D0) THEN
  519. CALL ERREUR(21)
  520. SEGDES MELEME
  521. RETURN
  522. ENDIF
  523. XVN=XVN/DVN
  524. YVN=YVN/DVN
  525. ZVN=ZVN/DVN
  526. XNORM(1)=XVN
  527. XNORM(2)=YVN
  528. XNORM(3)=ZVN
  529. C Recherche des points verifiant la condition d'appartenance (INDIC)
  530. SEGINI INDIC
  531. NBON=0
  532. DO IP=1,NBELEM
  533. IPoin=NUM(1,IP)
  534. IRef=(IPoin-1)*idimp1
  535. XPT=XCOOR(IRef+1)-XP1
  536. YPT=XCOOR(IRef+2)-YP1
  537. ZPT=0.D0
  538. IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1
  539. SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
  540. XV=SCAL*XVN-XPT
  541. YV=SCAL*YVN-YPT
  542. ZV=SCAL*ZVN-ZPT
  543. Val1=XV*XV+YV*YV+ZV*ZV
  544. IF (Val1.LE.CRIT) THEN
  545. INDIC(IP)=IPoin
  546. NBON=NBON+1
  547. ENDIF
  548. ENDDO
  549. C Intersection du maillage avec un PLAN
  550. C --------------------------
  551. ELSE IF (IMC.EQ.2) THEN
  552. CALL LIROBJ('POINT ',IPoin1,1,IRetou)
  553. CALL LIROBJ('POINT ',IPoin2,1,IRetou)
  554. CALL LIROBJ('POINT ',IPoin3,1,IRetou)
  555. IF (IERR.NE.0) THEN
  556. SEGDES MELEME
  557. RETURN
  558. ENDIF
  559. C Determination du vecteur normal au PLAN
  560. IRef1=(IPoin1-1)*idimp1
  561. IRef2=(IPoin2-1)*idimp1
  562. IRef3=(IPoin3-1)*idimp1
  563. XP1=XCOOR(IRef1+1)
  564. XP2=XCOOR(IRef2+1)
  565. XP3=XCOOR(IRef3+1)
  566. YP1=XCOOR(IRef1+2)
  567. YP2=XCOOR(IRef2+2)
  568. YP3=XCOOR(IRef3+2)
  569. ZP1=0.D0
  570. ZP2=0.D0
  571. ZP3=0.D0
  572. IF (IDIM.GE.3) THEN
  573. ZP1=XCOOR(IRef1+3)
  574. ZP2=XCOOR(IRef2+3)
  575. ZP3=XCOOR(IRef3+3)
  576. ENDIF
  577. XVN=(YP2-YP1)*(ZP3-ZP1)-(ZP2-ZP1)*(YP3-YP1)
  578. YVN=(ZP2-ZP1)*(XP3-XP1)-(XP2-XP1)*(ZP3-ZP1)
  579. ZVN=(XP2-XP1)*(YP3-YP1)-(YP2-YP1)*(XP3-XP1)
  580. DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
  581. IF (DVN.EQ.0.D0) THEN
  582. CALL ERREUR(21)
  583. SEGDES MELEME
  584. RETURN
  585. ENDIF
  586. XVN=XVN/DVN
  587. YVN=YVN/DVN
  588. ZVN=ZVN/DVN
  589. XNORM(1)=XVN
  590. XNORM(2)=YVN
  591. XNORM(3)=ZVN
  592. C Recherche des points verifiant la condition d'appartenance (INDIC)
  593. SEGINI INDIC
  594. NBON=0
  595. DO IP=1,NBELEM
  596. IPoin=NUM(1,IP)
  597. IRef=(IPoin-1)*idimp1
  598. XV=XCOOR(IRef+1)-XP1
  599. YV=XCOOR(IRef+2)-YP1
  600. ZV=0.D0
  601. IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1
  602. Val1=ABS(XV*XVN+YV*YVN+ZV*ZVN)
  603. IF (Val1.LE.CRIT) THEN
  604. INDIC(IP)=IPoin
  605. NBON=NBON+1
  606. ENDIF
  607. ENDDO
  608. C Intersection du maillage avec un CYLIndre
  609. C --------------------------
  610. ELSE IF (IMC.EQ.3) THEN
  611. CALL LIROBJ('POINT ',IPoin1,1,IRetou)
  612. CALL LIROBJ('POINT ',IPoin2,1,IRetou)
  613. CALL LIROBJ('POINT ',IPoin3,1,IRetou)
  614. IF (IERR.NE.0) THEN
  615. SEGDES MELEME
  616. RETURN
  617. ENDIF
  618. IRef1=(IPoin1-1)*idimp1
  619. IRef2=(IPoin2-1)*idimp1
  620. IRef3=(IPoin3-1)*idimp1
  621. XP1=XCOOR(IRef1+1)
  622. XVN=XCOOR(IRef2+1)-XP1
  623. XPT=XCOOR(IRef3+1)-XP1
  624. YP1=XCOOR(IRef1+2)
  625. YVN=XCOOR(IRef2+2)-YP1
  626. YPT=XCOOR(IRef3+2)-YP1
  627. IF (IDIM.GE.3) THEN
  628. ZP1=XCOOR(IRef1+3)
  629. ZVN=XCOOR(IRef2+3)-ZP1
  630. ZPT=XCOOR(IRef3+3)-ZP1
  631. ELSE
  632. ZP1=0.D0
  633. ZVN=0.D0
  634. ZPT=0.D0
  635. ENDIF
  636. C Determination de l'axe du cylindre (axe P1-P2)
  637. DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
  638. IF (DVN.EQ.0.D0) THEN
  639. CALL ERREUR(21)
  640. RETURN
  641. ENDIF
  642. XVN=XVN/DVN
  643. YVN=YVN/DVN
  644. ZVN=ZVN/DVN
  645. C Calcul du rayon du cylindre (distance P3 a axe P1-P2)
  646. SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
  647. XV=SCAL*XVN-XPT
  648. YV=SCAL*YVN-YPT
  649. ZV=SCAL*ZVN-ZPT
  650. Rayon=SQRT(XV*XV+YV*YV+ZV*ZV)
  651. IF (Rayon.EQ.0.D0) THEN
  652. CALL ERREUR(21)
  653. SEGDES MELEME
  654. RETURN
  655. ENDIF
  656. C Recherche des points verifiant la condition d'appartenance (INDIC)
  657. SEGINI INDIC
  658. NBON=0
  659. DO IP=1,NBELEM
  660. IPoin=NUM(1,IP)
  661. IRef=(IPoin-1)*idimp1
  662. XPT=XCOOR(IRef+1)-XP1
  663. YPT=XCOOR(IRef+2)-YP1
  664. ZPT=0.D0
  665. IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1
  666. SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
  667. XV=SCAL*XVN-XPT
  668. YV=SCAL*YVN-YPT
  669. ZV=SCAL*ZVN-ZPT
  670. Val1=SQRT(XV*XV+YV*YV+ZV*ZV)
  671. IF (ABS(Val1-Rayon).LE.CRIT) THEN
  672. INDIC(IP)=IPoin
  673. NBON=NBON+1
  674. ENDIF
  675. ENDDO
  676. C Intersection du maillage avec un CONE
  677. C --------------------------
  678. ELSE IF (IMC.EQ.4) THEN
  679. CALL LIROBJ('POINT ',IPoin1,1,IRetou)
  680. CALL LIROBJ('POINT ',IPoin2,1,IRetou)
  681. CALL LIROBJ('POINT ',IPoin3,1,IRetou)
  682. IF (IERR.NE.0) THEN
  683. SEGDES MELEME
  684. RETURN
  685. ENDIF
  686. IRef1=(IPoin1-1)*idimp1
  687. IRef2=(IPoin2-1)*idimp1
  688. IRef3=(IPoin3-1)*idimp1
  689. XP1=XCOOR(IRef1+1)
  690. XVN=XCOOR(IRef2+1)-XP1
  691. XV=XCOOR(IRef3+1)-XP1
  692. YP1=XCOOR(IRef1+2)
  693. YVN=XCOOR(IRef2+2)-YP1
  694. YV=XCOOR(IRef3+2)-YP1
  695. IF (IDIM.GE.3) THEN
  696. ZP1=XCOOR(IRef1+3)
  697. ZVN=XCOOR(IRef2+3)-ZP1
  698. ZV=XCOOR(IRef3+3)-ZP1
  699. ELSE
  700. ZP1=0.D0
  701. ZVN=0.D0
  702. ZV=0.D0
  703. ENDIF
  704. C Recherche du vecteur directeur de l'axe du cone
  705. DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
  706. IF (DVN.EQ.0.D0) THEN
  707. CALL ERREUR(21)
  708. SEGDES MELEME
  709. RETURN
  710. ENDIF
  711. XVN=XVN/DVN
  712. YVN=YVN/DVN
  713. ZVN=ZVN/DVN
  714. C Calcul de l'angle au sommet du cone
  715. Val1=XV*XV+YV*YV+ZV*ZV
  716. IF (Val1.EQ.0.D0) THEN
  717. CALL ERREUR(21)
  718. SEGDES MELEME
  719. RETURN
  720. ENDIF
  721. SCAL=XV*XVN+YV*YVN+ZV*ZVN
  722. XPT=SCAL*XVN-XV
  723. YPT=SCAL*YVN-YV
  724. ZPT=SCAL*ZVN-ZV
  725. SISOM=SQRT((XPT*XPT+YPT*YPT+ZPT*ZPT)/Val1)
  726. C Recherche des points verifiant la condition d'appartenance (INDIC)
  727. SEGINI INDIC
  728. NBON=0
  729. DO IP=1,NBELEM
  730. IPoin=NUM(1,IP)
  731. IRef=(IPoin-1)*idimp1
  732. XV=XCOOR(IRef+1)-XP1
  733. YV=XCOOR(IRef+2)-YP1
  734. ZV=0.D0
  735. IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1
  736. Val1=SQRT(XV*XV+YV*YV+ZV*ZV)*SISOM
  737. SCAL=XV*XVN+YV*YVN+ZV*ZVN
  738. XPT=SCAL*XVN-XV
  739. YPT=SCAL*YVN-YV
  740. ZPT=SCAL*ZVN-ZV
  741. Val2=SQRT(XPT*XPT+YPT*YPT+ZPT*ZPT)
  742. IF (ABS(Val2-Val1).LE.CRIT) THEN
  743. INDIC(IP)=IPoin
  744. NBON=NBON+1
  745. ENDIF
  746. ENDDO
  747. C Intersection du maillage avec une SPHERE
  748. C --------------------------
  749. ELSE IF (IMC.EQ.5) THEN
  750. CALL LIROBJ('POINT ',IPoin1,1,IRetou)
  751. CALL LIROBJ('POINT ',IPoin2,1,IRetou)
  752. IF (IERR.NE.0) THEN
  753. SEGDES MELEME
  754. RETURN
  755. ENDIF
  756. IRef1=(IPoin1-1)*idimp1
  757. IRef2=(IPoin2-1)*idimp1
  758. XP1=XCOOR(IRef1+1)
  759. XV=XCOOR(IRef2+1)-XP1
  760. YP1=XCOOR(IRef1+2)
  761. YV=XCOOR(IRef2+2)-YP1
  762. IF (IDIM.GE.3) THEN
  763. ZP1=XCOOR(IRef1+3)
  764. ZV=XCOOR(IRef2+3)-ZP1
  765. ELSE
  766. ZP1=0.D0
  767. ZV=0.D0
  768. ENDIF
  769. C Calcul du rayon de la sphere
  770. Rayon=SQRT(XV*XV+YV*YV+ZV*ZV)
  771. IF (Rayon.EQ.0.D0) THEN
  772. CALL ERREUR(21)
  773. SEGDES MELEME
  774. RETURN
  775. ENDIF
  776. C Recherche des points verifiant la condition d'appartenance (INDIC)
  777. SEGINI INDIC
  778. NBON=0
  779. DO IP=1,NBELEM
  780. IPoin=NUM(1,IP)
  781. IRef=(IPoin-1)*idimp1
  782. XV=XCOOR(IRef+1)-XP1
  783. YV=XCOOR(IRef+2)-YP1
  784. ZV=0.D0
  785. IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1
  786. Val1=SQRT(XV*XV+YV*YV+ZV*ZV)
  787. IF (ABS(Val1-Rayon).LE.CRIT) THEN
  788. INDIC(IP)=IPoin
  789. NBON=NBON+1
  790. ENDIF
  791. ENDDO
  792. C Intersection du maillage avec un TORE
  793. C --------------------------
  794. ELSE IF (IMC.EQ.6) THEN
  795. CALL LIROBJ('POINT ',IPoin1,1,IRetou)
  796. CALL LIROBJ('POINT ',IPoin2,1,IRetou)
  797. CALL LIROBJ('POINT ',IPoin3,1,IRetou)
  798. CALL LIROBJ('POINT ',IPoin4,1,IRetou)
  799. IF (IERR.NE.0) THEN
  800. SEGDES MELEME
  801. RETURN
  802. ENDIF
  803. IRef1=(IPoin1-1)*idimp1
  804. IRef2=(IPoin2-1)*idimp1
  805. IRef3=(IPoin3-1)*idimp1
  806. IRef4=(IPoin4-1)*idimp1
  807. XP1=XCOOR(IRef1+1)
  808. XVN=XCOOR(IRef2+1)-XP1
  809. XV=XCOOR(IRef3+1)-XP1
  810. XPT=XCOOR(IRef4+1)-XP1
  811. YP1=XCOOR(IRef1+2)
  812. YVN=XCOOR(IRef2+2)-YP1
  813. YV=XCOOR(IRef3+2)-YP1
  814. YPT=XCOOR(IRef4+2)-YP1
  815. IF (IDIM.GE.3) THEN
  816. ZP1=XCOOR(IRef1+3)
  817. ZVN=XCOOR(IRef2+3)-ZP1
  818. ZV=XCOOR(IRef3+3)-ZP1
  819. ZPT=XCOOR(IRef4+3)-ZP1
  820. ELSE
  821. ZP1=0.D0
  822. ZVN=0.D0
  823. ZV=0.D0
  824. ZPT=0.D0
  825. ENDIF
  826. C Determination du vecteur directeur de l'axe du tore
  827. DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
  828. IF (DVN.EQ.0.D0) THEN
  829. CALL ERREUR(21)
  830. SEGDES MELEME
  831. RETURN
  832. ENDIF
  833. XVN=XVN/DVN
  834. YVN=YVN/DVN
  835. ZVN=ZVN/DVN
  836. C Calcul du grand rayon du tore
  837. GRayon=XV*XV+YV*YV+ZV*ZV
  838. IF (GRayon.EQ.0.D0) THEN
  839. CALL ERREUR(21)
  840. SEGDES MELEME
  841. RETURN
  842. ENDIF
  843. C Calcul du petit rayon du tore
  844. SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
  845. XV=XPT-SCAL*XVN
  846. YV=YPT-SCAL*YVN
  847. ZV=ZPT-SCAL*ZVN
  848. Val1=SQRT(Grayon/(XV*XV+YV*YV+ZV*ZV))
  849. XV=XV*Val1-XPT
  850. YV=YV*Val1-YPT
  851. ZV=ZV*Val1-ZPT
  852. PRayon=SQRT(XV*XV+YV*YV+ZV*ZV)
  853. IF (PRayon.EQ.0.D0) THEN
  854. CALL ERREUR(21)
  855. SEGDES MELEME
  856. RETURN
  857. ENDIF
  858. C Recherche des points verifiant la condition d'appartenance (INDIC)
  859. SEGINI INDIC
  860. NBON=0
  861. DO IP=1,NBELEM
  862. IPoin=NUM(1,IP)
  863. IRef=(IPoin-1)*idimp1
  864. XPT=XCOOR(IRef+1)-XP1
  865. YPT=XCOOR(IRef+2)-YP1
  866. ZPT=0.D0
  867. IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1
  868. SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
  869. XV=XPT-SCAL*XVN
  870. YV=YPT-SCAL*YVN
  871. ZV=ZPT-SCAL*ZVN
  872. Val1=SQRT(GRayon/(XV*XV+YV*YV+ZV*ZV))
  873. XV=XV*Val1-XPT
  874. YV=YV*Val1-YPT
  875. ZV=ZV*Val1-ZPT
  876. Val2=SQRT(XV*XV+YV*YV+ZV*ZV)
  877. IF (ABS(Val2-PRayon).LE.CRIT) THEN
  878. INDIC(IP)=IPoin
  879. NBON=NBON+1
  880. ENDIF
  881. ENDDO
  882. ENDIF
  883. C Creation du maillage de sortie
  884. C ----------------------
  885. IF (NBON.EQ.0) THEN
  886. CALL ERREUR(18)
  887. SEGDES MELEME
  888. SEGSUP INDIC
  889. RETURN
  890. ENDIF
  891. NBEL=NBELEM
  892. NBNN=1
  893. NBELEM=NBON
  894. NBSOUS=0
  895. NBREF=0
  896. SEGINI IPT1
  897. IPT1.ITYPEL=1
  898. IPLAC=0
  899. DO IP=1,NBEL
  900. IF (INDIC(IP).NE.0) THEN
  901. IPLAC=IPLAC+1
  902. IPT1.NUM(1,IPLAC)=INDIC(IP)
  903. IPT1.ICOLOR(IPLAC)=IDCOUL
  904. ENDIF
  905. ENDDO
  906. SEGSUP INDIC
  907. SEGDES MELEME,IPT1
  908. CALL ECROBJ('MAILLAGE',IPT1)
  909.  
  910. RETURN
  911.  
  912.  
  913.  
  914. 300 CONTINUE
  915. c IF(IIMPI.EQ.333) WRITE(IOIMP,*) 'POINT syntaxe 3 :'
  916. C -------------------------------
  917. C Operateur POINT - syntaxe 3 :
  918. C -------------------------------
  919. C Extraction d'un CHAML ou d'un CHPO du(des) point(s) verifiant une
  920. C condition indiquee par MOT1.
  921. C Geo1 = POINT Che1 Mot1 ...
  922.  
  923. Val1=0.D0
  924. Val2=0.D0
  925. IPList=0
  926. IF (IMM.GT.2) THEN
  927. CALL LIRREE(Val1,1,IRetou)
  928. IF (IERR.NE.0) RETURN
  929. IF (IMM.EQ.9) THEN
  930. CALL LIRREE(Val2,1,IRetou)
  931. IF (IERR.NE.0) RETURN
  932. ENDIF
  933. ENDIF
  934. CALL LIRMOT(MOABS,1,IAB,0)
  935. CALL LIRMOT(MOTAV,2,IAV,0)
  936. IF (IAV.EQ.0) IAV=1
  937. CALL LIROBJ('LISTMOTS',IPList,0,IRetou)
  938. IF (IERR.NE.0) RETURN
  939. CALL LIROBJ('CHPOINT ',IPCHP,0,IRetou)
  940. IF (IERR.NE.0) RETURN
  941. IF (IRetou.NE.0) THEN
  942. CALL EXPCHP(IPCHP,IMM,IAB,IAV,IPList,Val1,Val2,IPMAIL)
  943. ELSE
  944. CALL LIROBJ('MCHAML ',IPCHE,1,IRetou)
  945. IF (IERR.NE.0) RETURN
  946. CALL EXPCHE(IPCHE,IMM,IAB,IAV,IPList,Val1,Val2,IPMAIL)
  947. ENDIF
  948. IF (IERR.EQ.0) CALL ECROBJ('MAILLAGE',IPMAIL)
  949.  
  950. RETURN
  951.  
  952. END
  953.  
  954.  
  955.  
  956.  

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