Télécharger la plaquette Cast3M


fr en

Télécharger poiext.eso

Retour à la liste

Numérotation des lignes : 

poiext

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

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