Esope source files

   1 : C POIEXT    SOURCE    GOUNAND   13/08/01    21:15:07     7810           
   2 : 
   3 : C  Sert a extraire un point d'un objet MAILLAGE
   4 : C  Novembre 1985 :
   5 : C    Extrait d'un objet tous les points appartenant a une DROIte ou a
   6 : C    un PLAN (possibilite de critere de distance)
   7 : C  Fevrier 1986 :
   8 : C    Extrait d'un objet tous les points appartenant a un CYLIndre, un
   9 : C    CONE, une SPHEre ou un TORE
  10 : C  09/2003 :
  11 : C    Introduction de la DIMEnsion 1. Ajout de la syntaxe pour creer un
  12 : C    point a partir de ses coordonnees. Inhibation d'options non compa-
  13 : C    tibles avec la dimension 1. Option PROC adapte au cas IDIM=1.
  14 : C    Ajout de tests sur le type d'element pour POINT Entier,
  15 : C    POINT INITIAL et POINT FINAL (conforme a la notice).
  16 : C  01/2004 :
  17 : C    Modifications des methodes de recherche des points appartenant a
  18 : C    une DROIte, un PLAN, un CYLIndre, un CONE. Methodes moins sensibles
  19 : C    aux imprecisions numeriques.
  20 : C  07/2013 GOUNAND
  21 : C    Les options INIT FINA et PROC détectent les maillages vides
  22 : C    en entrée et provoquent une erreur
  23 : C
  24 :       SUBROUTINE POIEXT
  25 : 
  26 :       IMPLICIT INTEGER(I-N)
  27 :       IMPLICIT REAL*8 (A-H,O-Z)
  28 : 
  29 : -INC CCOPTIO
  30 : -INC SMELEME
  31 : -INC SMLREEL
  32 : -INC SMCOORD
  33 : -INC CCGEOME
  34 : 
  35 : C-    DIMENSION XNORM(3)
  36 :       CHARACTER*4 MOT,IPRO
  37 :       CHARACTER*4 MCLE(9),MOTM(9),MOABS(1),MOTAV(2)
  38 :       SEGMENT INDIC(NBELEM)
  39 : 
  40 : C  COMMON de sauvegarde lors d'un appel par SYMTRI
  41 :       COMMON / CSYMTR / XNORM(3)
  42 : 
  43 :       DATA MCLE / 'INIT','FINA','PROC',
  44 :      .     'DROI','PLAN','CYLI','CONE','SPHE','TORE' /
  45 :       DATA MOTM / 'MAXI','MINI','SUPE','EGSU','EGAL','EGIN','INFE',
  46 :      .     'DIFF','COMP' /
  47 :       DATA MOTAV / 'AVEC','SANS' /
  48 :       DATA MOABS / 'ABS ' /
  49 : 
  49 :       CALL LIRMOT(MOTM,9,IMM,0)
  50 :       IF (IMM.EQ.0) THEN
  51 :          idimp1=IDIM+1
  51 :          CALL LIRMOT(MCLE(4),6,IMC,0)
  52 :          IF (IMC.EQ.0) THEN
  52 :             CALL LIROBJ('MAILLAGE',MELEME,0,IRetou)
  53 :             IF (IERR.NE.0) RETURN
  54 : C  Operateur POINT - syntaxe 1 :
  55 : C -------------------------------
  56 : C  Creation du POINT1 a partir des IDIM coordonnees fournies.
  57 : C  Pour IDIM = 1, seul moyen de creer un point.
  58 : C    Point1 = POINT Val1 (Val2) (Val3) ;
  59 :             IF (IRetou.EQ.0) THEN
  59 :                CALL LIROBJ('LISTREEL',MLREEL,0,IRETOU)
  60 :                IF (IRETOU.NE.0) THEN
  61 :                   SEGACT MLREEL
  62 :                   NR=PROG(/1)
  63 :                   SEGACT MCOORD*MOD
  64 :                   NbPts=1+(XCOOR(/1)/idimp1)
  65 :                   SEGADJ MCOORD
  66 :                   IRef=(NbPts-1)*idimp1
  67 :                   DO IIDIM=1,IDIM
  68 :                      IF (IIDIM.LE.NR) THEN 
  69 :                         XVAL=PROG(IIDIM)
  70 :                      ELSE 
  71 :                         XVAL=0.D0
  72 :                      ENDIF
  73 :                      XCOOR(IRef+IIDIM)=XVAL
  74 :                   ENDDO
  75 :                   SEGDES MLREEL
  76 :                ELSE
  76 :                   CALL MESLIR(-149)
  76 :                   CALL LIRREE(Val1,1,IRetou)
  76 :                   CALL MESLIR(-150)
  77 :                   IF (IDIM.EQ.2) THEN
  77 :                      CALL LIRREE(Val2,1,IRetou)
  78 :                   ELSE IF (IDIM.EQ.3) THEN
  78 :                      CALL LIRREE(Val2,1,IRetou)
  79 :                      IF (IERR.NE.0) RETURN
  79 :                      CALL MESLIR(-151)
  79 :                      CALL LIRREE(Val3,1,IRetou)
  80 :                   ENDIF
  81 :                   IF (IERR.NE.0) RETURN
  82 :                   SEGACT MCOORD*MOD
  83 :                   NbPts=1+(XCOOR(/1)/idimp1)
  84 :                   SEGADJ MCOORD
  85 :                   IRef=(NbPts-1)*idimp1
  86 :                   XCOOR(IRef+1)=Val1
  87 :                   IF (IDIM.EQ.2) THEN
  88 :                      XCOOR(IRef+2)=Val2
  89 :                   ELSE IF (IDIM.EQ.3) THEN
  90 :                      XCOOR(IRef+2)=Val2
  91 :                      XCOOR(IRef+3)=Val3
  92 :                   ENDIF
  93 :                ENDIF
  94 :                XCOOR(IRef+idimp1)=DENSIT
  94 :                CALL ECROBJ('POINT   ',NbPts)
  95 :                RETURN
  96 : C  Operateur POINT - syntaxe 2 :
  97 : C -------------------------------
  98 : C  1) Recuperation du POINT1 du maillage GEO1 par son numero
  99 : C       Point1 = POINT Geo1 Entier1 ;
 100 : C     Note : le maillage est constitue uniquement de POI1 ou SEG2.
 101 : C  2) Recuperation du point INITIAL ou FINAL du maillage GEO1
 102 : C       Point1 = POINT Geo1 'INITIAL' | 'FINAL' ;
 103 : C     Note : le maillage est constitue uniquement de POI1, SEG2 ou SEG3.
 104 : C  3) Recuperation du POINT1 du maillage GEO1 le plus proche du POINT2.
 105 : C       Point1 = POINT Geo1 'PROC' Point2 ;
 106 :             ELSE
 107 : C  Lecture facultative d'un entier (numero du point a trouver)
 107 :                CALL LIRENT(JPTR,0,IRetou)
 108 :                IF (IRetou.NE.0) THEN
 109 :                   IF (JPTR.LE.0) THEN
 109 :                      CALL ERREUR(262)
 110 :                      RETURN
 111 :                   ENDIF
 112 :                   MOT='    '
 113 :                   IPRO='    '
 114 :                ELSE
 115 : C  Lecture d'un des mots cles ('INIT','FINA','PROC')
 115 :                   CALL LIRCHA(MOT,0,IRetou)
 116 :                   IF (IRetou.EQ.0) THEN
 117 : C  Le type de l'operande est incorrect
 117 :                      CALL QUETYP(MOTERR(1:8),0,IRetou)
 118 :                      IF (IRetou.NE.0) THEN
 118 :                         CALL ERREUR (39)
 119 :                      ELSE
 119 :                         CALL ERREUR(533)
 120 :                      ENDIF
 121 :                      RETURN
 122 :                   ENDIF
 123 : C  Recuperation du numero du point INIT ou FINA
 124 :                   IF ((MOT.EQ.MCLE(1)).OR.(MOT.EQ.MCLE(2))) THEN
 125 :                      IPRO='    '
 126 :                   ELSE IF (MOT.EQ.MCLE(3)) THEN
 127 :                      IPRO=MOT
 128 :                   ELSE
 128 :                      CALL ERREUR(7)
 129 :                      RETURN
 130 :                   ENDIF
 131 :                ENDIF
 132 :                SEGACT MELEME
 133 : C  Traitement des syntaxes POINT Geo1 Entier1 | 'INIT' | 'FINA'
 134 :                IF (IPRO.EQ.'    ') THEN
 135 : C  Erreur si objet maillage complexe
 136 :                   IF (LISOUS(/1).NE.0) THEN
 136 :                      CALL ERREUR(25)
 137 :                      SEGDES MELEME
 138 :                      RETURN
 139 :                   ENDIF
 140 :                   IF (NUM(/1).LE.0.OR.NUM(/2).LE.0) THEN
 141 :                      MOTERR(1:8)='MAILLAGE'
 142 : C Une donnée de type %m1:8 est vide
 142 :                      CALL ERREUR(1027)
 143 :                      SEGDES MELEME
 144 :                      RETURN
 145 :                   ENDIF
 146 : C                     
 147 : C  Erreur si objet maillage n'est pas de type POI1,SEG2,SEG3
 148 : Cde           IF ((ITYPEL.NE.1).AND.(ITYPEL.NE.2).AND.
 149 : Cde  .            (ITYPEL.NE.3)) THEN
 149 : Cde             CALL ERREUR(16)
 150 : Cde             SEGDES MELEME
 151 : Cde             RETURN
 152 : Cde           ENDIF
 153 :                   IF (MOT.EQ.'    ') THEN
 154 : Cde             IF (ITYPEL.EQ.3) THEN
 154 : Cde               CALL ERREUR(16)
 155 : Cde               SEGDES MELEME
 156 : Cde               RETURN
 157 : Cde             ENDIF
 158 : Cde             IF ((ITYPEL.EQ.2).AND.(JPTR.EQ.(NUM(/2)+1))) THEN
 159 :                      IF ((ITYPEL.NE.1).AND.(JPTR.EQ.(NUM(/2)+1))) THEN
 160 :                         IPTR=NUM(/1)
 161 :                         JPTR=NUM(/2)
 162 :                      ELSE
 163 :                         IPTR=1
 164 :                      ENDIF
 165 :                      IF (JPTR.GT.NUM(/2)) THEN
 165 :                         CALL ERREUR(262)
 166 :                         SEGDES MELEME
 167 :                         RETURN
 168 :                      ENDIF
 169 :                   ELSE IF (MOT.EQ.MCLE(1)) THEN
 170 :                      IPTR=1
 171 :                      JPTR=1
 172 :                   ELSE IF (MOT.EQ.MCLE(2)) THEN
 173 :                      IPTR=NUM(/1)
 174 :                      JPTR=NUM(/2)
 175 :                   ENDIF
 176 :                   IPoin=NUM(IPTR,JPTR)
 176 :                   CALL ECROBJ('POINT ',IPoin)
 177 :                   SEGDES MELEME
 178 :                   RETURN
 179 :                ELSE
 180 : C  Traitement de l'option PROC
 180 :                   CALL LIROBJ('POINT ',IPoin1,1,IRetou)
 181 :                   IF (IERR.NE.0) THEN
 182 :                      SEGDES MELEME
 183 :                      RETURN
 184 :                   ENDIF
 185 :                   Val1=1.E30
 186 :                   IRef=(IPoin1-1)*idimp1
 187 :                   XP1=XCOOR(IRef+1)
 188 :                   YP1=XCOOR(IRef+2)
 189 :                   ZP1=0.D0
 190 :                   IF (IDIM.GE.3) ZP1=XCOOR(IRef+3)
 191 :                   IPT1=MELEME
 192 :                   Ipoin=0
 193 :                   DO IOB=1,MAX(1,LISOUS(/1))
 194 :                      IF (LISOUS(/1).NE.0) THEN
 195 :                         IPT1=LISOUS(IOB)
 196 :                         SEGACT IPT1
 197 :                      ENDIF
 198 :                      DO i=1,IPT1.NUM(/2)
 199 :                         DO j=1,IPT1.NUM(/1)
 200 :                            IRef=(IPT1.NUM(j,i)-1)*idimp1
 201 :                            IF (IDIM.EQ.1) THEN
 202 :                               DVN=ABS(XP1-XCOOR(IRef+1))
 203 :                            ELSE
 204 :                               XVN=XP1-XCOOR(IRef+1)
 205 :                               YVN=YP1-XCOOR(IRef+2)
 206 :                               ZVN=0.D0
 207 :                               IF (IDIM.GE.3) ZVN=ZP1-XCOOR(IRef+3)
 208 :                               DVN=XVN*XVN+YVN*YVN+ZVN*ZVN
 209 :                            ENDIF
 210 :                            IF (DVN.LE.Val1) THEN
 211 :                               Val1=DVN
 212 :                               IPoin=IPT1.NUM(j,i)
 213 :                            ENDIF
 214 :                         ENDDO
 215 :                      ENDDO
 216 :                      IF (LISOUS(/1).NE.0) SEGDES IPT1
 217 :                   ENDDO
 218 :                   IF (IPoin.EQ.0) THEN
 219 :                      MOTERR(1:8)='MAILLAGE'
 220 : C Une donnée de type %m1:8 est vide
 220 :                      CALL ERREUR(1027)
 221 :                   ELSE
 221 :                      IF (IERR.EQ.0) CALL ECROBJ('POINT ',IPoin)
 222 :                   ENDIF
 223 :                   SEGDES MELEME
 224 :                   RETURN
 225 :                ENDIF
 226 :             ENDIF
 227 :          ELSE
 228 : C  Operateur POINT - syntaxe 2 :
 229 : C -------------------------------
 230 : C  Extraction du(des) point(s) du maillage GEO1 situe(s) sur une
 231 : C  DROIte, un PLAN, un CYLIndre, un CONE, une SPHEre ou un TORE.
 232 : C  L'appartenance du point a l'objet est definie suivant un critere
 233 : C  de distance non nul (appartenance stricte : critere tres faible).
 234 : C    Geof = POINT Geo1 MOT1 Points ... (Critere)
 235 : C  En dimension 1, ces criteres n'ont aucun sens.
 236 :             IF (IDIM.EQ.1) THEN
 237 :                MOTERR(1:4)=MCLE(3+IMC)
 238 :                INTERR(1)=IDIM
 238 :                CALL ERREUR(971)
 239 :                RETURN
 240 :             ENDIF
 241 : C  Lectures obligatoire du maillage et facultative du critere
 241 :             CALL LIROBJ('MAILLAGE',MELEME,1,IRetou)
 242 :             IF (IERR.NE.0) RETURN
 243 :             CRIT=DENSIT/10.D0
 243 :             CALL LIRREE(CRIT,0,IRetou)
 243 :             IF (CRIT.LE.0.D0) CALL ERREUR(17)
 244 :             IF (IERR.NE.0) RETURN
 245 : C  Transformation du maillage en 'POI1'
 246 :             SEGACT MELEME
 246 :             CALL CHANGE(MELEME,1)
 247 :             SEGACT MELEME
 248 :             NBELEM=NUM(/2)
 249 : C  Intersection du maillage avec une DROIte
 250 : C --------------------------
 251 :             IF (IMC.EQ.1) THEN
 251 :                CALL LIROBJ('POINT ',IPoin1,1,IRetou)
 251 :                CALL LIROBJ('POINT ',IPoin2,1,IRetou)
 252 :                IF (IERR.NE.0) THEN
 253 :                   SEGDES MELEME
 254 :                   RETURN
 255 :                ENDIF
 256 : C  Recherche du vecteur directeur de la DROIte
 257 :                IRef1=(IPoin1-1)*idimp1
 258 :                IRef2=(IPoin2-1)*idimp1
 259 :                XP1=XCOOR(IRef1+1)
 260 :                XVN=XCOOR(IRef2+1)-XP1
 261 :                YP1=XCOOR(IRef1+2)
 262 :                YVN=XCOOR(IRef2+2)-YP1
 263 :                ZP1=0.D0
 264 :                ZVN=0.D0
 265 :                IF (IDIM.GE.3) THEN
 266 :                   ZP1=XCOOR(IRef1+3)
 267 :                   ZVN=XCOOR(IRef2+3)-ZP1
 268 :                ENDIF
 269 :                DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
 270 :                IF (DVN.EQ.0.D0) THEN
 270 :                   CALL ERREUR(21)
 271 :                   SEGDES MELEME
 272 :                   RETURN
 273 :                ENDIF
 274 :                XVN=XVN/DVN
 275 :                YVN=YVN/DVN
 276 :                ZVN=ZVN/DVN
 277 :                XNORM(1)=XVN
 278 :                XNORM(2)=YVN
 279 :                XNORM(3)=ZVN
 280 : C  Recherche des points verifiant la condition d'appartenance (INDIC)
 281 :                CRIT=CRIT*CRIT
 282 :                SEGINI INDIC
 283 :                NBON=0
 284 :                DO IP=1,NBELEM
 285 :                   IPoin=NUM(1,IP)
 286 :                   IRef=(IPoin-1)*idimp1
 287 :                   XPT=XCOOR(IRef+1)-XP1
 288 :                   YPT=XCOOR(IRef+2)-YP1
 289 :                   ZPT=0.D0
 290 :                   IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1
 291 :                   SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
 292 :                   XV=SCAL*XVN-XPT
 293 :                   YV=SCAL*YVN-YPT
 294 :                   ZV=SCAL*ZVN-ZPT
 295 :                   Val1=XV*XV+YV*YV+ZV*ZV
 296 :                   IF (Val1.LE.CRIT) THEN
 297 :                      INDIC(IP)=IPoin
 298 :                      NBON=NBON+1
 299 :                   ENDIF
 300 :                ENDDO
 301 : C  Intersection du maillage avec un PLAN
 302 : C --------------------------
 303 :             ELSE IF (IMC.EQ.2) THEN
 303 :                CALL LIROBJ('POINT ',IPoin1,1,IRetou)
 303 :                CALL LIROBJ('POINT ',IPoin2,1,IRetou)
 303 :                CALL LIROBJ('POINT ',IPoin3,1,IRetou)
 304 :                IF (IERR.NE.0) THEN
 305 :                   SEGDES MELEME
 306 :                   RETURN
 307 :                ENDIF
 308 : C  Determination du vecteur normal au PLAN
 309 :                IRef1=(IPoin1-1)*idimp1
 310 :                IRef2=(IPoin2-1)*idimp1
 311 :                IRef3=(IPoin3-1)*idimp1
 312 :                XP1=XCOOR(IRef1+1)
 313 :                XP2=XCOOR(IRef2+1)
 314 :                XP3=XCOOR(IRef3+1)
 315 :                YP1=XCOOR(IRef1+2)
 316 :                YP2=XCOOR(IRef2+2)
 317 :                YP3=XCOOR(IRef3+2)
 318 :                ZP1=0.D0
 319 :                ZP2=0.D0
 320 :                ZP3=0.D0
 321 :                IF (IDIM.GE.3) THEN
 322 :                   ZP1=XCOOR(IRef1+3)
 323 :                   ZP2=XCOOR(IRef2+3)
 324 :                   ZP3=XCOOR(IRef3+3)
 325 :                ENDIF
 326 :                XVN=(YP2-YP1)*(ZP3-ZP1)-(ZP2-ZP1)*(YP3-YP1)
 327 :                YVN=(ZP2-ZP1)*(XP3-XP1)-(XP2-XP1)*(ZP3-ZP1)
 328 :                ZVN=(XP2-XP1)*(YP3-YP1)-(YP2-YP1)*(XP3-XP1)
 329 :                DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
 330 :                IF (DVN.EQ.0.D0) THEN
 330 :                   CALL ERREUR(21)
 331 :                   SEGDES MELEME
 332 :                   RETURN
 333 :                ENDIF
 334 :                XVN=XVN/DVN
 335 :                YVN=YVN/DVN
 336 :                ZVN=ZVN/DVN
 337 :                XNORM(1)=XVN
 338 :                XNORM(2)=YVN
 339 :                XNORM(3)=ZVN
 340 : C  Recherche des points verifiant la condition d'appartenance (INDIC)
 341 :                SEGINI INDIC
 342 :                NBON=0
 343 :                DO IP=1,NBELEM
 344 :                   IPoin=NUM(1,IP)
 345 :                   IRef=(IPoin-1)*idimp1
 346 :                   XV=XCOOR(IRef+1)-XP1
 347 :                   YV=XCOOR(IRef+2)-YP1
 348 :                   ZV=0.D0
 349 :                   IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1
 350 :                   Val1=ABS(XV*XVN+YV*YVN+ZV*ZVN)
 351 :                   IF (Val1.LE.CRIT) THEN
 352 :                      INDIC(IP)=IPoin
 353 :                      NBON=NBON+1
 354 :                   ENDIF
 355 :                ENDDO
 356 : C  Intersection du maillage avec un CYLIndre
 357 : C --------------------------
 358 :             ELSE IF (IMC.EQ.3) THEN
 358 :                CALL LIROBJ('POINT ',IPoin1,1,IRetou)
 358 :                CALL LIROBJ('POINT ',IPoin2,1,IRetou)
 358 :                CALL LIROBJ('POINT ',IPoin3,1,IRetou)
 359 :                IF (IERR.NE.0) THEN
 360 :                   SEGDES MELEME
 361 :                   RETURN
 362 :                ENDIF
 363 :                IRef1=(IPoin1-1)*idimp1
 364 :                IRef2=(IPoin2-1)*idimp1
 365 :                IRef3=(IPoin3-1)*idimp1
 366 :                XP1=XCOOR(IRef1+1)
 367 :                XVN=XCOOR(IRef2+1)-XP1
 368 :                XPT=XCOOR(IRef3+1)-XP1
 369 :                YP1=XCOOR(IRef1+2)
 370 :                YVN=XCOOR(IRef2+2)-YP1
 371 :                YPT=XCOOR(IRef3+2)-YP1
 372 :                IF (IDIM.GE.3) THEN
 373 :                   ZP1=XCOOR(IRef1+3)
 374 :                   ZVN=XCOOR(IRef2+3)-ZP1
 375 :                   ZPT=XCOOR(IRef3+3)-ZP1
 376 :                ELSE
 377 :                   ZP1=0.D0
 378 :                   ZVN=0.D0
 379 :                   ZPT=0.D0
 380 :                ENDIF
 381 : C  Determination de l'axe du cylindre (axe P1-P2)
 382 :                DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
 383 :                IF (DVN.EQ.0.D0) THEN
 383 :                   CALL ERREUR(21)
 384 :                   RETURN
 385 :                ENDIF
 386 :                XVN=XVN/DVN
 387 :                YVN=YVN/DVN
 388 :                ZVN=ZVN/DVN
 389 : C  Calcul du rayon du cylindre (distance P3 a axe P1-P2)
 390 :                SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
 391 :                XV=SCAL*XVN-XPT
 392 :                YV=SCAL*YVN-YPT
 393 :                ZV=SCAL*ZVN-ZPT
 394 :                Rayon=SQRT(XV*XV+YV*YV+ZV*ZV)
 395 :                IF (Rayon.EQ.0.D0) THEN
 395 :                   CALL ERREUR(21)
 396 :                   SEGDES MELEME
 397 :                   RETURN
 398 :                ENDIF
 399 : C  Recherche des points verifiant la condition d'appartenance (INDIC)
 400 :                SEGINI INDIC
 401 :                NBON=0
 402 :                DO IP=1,NBELEM
 403 :                   IPoin=NUM(1,IP)
 404 :                   IRef=(IPoin-1)*idimp1
 405 :                   XPT=XCOOR(IRef+1)-XP1
 406 :                   YPT=XCOOR(IRef+2)-YP1
 407 :                   ZPT=0.D0
 408 :                   IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1
 409 :                   SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
 410 :                   XV=SCAL*XVN-XPT
 411 :                   YV=SCAL*YVN-YPT
 412 :                   ZV=SCAL*ZVN-ZPT
 413 :                   Val1=SQRT(XV*XV+YV*YV+ZV*ZV)
 414 :                   IF (ABS(Val1-Rayon).LE.CRIT) THEN
 415 :                      INDIC(IP)=IPoin
 416 :                      NBON=NBON+1
 417 :                   ENDIF
 418 :                ENDDO
 419 : C  Intersection du maillage avec un CONE
 420 : C --------------------------
 421 :             ELSE IF (IMC.EQ.4) THEN
 421 :                CALL LIROBJ('POINT ',IPoin1,1,IRetou)
 421 :                CALL LIROBJ('POINT ',IPoin2,1,IRetou)
 421 :                CALL LIROBJ('POINT ',IPoin3,1,IRetou)
 422 :                IF (IERR.NE.0) THEN
 423 :                   SEGDES MELEME
 424 :                   RETURN
 425 :                ENDIF
 426 :                IRef1=(IPoin1-1)*idimp1
 427 :                IRef2=(IPoin2-1)*idimp1
 428 :                IRef3=(IPoin3-1)*idimp1
 429 :                XP1=XCOOR(IRef1+1)
 430 :                XVN=XCOOR(IRef2+1)-XP1
 431 :                XV=XCOOR(IRef3+1)-XP1
 432 :                YP1=XCOOR(IRef1+2)
 433 :                YVN=XCOOR(IRef2+2)-YP1
 434 :                YV=XCOOR(IRef3+2)-YP1
 435 :                IF (IDIM.GE.3) THEN
 436 :                   ZP1=XCOOR(IRef1+3)
 437 :                   ZVN=XCOOR(IRef2+3)-ZP1
 438 :                   ZV=XCOOR(IRef3+3)-ZP1
 439 :                ELSE
 440 :                   ZP1=0.D0
 441 :                   ZVN=0.D0
 442 :                   ZV=0.D0
 443 :                ENDIF
 444 : C  Recherche du vecteur directeur de l'axe du cone
 445 :                DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
 446 :                IF (DVN.EQ.0.D0) THEN
 446 :                   CALL ERREUR(21)
 447 :                   SEGDES MELEME
 448 :                   RETURN
 449 :                ENDIF
 450 :                XVN=XVN/DVN
 451 :                YVN=YVN/DVN
 452 :                ZVN=ZVN/DVN
 453 : C  Calcul de l'angle au sommet du cone
 454 :                Val1=XV*XV+YV*YV+ZV*ZV
 455 :                IF (Val1.EQ.0.D0) THEN
 455 :                   CALL ERREUR(21)
 456 :                   SEGDES MELEME
 457 :                   RETURN
 458 :                ENDIF
 459 :                SCAL=XV*XVN+YV*YVN+ZV*ZVN
 460 :                XPT=SCAL*XVN-XV
 461 :                YPT=SCAL*YVN-YV
 462 :                ZPT=SCAL*ZVN-ZV
 463 :                SISOM=SQRT((XPT*XPT+YPT*YPT+ZPT*ZPT)/Val1)
 464 : C  Recherche des points verifiant la condition d'appartenance (INDIC)
 465 :                SEGINI INDIC
 466 :                NBON=0
 467 :                DO IP=1,NBELEM
 468 :                   IPoin=NUM(1,IP)
 469 :                   IRef=(IPoin-1)*idimp1
 470 :                   XV=XCOOR(IRef+1)-XP1
 471 :                   YV=XCOOR(IRef+2)-YP1
 472 :                   ZV=0.D0
 473 :                   IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1
 474 :                   Val1=SQRT(XV*XV+YV*YV+ZV*ZV)*SISOM
 475 :                   SCAL=XV*XVN+YV*YVN+ZV*ZVN
 476 :                   XPT=SCAL*XVN-XV
 477 :                   YPT=SCAL*YVN-YV
 478 :                   ZPT=SCAL*ZVN-ZV
 479 :                   Val2=SQRT(XPT*XPT+YPT*YPT+ZPT*ZPT)
 480 :                   IF (ABS(Val2-Val1).LE.CRIT) THEN
 481 :                      INDIC(IP)=IPoin
 482 :                      NBON=NBON+1
 483 :                   ENDIF
 484 :                ENDDO
 485 : C  Intersection du maillage avec une SPHERE
 486 : C --------------------------
 487 :             ELSE IF (IMC.EQ.5) THEN
 487 :                CALL LIROBJ('POINT ',IPoin1,1,IRetou)
 487 :                CALL LIROBJ('POINT ',IPoin2,1,IRetou)
 488 :                IF (IERR.NE.0) THEN
 489 :                   SEGDES MELEME
 490 :                   RETURN
 491 :                ENDIF
 492 :                IRef1=(IPoin1-1)*idimp1
 493 :                IRef2=(IPoin2-1)*idimp1
 494 :                XP1=XCOOR(IRef1+1)
 495 :                XV=XCOOR(IRef2+1)-XP1
 496 :                YP1=XCOOR(IRef1+2)
 497 :                YV=XCOOR(IRef2+2)-YP1
 498 :                IF (IDIM.GE.3) THEN
 499 :                   ZP1=XCOOR(IRef1+3)
 500 :                   ZV=XCOOR(IRef2+3)-ZP1
 501 :                ELSE
 502 :                   ZP1=0.D0
 503 :                   ZV=0.D0
 504 :                ENDIF
 505 : C  Calcul du rayon de la sphere
 506 :                Rayon=SQRT(XV*XV+YV*YV+ZV*ZV)
 507 :                IF (Rayon.EQ.0.D0) THEN
 507 :                   CALL ERREUR(21)
 508 :                   SEGDES MELEME
 509 :                   RETURN
 510 :                ENDIF
 511 : C  Recherche des points verifiant la condition d'appartenance (INDIC)
 512 :                SEGINI INDIC
 513 :                NBON=0
 514 :                DO IP=1,NBELEM
 515 :                   IPoin=NUM(1,IP)
 516 :                   IRef=(IPoin-1)*idimp1
 517 :                   XV=XCOOR(IRef+1)-XP1
 518 :                   YV=XCOOR(IRef+2)-YP1
 519 :                   ZV=0.D0
 520 :                   IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1
 521 :                   Val1=SQRT(XV*XV+YV*YV+ZV*ZV)
 522 :                   IF (ABS(Val1-Rayon).LE.CRIT) THEN
 523 :                      INDIC(IP)=IPoin
 524 :                      NBON=NBON+1
 525 :                   ENDIF
 526 :                ENDDO
 527 : C  Intersection du maillage avec un TORE
 528 : C --------------------------
 529 :             ELSE IF (IMC.EQ.6) THEN
 529 :                CALL LIROBJ('POINT ',IPoin1,1,IRetou)
 529 :                CALL LIROBJ('POINT ',IPoin2,1,IRetou)
 529 :                CALL LIROBJ('POINT ',IPoin3,1,IRetou)
 529 :                CALL LIROBJ('POINT ',IPoin4,1,IRetou)
 530 :                IF (IERR.NE.0) THEN
 531 :                   SEGDES MELEME
 532 :                   RETURN
 533 :                ENDIF
 534 :                IRef1=(IPoin1-1)*idimp1
 535 :                IRef2=(IPoin2-1)*idimp1
 536 :                IRef3=(IPoin3-1)*idimp1
 537 :                IRef4=(IPoin4-1)*idimp1
 538 :                XP1=XCOOR(IRef1+1)
 539 :                XVN=XCOOR(IRef2+1)-XP1
 540 :                XV=XCOOR(IRef3+1)-XP1
 541 :                XPT=XCOOR(IRef4+1)-XP1
 542 :                YP1=XCOOR(IRef1+2)
 543 :                YVN=XCOOR(IRef2+2)-YP1
 544 :                YV=XCOOR(IRef3+2)-YP1
 545 :                YPT=XCOOR(IRef4+2)-YP1
 546 :                IF (IDIM.GE.3) THEN
 547 :                   ZP1=XCOOR(IRef1+3)
 548 :                   ZVN=XCOOR(IRef2+3)-ZP1
 549 :                   ZV=XCOOR(IRef3+3)-ZP1
 550 :                   ZPT=XCOOR(IRef4+3)-ZP1
 551 :                ELSE
 552 :                   ZP1=0.D0
 553 :                   ZVN=0.D0
 554 :                   ZV=0.D0
 555 :                   ZPT=0.D0
 556 :                ENDIF
 557 : C  Determination du vecteur directeur de l'axe du tore
 558 :                DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN)
 559 :                IF (DVN.EQ.0.D0) THEN
 559 :                   CALL ERREUR(21)
 560 :                   SEGDES MELEME
 561 :                   RETURN
 562 :                ENDIF
 563 :                XVN=XVN/DVN
 564 :                YVN=YVN/DVN
 565 :                ZVN=ZVN/DVN
 566 : C  Calcul du grand rayon du tore
 567 :                GRayon=XV*XV+YV*YV+ZV*ZV
 568 :                IF (GRayon.EQ.0.D0) THEN
 568 :                   CALL ERREUR(21)
 569 :                   SEGDES MELEME
 570 :                   RETURN
 571 :                ENDIF
 572 : C  Calcul du petit rayon du tore
 573 :                SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
 574 :                XV=XPT-SCAL*XVN
 575 :                YV=YPT-SCAL*YVN
 576 :                ZV=ZPT-SCAL*ZVN
 577 :                Val1=SQRT(Grayon/(XV*XV+YV*YV+ZV*ZV))
 578 :                XV=XV*Val1-XPT
 579 :                YV=YV*Val1-YPT
 580 :                ZV=ZV*Val1-ZPT
 581 :                PRayon=SQRT(XV*XV+YV*YV+ZV*ZV)
 582 :                IF (PRayon.EQ.0.D0) THEN
 582 :                   CALL ERREUR(21)
 583 :                   SEGDES MELEME
 584 :                   RETURN
 585 :                ENDIF
 586 : C  Recherche des points verifiant la condition d'appartenance (INDIC)
 587 :                SEGINI INDIC
 588 :                NBON=0
 589 :                DO IP=1,NBELEM
 590 :                   IPoin=NUM(1,IP)
 591 :                   IRef=(IPoin-1)*idimp1
 592 :                   XPT=XCOOR(IRef+1)-XP1
 593 :                   YPT=XCOOR(IRef+2)-YP1
 594 :                   ZPT=0.D0
 595 :                   IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1
 596 :                   SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN
 597 :                   XV=XPT-SCAL*XVN
 598 :                   YV=YPT-SCAL*YVN
 599 :                   ZV=ZPT-SCAL*ZVN
 600 :                   Val1=SQRT(GRayon/(XV*XV+YV*YV+ZV*ZV))
 601 :                   XV=XV*Val1-XPT
 602 :                   YV=YV*Val1-YPT
 603 :                   ZV=ZV*Val1-ZPT
 604 :                   Val2=SQRT(XV*XV+YV*YV+ZV*ZV)
 605 :                   IF (ABS(Val2-PRayon).LE.CRIT) THEN
 606 :                      INDIC(IP)=IPoin
 607 :                      NBON=NBON+1
 608 :                   ENDIF
 609 :                ENDDO
 610 :             ENDIF
 611 : C  Creation du maillage de sortie
 612 : C ----------------------
 613 :             IF (NBON.EQ.0) THEN
 613 :                CALL ERREUR(18)
 614 :                SEGDES MELEME
 615 :                SEGSUP INDIC
 616 :                RETURN
 617 :             ENDIF
 618 :             NBEL=NBELEM
 619 :             NBNN=1
 620 :             NBELEM=NBON
 621 :             NBSOUS=0
 622 :             NBREF=0
 623 :             SEGINI IPT1
 624 :             IPT1.ITYPEL=1
 625 :             IPLAC=0
 626 :             DO IP=1,NBEL
 627 :                IF (INDIC(IP).NE.0) THEN
 628 :                   IPLAC=IPLAC+1
 629 :                   IPT1.NUM(1,IPLAC)=INDIC(IP)
 630 :                   IPT1.ICOLOR(IPLAC)=IDCOUL
 631 :                ENDIF
 632 :             ENDDO
 633 :             SEGSUP INDIC
 634 :             SEGDES MELEME,IPT1
 634 :             CALL ECROBJ('MAILLAGE',IPT1)
 635 :          ENDIF
 636 : 
 637 : C  Operateur POINT - syntaxe 3 :
 638 : C -------------------------------
 639 : C  Extraction d'un CHAML ou d'un CHPO du(des) point(s) verifiant une
 640 : C  condition indiquee par MOT1.
 641 : C    Geo1 = POINT Che1 Mot1 ...
 642 :       ELSE
 643 :          Val1=0.D0
 644 :          Val2=0.D0
 645 :          IPList=0
 646 :          IF (IMM.GT.2) THEN
 646 :             CALL LIRREE(Val1,1,IRetou)
 647 :             IF (IERR.NE.0) RETURN
 648 :             IF (IMM.EQ.9) THEN
 648 :                CALL LIRREE(Val2,1,IRetou)
 649 :                IF (IERR.NE.0) RETURN
 650 :             ENDIF
 651 :          ENDIF
 651 :          CALL LIRMOT(MOABS,1,IAB,0)
 651 :          CALL LIRMOT(MOTAV,2,IAV,0)
 652 :          IF (IAV.EQ.0) IAV=1
 652 :          CALL LIROBJ('LISTMOTS',IPList,0,IRetou)
 653 :          IF (IERR.NE.0) RETURN
 653 :          CALL LIROBJ('CHPOINT ',IPCHP,0,IRetou)
 654 :          IF (IERR.NE.0) RETURN
 655 :          IF (IRetou.NE.0) THEN
 655 :             CALL EXPCHP(IPCHP,IMM,IAB,IAV,IPList,Val1,Val2,IPMAIL)
 656 :          ELSE
 656 :             CALL LIROBJ('MCHAML  ',IPCHE,1,IRetou)
 657 :             IF (IERR.NE.0) RETURN
 657 :             CALL EXPCHE(IPCHE,IMM,IAB,IAV,IPList,Val1,Val2,IPMAIL)
 658 :          ENDIF
 658 :          IF (IERR.EQ.0) CALL ECROBJ('MAILLAGE',IPMAIL)
 659 :       ENDIF
 660 : 
 661 :       RETURN
 662 :       END
 663 : 
 664 :  
 665 :  
 666 :  
 667 :  
 668 :  
 669 :