C POIEXT SOURCE CB215821 23/07/12 21:15:10 11704 c======================================================================= C Sert a extraire un point d'un objet MAILLAGE C C Novembre 1985 : C Extrait d'un objet tous les points appartenant a une DROIte ou a C un PLAN (possibilite de critere de distance) C C Fevrier 1986 : C Extrait d'un objet tous les points appartenant a un CYLIndre, un C CONE, une SPHEre ou un TORE C C 09/2003 : C Introduction de la DIMEnsion 1. Ajout de la syntaxe pour creer un C point a partir de ses coordonnees. Inhibation d'options non compa- C tibles avec la dimension 1. Option PROC adapte au cas IDIM=1. C Ajout de tests sur le type d'element pour POINT Entier, C POINT INITIAL et POINT FINAL (conforme a la notice). C C 01/2004 : C Modifications des methodes de recherche des points appartenant a C une DROIte, un PLAN, un CYLIndre, un CONE. Methodes moins sensibles C aux imprecisions numeriques. C C 07/2013 GOUNAND C Les options INIT FINA et PROC détectent les maillages vides C en entrée et provoquent une erreur C C 04/2015 CB215821 C Ajout de la possibilite de generer un MAILLAGE de POI1 en fournissant C des LISTREELS de coordonnees et optionnellement un LISTREEL de densite C C 11/2016 JCARDO C Ajout de l'option 'JONC' (en liaison avec les developpements dans C l'operateur PART) C c 11/2016 BP c ajout de la syntaxe : ptij = POIN ieme_noeud jeme_element maillage; c c======================================================================= SUBROUTINE POIEXT IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) -INC PPARAM -INC CCOPTIO -INC SMELEME -INC SMLREEL -INC SMLMOTS -INC SMCOORD -INC CCGEOME C Tableau fixe pour stocker les MLREEL (3 en 3D + densite) INTEGER LISLRE(4) C- DIMENSION XNORM(3) CHARACTER*4 MOT,IPRO CHARACTER*4 MCLE(10),MOTM(9),MOABS(1),MOTAV(2),MODENS(1) CHARACTER*8 MOTTYP SEGMENT INDIC(NBELEM) SEGMENT NETN(nbpts+1) C COMMON de sauvegarde lors d'un appel par SYMTRI COMMON / CSYMTR / XNORM(3) DATA MCLE / 'INIT','FINA','JONC','PROC', . 'DROI','PLAN','CYLI','CONE','SPHE','TORE' / DATA MOTM / 'MAXI','MINI','SUPE','EGSU','EGAL','EGIN','INFE', . 'DIFF','COMP' / DATA MOTAV / 'AVEC','SANS' / DATA MOABS / 'ABS ' / DATA MODENS /'DENS'/ segact mcoord c Lecture des mots cles correspondant a la 3eme syntaxe (CHPOINT) CALL LIRMOT(MOTM,9,IMM,0) IF (IMM.NE.0) GOTO 300 c syntaxe 3 ==> GOTO 300 idimp1=IDIM+1 c Lecture des mots cles geometrique CALL LIRMOT(MCLE(5),6,IMC,0) IF (IMC.NE.0) GOTO 250 c syntaxe 2.b ==> GOTO 250 c Lecture d'un maillage en entree CALL LIROBJ('MAILLAGE',MELEME,0,IRetou) IF (IERR.NE.0) RETURN IF (IRetou.NE.0) GOTO 200 c syntaxe 2.a ==> GOTO 200 100 CONTINUE c IF(IIMPI.EQ.333) WRITE(IOIMP,*) 'POINT syntaxe 1:' C ------------------------------- C Operateur POINT - syntaxe 1 : C ------------------------------- C CB215821 Ajout de la possibilité de générer un MAILLAGE de POI1 en C fournissant directement un LISTREEL par COORDONNEES plus C optionnellement un LISTREEL pour la densité CALL LIROBJ('LISTREEL',MLREEL,0,IRETOU) IF (IRETOU.NE.0) THEN IF(IDIM .EQ. 0) THEN C Dimension non définie CALL ERREUR(14) RETURN ENDIF SEGACT,MLREEL LISLRE(1)=MLREEL NBVAL = PROG(/1) C Lecture OBLIGATOIRE d'un nombre de LISTREEL egal a IDIM DO I=2,IDIM CALL LIROBJ('LISTREEL',MLREEL,1,IRETOU) IF (IERR .NE. 0) RETURN SEGACT,MLREEL JG=MLREEL.PROG(/1) C Verification de la longeur des LISTREELS lus IF (JG .NE. NBVAL) THEN CALL ERREUR(217) RETURN ENDIF LISLRE(I)=MLREEL ENDDO C Lecture mot-cle 'DENS' CALL LIRMOT(MODENS,1,IPOS,0) IF (IERR.NE.0) RETURN IF(IPOS .EQ. 1)THEN CALL LIROBJ('LISTREEL',MLREEL,1,IRETOU) IF (IERR.NE.0) RETURN SEGACT,MLREEL JG=MLREEL.PROG(/1) C Verification de la longeur du LISTREELS lu IF (JG .NE. NBVAL) THEN CALL ERREUR(217) RETURN ENDIF LISLRE(IDIM+1)=MLREEL ELSE LISLRE(IDIM+1)=0 ENDIF C On remplit le MELEME de POI1 NBNN = 1 NBELEM = NBVAL NBSOUS = 0 NBREF = 0 SEGINI,MELEME ITYPEL = 1 IF (NBVAL .GT. 0) THEN C Ajout des NOEUDS au SEGMENT MCOORD et ajustement SEGACT,MCOORD*MOD NBANC=nbpts NBPTS=NBANC+NBVAL SEGADJ,MCOORD NBAD = 0 IDINI = NBANC*(IDIM+1) DO I=1,NBVAL DO J=1,IDIM+1 C Le cas (IDIM+1) sert à remplir la densité NBAD = NBAD+1 MLREEL=LISLRE(J) IF(J .EQ. IDIM+1)THEN IF (MLREEL .EQ. 0) XVAL=DENSIT ELSE XVAL=MLREEL.PROG(I) ENDIF XCOOR(IDINI+NBAD)=XVAL ENDDO NUM(1,I) =NBANC+I ICOLOR(I)=IDCOUL ENDDO ENDIF SEGDES,MCOORD CALL ACTOBJ('MAILLAGE',MELEME,1) CALL ECROBJ('MAILLAGE',MELEME) RETURN ELSE C Creation d'un seul POINT par trois flottants C Pour IDIM = 1, seul moyen de creer un point. CALL MESLIR(-149) CALL LIRREE(Val1,1,IRetou) IF (IERR.NE.0) RETURN CALL MESLIR(-150) IF (IDIM.EQ.2) THEN CALL LIRREE(Val2,1,IRetou) ELSE IF (IDIM.EQ.3) THEN CALL LIRREE(Val2,1,IRetou) IF (IERR.NE.0) RETURN CALL MESLIR(-151) CALL LIRREE(Val3,1,IRetou) ENDIF IF (IERR.NE.0) RETURN C Lecture mot-cle 'DENS' CALL LIRMOT(MODENS,1,IPOS,0) IF (IERR.NE.0) RETURN IF(IPOS .EQ. 1)THEN CALL MESLIR(-238) CALL LIRREE(ValDen,1,IRetou) IF (IERR.NE.0) RETURN ELSE ValDen=DENSIT ENDIF SEGACT MCOORD*MOD NbPts=1+nbpts SEGADJ MCOORD IRef=(NbPts-1)*idimp1 XCOOR(IRef+1)=Val1 IF (IDIM.EQ.2) THEN XCOOR(IRef+2)=Val2 ELSE IF (IDIM.EQ.3) THEN XCOOR(IRef+2)=Val2 XCOOR(IRef+3)=Val3 ENDIF ENDIF XCOOR(IRef+idimp1)=ValDen SEGDES,MCOORD CALL ECROBJ('POINT ',NbPts) RETURN 200 CONTINUE c IF(IIMPI.EQ.333) WRITE(IOIMP,*) 'POINT syntaxe 2:' C ------------------------------- C Operateur POINT - syntaxe 2 : C ------------------------------- C 1) Recuperation du POINT1 du maillage GEO1 par son numero C Point1 = POINT Geo1 Entier1 ; C Note : le maillage est constitue uniquement de POI1 ou SEG2. C 1') Recuperation de POINT1 = Entier1^eme noeud du Entier2^eme element de GEO1 C Point1 = POINT Geo1 Entier1 Entier2; C Note : le maillage est constitue uniquement de POI1 ou SEG2. C 2) Recuperation du point INITIAL ou FINAL du maillage GEO1 C Point1 = POINT Geo1 'INITIAL' | 'FINAL' ; C Note : le maillage est constitue uniquement de POI1, SEG2 ou SEG3. C 3) Recuperation du POINT1 du maillage GEO1 le plus proche du POINT2. C Point1 = POINT Geo1 'PROC' Point2 ; CALL ACTOBJ('MAILLAGE',MELEME,1) C Lecture facultative d'un entier (numero du point a trouver) CALL LIRENT(JPTR,0,IRetou) c cas ou on a bien lu un entier IF (IRetou.NE.0) THEN IF (JPTR.LE.0) THEN CALL ERREUR(262) RETURN ENDIF MOT=' ' IPRO=' ' C Lecture d'un 2eme entier CALL LIRENT(KPTR,0,IRetou) IF (IRetou.NE.0) THEN IF (KPTR.LE.0) THEN CALL ERREUR(262) RETURN ENDIF MOT='MELE' ENDIF ELSE C Lecture d'un des mots cles ('INIT','FINA','JONC','PROC') CALL LIRCHA(MOT,0,IRetou) IF (IRetou.EQ.0) THEN C Le type de l'operande est incorrect CALL QUETYP(MOTERR(1:8),0,IRetou) IF (IRetou.NE.0) THEN CALL ERREUR (39) ELSE CALL ERREUR(533) ENDIF RETURN ENDIF C point INIT ou FINA => IPRO=' ' IF ((MOT.EQ.MCLE(1)).OR.(MOT.EQ.MCLE(2))) THEN IPRO=' ' C point JONC ou PROCH ELSEIF (MOT.EQ.MCLE(3)) THEN IPRO=MOT GOTO 203 ELSEIF (MOT.EQ.MCLE(4)) THEN IPRO=MOT GOTO 204 ELSE MOTERR=MOT CALL ERREUR(7) RETURN ENDIF ENDIF C Traitement des syntaxes POINT Geo1 Entier1 | 'INIT' | 'FINA' C Erreur si objet maillage complexe IF (LISOUS(/1).NE.0) THEN CALL ERREUR(25) RETURN ENDIF IF (NUM(/1).LE.0.OR.NUM(/2).LE.0) THEN MOTERR='MAILLAGE' C Une donnée de type %m1:8 est vide CALL ERREUR(1027) RETURN ENDIF c --- syntaxe POINT Geo1 Entier1 Entier2 --- IF (MOT.EQ.'MELE') THEN c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option Entier1 Entier2' IPTR=JPTR JPTR=KPTR IF (IPTR.GT.NUM(/1).OR.JPTR.GT.NUM(/2)) THEN CALL ERREUR(262) RETURN ENDIF c --- syntaxe POINT Geo1 Entier1 --- ELSEIF (MOT.EQ.' ') THEN c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option Entier1' IF ((ITYPEL.NE.1).AND.(JPTR.EQ.(NUM(/2)+1))) THEN IPTR=NUM(/1) JPTR=NUM(/2) ELSE IPTR=1 ENDIF IF (JPTR.GT.NUM(/2)) THEN CALL ERREUR(262) RETURN ENDIF c --- syntaxe POINT Geo1 'INIT' --- ELSE IF (MOT.EQ.MCLE(1)) THEN c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option INIT' IPTR=1 JPTR=1 c --- syntaxe POINT Geo1 'FINA' --- c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option FINA' ELSE IF (MOT.EQ.MCLE(2)) THEN IPTR=NUM(/1) JPTR=NUM(/2) ENDIF c ---partie commune aux syntaxes Entier(s) INIT FINA --- IPoin=NUM(IPTR,JPTR) CALL ECROBJ('POINT ',IPoin) RETURN C --- Traitement de l'option POIN 'JONC' --- 203 CONTINUE c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option JONC' NBSOUS=LISOUS(/1) NBS=MAX(1,NBSOUS) IPT1=MELEME NBELEM=0 SEGINI,NETN DO IO=1,NBS IF (NBSOUS.GT.0) THEN IPT1=LISOUS(IO) ENDIF DO 222 J=1,IPT1.NUM(/2) DO 223 I=1,IPT1.NUM(/1) IF (NETN(IPT1.NUM(I,J)).EQ.2) NBELEM=NBELEM+1 NETN(IPT1.NUM(I,J))=NETN(IPT1.NUM(I,J))+1 223 CONTINUE 222 CONTINUE ENDDO NBNN=1 NBREF=0 NBSOUS=0 SEGINI,IPT1 IPT1.ITYPEL=1 K=0 DO I=1,NETN(/1) IF (NETN(I).GT.2) THEN K=K+1 IPT1.NUM(1,K)=I ENDIF ENDDO SEGSUP,NETN CALL ECROBJ('MAILLAGE',IPT1) RETURN C --- Traitement de l'option POIN 'PROC' --- 204 CONTINUE c IF(IIMPI.EQ.333) WRITE(IOIMP,*) '+ option PROCH' CALL LIROBJ('POINT ',IPoin1,1,IRetou) IF (IERR.NE.0) THEN RETURN ENDIF Val1=1.E30 IRef=(IPoin1-1)*idimp1 XP1=XCOOR(IRef+1) YP1=XCOOR(IRef+2) ZP1=0.D0 IF (IDIM.GE.3) ZP1=XCOOR(IRef+3) IPT1=MELEME Ipoin=0 DO IOB=1,MAX(1,LISOUS(/1)) IF (LISOUS(/1).NE.0) THEN IPT1=LISOUS(IOB) ENDIF DO i=1,IPT1.NUM(/2) DO j=1,IPT1.NUM(/1) IRef=(IPT1.NUM(j,i)-1)*idimp1 IF (IDIM.EQ.1) THEN DVN=ABS(XP1-XCOOR(IRef+1)) ELSE XVN=XP1-XCOOR(IRef+1) YVN=YP1-XCOOR(IRef+2) ZVN=0.D0 IF (IDIM.GE.3) ZVN=ZP1-XCOOR(IRef+3) DVN=XVN*XVN+YVN*YVN+ZVN*ZVN ENDIF IF (DVN.LE.Val1) THEN Val1=DVN IPoin=IPT1.NUM(j,i) ENDIF ENDDO ENDDO ENDDO IF (IPoin.EQ.0) THEN MOTERR='MAILLAGE' C Une donnée de type %m1:8 est vide CALL ERREUR(1027) ELSE IF (IERR.EQ.0) CALL ECROBJ('POINT ',IPoin) ENDIF RETURN 250 CONTINUE c IF(IIMPI.EQ.333) WRITE(IOIMP,*) 'POINT syntaxe 2.b :' C ------------------------------- C Operateur POINT - syntaxe 2.b : C ------------------------------- C Extraction du(des) point(s) du maillage GEO1 situe(s) sur une C DROIte, un PLAN, un CYLIndre, un CONE, une SPHEre ou un TORE. C L'appartenance du point a l'objet est definie suivant un critere C de distance non nul (appartenance stricte : critere tres faible). C Geof = POINT Geo1 MOT1 Points ... (Critere) C En dimension 1, ces criteres n'ont aucun sens. IF (IDIM.EQ.1) THEN MOTERR=MCLE(3+IMC) INTERR(1)=IDIM CALL ERREUR(971) RETURN ENDIF C Lectures obligatoire du maillage et facultative du critere CALL LIROBJ('MAILLAGE',MELEME,1,IRetou) IF (IERR.NE.0) RETURN CRIT=DENSIT/10.D0 CALL LIRREE(CRIT,0,IRetou) IF (CRIT.LE.0.D0) CALL ERREUR(17) IF (IERR.NE.0) RETURN C Transformation du maillage en 'POI1' CALL CHANGE(MELEME,1) NBELEM=NUM(/2) C Intersection du maillage avec une DROIte C -------------------------- IF (IMC.EQ.1) THEN CALL LIROBJ('POINT ',IPoin1,1,IRetou) CALL LIROBJ('POINT ',IPoin2,1,IRetou) IF (IERR.NE.0) THEN RETURN ENDIF C Recherche du vecteur directeur de la DROIte IRef1=(IPoin1-1)*idimp1 IRef2=(IPoin2-1)*idimp1 XP1=XCOOR(IRef1+1) XVN=XCOOR(IRef2+1)-XP1 YP1=XCOOR(IRef1+2) YVN=XCOOR(IRef2+2)-YP1 ZP1=0.D0 ZVN=0.D0 IF (IDIM.GE.3) THEN ZP1=XCOOR(IRef1+3) ZVN=XCOOR(IRef2+3)-ZP1 ENDIF DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN) IF (DVN.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF XVN=XVN/DVN YVN=YVN/DVN ZVN=ZVN/DVN XNORM(1)=XVN XNORM(2)=YVN XNORM(3)=ZVN C Recherche des points verifiant la condition d'appartenance (INDIC) CRIT=CRIT*CRIT SEGINI INDIC NBON=0 DO IP=1,NBELEM IPoin=NUM(1,IP) IRef=(IPoin-1)*idimp1 XPT=XCOOR(IRef+1)-XP1 YPT=XCOOR(IRef+2)-YP1 ZPT=0.D0 IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1 SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN XV=SCAL*XVN-XPT YV=SCAL*YVN-YPT ZV=SCAL*ZVN-ZPT Val1=XV*XV+YV*YV+ZV*ZV IF (Val1.LE.CRIT) THEN INDIC(IP)=IPoin NBON=NBON+1 ENDIF ENDDO C Intersection du maillage avec un PLAN C -------------------------- ELSE IF (IMC.EQ.2) THEN CALL LIROBJ('POINT ',IPoin1,1,IRetou) CALL LIROBJ('POINT ',IPoin2,1,IRetou) CALL LIROBJ('POINT ',IPoin3,1,IRetou) IF (IERR.NE.0) THEN RETURN ENDIF C Determination du vecteur normal au PLAN IRef1=(IPoin1-1)*idimp1 IRef2=(IPoin2-1)*idimp1 IRef3=(IPoin3-1)*idimp1 XP1=XCOOR(IRef1+1) XP2=XCOOR(IRef2+1) XP3=XCOOR(IRef3+1) YP1=XCOOR(IRef1+2) YP2=XCOOR(IRef2+2) YP3=XCOOR(IRef3+2) ZP1=0.D0 ZP2=0.D0 ZP3=0.D0 IF (IDIM.GE.3) THEN ZP1=XCOOR(IRef1+3) ZP2=XCOOR(IRef2+3) ZP3=XCOOR(IRef3+3) ENDIF XVN=(YP2-YP1)*(ZP3-ZP1)-(ZP2-ZP1)*(YP3-YP1) YVN=(ZP2-ZP1)*(XP3-XP1)-(XP2-XP1)*(ZP3-ZP1) ZVN=(XP2-XP1)*(YP3-YP1)-(YP2-YP1)*(XP3-XP1) DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN) IF (DVN.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF XVN=XVN/DVN YVN=YVN/DVN ZVN=ZVN/DVN XNORM(1)=XVN XNORM(2)=YVN XNORM(3)=ZVN C Recherche des points verifiant la condition d'appartenance (INDIC) SEGINI INDIC NBON=0 DO IP=1,NBELEM IPoin=NUM(1,IP) IRef=(IPoin-1)*idimp1 XV=XCOOR(IRef+1)-XP1 YV=XCOOR(IRef+2)-YP1 ZV=0.D0 IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1 Val1=ABS(XV*XVN+YV*YVN+ZV*ZVN) IF (Val1.LE.CRIT) THEN INDIC(IP)=IPoin NBON=NBON+1 ENDIF ENDDO C Intersection du maillage avec un CYLIndre C -------------------------- ELSE IF (IMC.EQ.3) THEN CALL LIROBJ('POINT ',IPoin1,1,IRetou) CALL LIROBJ('POINT ',IPoin2,1,IRetou) CALL LIROBJ('POINT ',IPoin3,1,IRetou) IF (IERR.NE.0) THEN RETURN ENDIF IRef1=(IPoin1-1)*idimp1 IRef2=(IPoin2-1)*idimp1 IRef3=(IPoin3-1)*idimp1 XP1=XCOOR(IRef1+1) XVN=XCOOR(IRef2+1)-XP1 XPT=XCOOR(IRef3+1)-XP1 YP1=XCOOR(IRef1+2) YVN=XCOOR(IRef2+2)-YP1 YPT=XCOOR(IRef3+2)-YP1 IF (IDIM.GE.3) THEN ZP1=XCOOR(IRef1+3) ZVN=XCOOR(IRef2+3)-ZP1 ZPT=XCOOR(IRef3+3)-ZP1 ELSE ZP1=0.D0 ZVN=0.D0 ZPT=0.D0 ENDIF C Determination de l'axe du cylindre (axe P1-P2) DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN) IF (DVN.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF XVN=XVN/DVN YVN=YVN/DVN ZVN=ZVN/DVN C Calcul du rayon du cylindre (distance P3 a axe P1-P2) SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN XV=SCAL*XVN-XPT YV=SCAL*YVN-YPT ZV=SCAL*ZVN-ZPT Rayon=SQRT(XV*XV+YV*YV+ZV*ZV) IF (Rayon.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF C Recherche des points verifiant la condition d'appartenance (INDIC) SEGINI INDIC NBON=0 DO IP=1,NBELEM IPoin=NUM(1,IP) IRef=(IPoin-1)*idimp1 XPT=XCOOR(IRef+1)-XP1 YPT=XCOOR(IRef+2)-YP1 ZPT=0.D0 IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1 SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN XV=SCAL*XVN-XPT YV=SCAL*YVN-YPT ZV=SCAL*ZVN-ZPT Val1=SQRT(XV*XV+YV*YV+ZV*ZV) IF (ABS(Val1-Rayon).LE.CRIT) THEN INDIC(IP)=IPoin NBON=NBON+1 ENDIF ENDDO C Intersection du maillage avec un CONE C -------------------------- ELSE IF (IMC.EQ.4) THEN CALL LIROBJ('POINT ',IPoin1,1,IRetou) CALL LIROBJ('POINT ',IPoin2,1,IRetou) CALL LIROBJ('POINT ',IPoin3,1,IRetou) IF (IERR.NE.0) THEN RETURN ENDIF IRef1=(IPoin1-1)*idimp1 IRef2=(IPoin2-1)*idimp1 IRef3=(IPoin3-1)*idimp1 XP1=XCOOR(IRef1+1) XVN=XCOOR(IRef2+1)-XP1 XV=XCOOR(IRef3+1)-XP1 YP1=XCOOR(IRef1+2) YVN=XCOOR(IRef2+2)-YP1 YV=XCOOR(IRef3+2)-YP1 IF (IDIM.GE.3) THEN ZP1=XCOOR(IRef1+3) ZVN=XCOOR(IRef2+3)-ZP1 ZV=XCOOR(IRef3+3)-ZP1 ELSE ZP1=0.D0 ZVN=0.D0 ZV=0.D0 ENDIF C Recherche du vecteur directeur de l'axe du cone DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN) IF (DVN.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF XVN=XVN/DVN YVN=YVN/DVN ZVN=ZVN/DVN C Calcul de l'angle au sommet du cone Val1=XV*XV+YV*YV+ZV*ZV IF (Val1.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF SCAL=XV*XVN+YV*YVN+ZV*ZVN XPT=SCAL*XVN-XV YPT=SCAL*YVN-YV ZPT=SCAL*ZVN-ZV SISOM=SQRT((XPT*XPT+YPT*YPT+ZPT*ZPT)/Val1) C Recherche des points verifiant la condition d'appartenance (INDIC) SEGINI INDIC NBON=0 DO IP=1,NBELEM IPoin=NUM(1,IP) IRef=(IPoin-1)*idimp1 XV=XCOOR(IRef+1)-XP1 YV=XCOOR(IRef+2)-YP1 ZV=0.D0 IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1 Val1=SQRT(XV*XV+YV*YV+ZV*ZV)*SISOM SCAL=XV*XVN+YV*YVN+ZV*ZVN XPT=SCAL*XVN-XV YPT=SCAL*YVN-YV ZPT=SCAL*ZVN-ZV Val2=SQRT(XPT*XPT+YPT*YPT+ZPT*ZPT) IF (ABS(Val2-Val1).LE.CRIT) THEN INDIC(IP)=IPoin NBON=NBON+1 ENDIF ENDDO C Intersection du maillage avec une SPHERE C -------------------------- ELSE IF (IMC.EQ.5) THEN CALL LIROBJ('POINT ',IPoin1,1,IRetou) CALL LIROBJ('POINT ',IPoin2,1,IRetou) IF (IERR.NE.0) THEN RETURN ENDIF IRef1=(IPoin1-1)*idimp1 IRef2=(IPoin2-1)*idimp1 XP1=XCOOR(IRef1+1) XV=XCOOR(IRef2+1)-XP1 YP1=XCOOR(IRef1+2) YV=XCOOR(IRef2+2)-YP1 IF (IDIM.GE.3) THEN ZP1=XCOOR(IRef1+3) ZV=XCOOR(IRef2+3)-ZP1 ELSE ZP1=0.D0 ZV=0.D0 ENDIF C Calcul du rayon de la sphere Rayon=SQRT(XV*XV+YV*YV+ZV*ZV) IF (Rayon.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF C Recherche des points verifiant la condition d'appartenance (INDIC) SEGINI INDIC NBON=0 DO IP=1,NBELEM IPoin=NUM(1,IP) IRef=(IPoin-1)*idimp1 XV=XCOOR(IRef+1)-XP1 YV=XCOOR(IRef+2)-YP1 ZV=0.D0 IF (IDIM.GE.3) ZV=XCOOR(IRef+3)-ZP1 Val1=SQRT(XV*XV+YV*YV+ZV*ZV) IF (ABS(Val1-Rayon).LE.CRIT) THEN INDIC(IP)=IPoin NBON=NBON+1 ENDIF ENDDO C Intersection du maillage avec un TORE C -------------------------- ELSE IF (IMC.EQ.6) THEN CALL LIROBJ('POINT ',IPoin1,1,IRetou) CALL LIROBJ('POINT ',IPoin2,1,IRetou) CALL LIROBJ('POINT ',IPoin3,1,IRetou) CALL LIROBJ('POINT ',IPoin4,1,IRetou) IF (IERR.NE.0) THEN RETURN ENDIF IRef1=(IPoin1-1)*idimp1 IRef2=(IPoin2-1)*idimp1 IRef3=(IPoin3-1)*idimp1 IRef4=(IPoin4-1)*idimp1 XP1=XCOOR(IRef1+1) XVN=XCOOR(IRef2+1)-XP1 XV=XCOOR(IRef3+1)-XP1 XPT=XCOOR(IRef4+1)-XP1 YP1=XCOOR(IRef1+2) YVN=XCOOR(IRef2+2)-YP1 YV=XCOOR(IRef3+2)-YP1 YPT=XCOOR(IRef4+2)-YP1 IF (IDIM.GE.3) THEN ZP1=XCOOR(IRef1+3) ZVN=XCOOR(IRef2+3)-ZP1 ZV=XCOOR(IRef3+3)-ZP1 ZPT=XCOOR(IRef4+3)-ZP1 ELSE ZP1=0.D0 ZVN=0.D0 ZV=0.D0 ZPT=0.D0 ENDIF C Determination du vecteur directeur de l'axe du tore DVN=SQRT(XVN*XVN+YVN*YVN+ZVN*ZVN) IF (DVN.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF XVN=XVN/DVN YVN=YVN/DVN ZVN=ZVN/DVN C Calcul du grand rayon du tore GRayon=XV*XV+YV*YV+ZV*ZV IF (GRayon.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF C Calcul du petit rayon du tore SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN XV=XPT-SCAL*XVN YV=YPT-SCAL*YVN ZV=ZPT-SCAL*ZVN Val1=SQRT(Grayon/(XV*XV+YV*YV+ZV*ZV)) XV=XV*Val1-XPT YV=YV*Val1-YPT ZV=ZV*Val1-ZPT PRayon=SQRT(XV*XV+YV*YV+ZV*ZV) IF (PRayon.EQ.0.D0) THEN CALL ERREUR(21) RETURN ENDIF C Recherche des points verifiant la condition d'appartenance (INDIC) SEGINI INDIC NBON=0 DO IP=1,NBELEM IPoin=NUM(1,IP) IRef=(IPoin-1)*idimp1 XPT=XCOOR(IRef+1)-XP1 YPT=XCOOR(IRef+2)-YP1 ZPT=0.D0 IF (IDIM.GE.3) ZPT=XCOOR(IRef+3)-ZP1 SCAL=XPT*XVN+YPT*YVN+ZPT*ZVN XV=XPT-SCAL*XVN YV=YPT-SCAL*YVN ZV=ZPT-SCAL*ZVN Val1=SQRT(GRayon/(XV*XV+YV*YV+ZV*ZV)) XV=XV*Val1-XPT YV=YV*Val1-YPT ZV=ZV*Val1-ZPT Val2=SQRT(XV*XV+YV*YV+ZV*ZV) IF (ABS(Val2-PRayon).LE.CRIT) THEN INDIC(IP)=IPoin NBON=NBON+1 ENDIF ENDDO ENDIF C Creation du maillage de sortie C ---------------------- IF (NBON.EQ.0) THEN CALL ERREUR(18) SEGSUP INDIC RETURN ENDIF NBEL=NBELEM NBNN=1 NBELEM=NBON NBSOUS=0 NBREF=0 SEGINI IPT1 IPT1.ITYPEL=1 IPLAC=0 DO IP=1,NBEL IF (INDIC(IP).NE.0) THEN IPLAC=IPLAC+1 IPT1.NUM(1,IPLAC)=INDIC(IP) IPT1.ICOLOR(IPLAC)=IDCOUL ENDIF ENDDO SEGSUP INDIC CALL ACTOBJ('MAILLAGE',IPT1,1) CALL ECROBJ('MAILLAGE',IPT1) RETURN 300 CONTINUE c IF(IIMPI.EQ.333) WRITE(IOIMP,*) 'POINT syntaxe 3 :' C ------------------------------- C Operateur POINT - syntaxe 3 : C ------------------------------- C Extraction d'un CHAML ou d'un CHPO du(des) point(s) verifiant une C condition indiquee par MOT1. C Geo1 = POINT Che1 Mot1 ... Val1=0.D0 Val2=0.D0 IPList=0 IF (IMM.GT.2) THEN CALL LIRREE(Val1,1,IRetou) IF (IERR.NE.0) RETURN IF (IMM.EQ.9) THEN CALL LIRREE(Val2,1,IRetou) IF (IERR.NE.0) RETURN ENDIF ENDIF CALL LIRMOT(MOABS,1,IAB,0) CALL LIRMOT(MOTAV,2,IAV,0) IF (IAV.EQ.0) IAV=1 CALL LIROBJ('LISTMOTS',IPList,0,IRetou) IF (IERR.NE.0) RETURN IF(IRETOU .EQ. 1) THEN MLMOTS=IPList SEGACT,MLMOTS ENDIF CALL QUETYP(MOTTYP,1,IRetou) CALL LIROBJ(MOTTYP,IPOIN,1,IRetou) IF (IERR.NE.0) RETURN CALL ACTOBJ(MOTTYP,IPOIN,1) IF (IERR.NE.0) RETURN IF (MOTTYP .EQ. 'CHPOINT ') THEN CALL EXPCHP(IPOIN,IMM,IAB,IAV,IPList,Val1,Val2,IPMAIL) ELSEIF (MOTTYP .EQ. 'MCHAML ') THEN CALL EXPCHE(IPOIN,IMM,IAB,IAV,IPList,Val1,Val2,IPMAIL) ELSE CALL ERREUR(686) RETURN ENDIF IF (IERR.NE.0) RETURN CALL ACTOBJ('MAILLAGE',IPMAIL,1) CALL ECROBJ('MAILLAGE',IPMAIL) END