Sources Esope | Cast3M

arc
C ARC       SOURCE    CB215821  23/01/25    21:15:05     11573          
C     1234567890124567890124567898012456789012345678901234567890123456
 
C     SOURCE   L. DI VALENTIN    20/06/97
 
      SUBROUTINE ARC(POINT,CENTRE,NORMAL,ANGLE,INBR,
     &     DEN1OX,DEN2OX,POINT2)
 
 
C     Cette routine, appele par prtran.eso, est chargee de construire
C     un arc de cercle en accord avec les elements utilises, de centre
C     CENTRE, de normale NORMAL, d'extremites POINT et POINT2 et
C     d'angle ANGLE.
C     Il y a une redondance entre ces 3 derniers parametres. En fait,
C     l'ANGLE est prioritaire devant POINT2 a moins que POINT2 = POINT
C     (cf plus loin). La procedure suppose que POINT2 est la 2eme
C     extremite de l'arc de cercle. Il est donc important que POINT2
C     reference ce point geometrique.
 
C     INBR est le nombre de decoupage de l'arc de cercle.
C     Si INBR = 0, les densites DEN1OX et DEN2OX sont prises en
C     compte pour la decoupe.
C
C     Si POINT2 = POINT alors quelque soit l'angle, la subroutine
C     construira un cercle. Dans ce cas, ANGLE sera change egal a
C     2*pi (=(2.D0*XPI))
C
C     Le maillage resultant est sauve dans la pile a la fin.
 
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
-INC CCREEL
 
-INC PPARAM
-INC CCOPTIO
-INC CCGEOME
-INC SMELEME
-INC SMCOORD
 
 
      SEGMENT TABPAR(NBELEM)
 
      INTEGER POINT,CENTRE,NORMAL,INBR,POINT2
      REAL*8    ALPHA,ANGLE,DEN1,DEN2
 
 
      DEN1 = DEN1OX
      DEN2 = DEN2OX
 
C     Vecteur Normal :
 
      IF (IDIM.EQ.3) THEN
         Xn = XCOOR((NORMAL-1)*(IDIM+1) + 1)
         Yn = XCOOR((NORMAL-1)*(IDIM+1) + 2)
         Zn = XCOOR((NORMAL-1)*(IDIM+1) + 3)
 
         ALPHA = SQRT(Xn*Xn+Yn*Yn+Zn*Zn)
         IF(ALPHA.LT.XPETIT) ALPHA = 1.D0
         Xn = Xn / ALPHA
         Yn = Yn / ALPHA
         Zn = Zn / ALPHA
 
      ELSE
         Xn = 0.D0
         Yn = 0.D0
         Zn = 1.D0
 
      ENDIF
 
 
C     On recupere les coordonnees du CENTRE
 
      XCEN = XCOOR((CENTRE-1)*(IDIM+1) + 1)
      YCEN = XCOOR((CENTRE-1)*(IDIM+1) + 2)
      IF (IDIM.EQ.3) THEN
         ZCEN = XCOOR((CENTRE-1)*(IDIM+1) + 3)
      ENDIF
 
C     On recupere les coordonnees du vecteur (CENTRE -> POINT)
 
      X1 = XCOOR((POINT-1)*(IDIM+1) + 1) - XCEN
      Y1 = XCOOR((POINT-1)*(IDIM+1) + 2) - YCEN
 
      IF (IDIM.EQ.3) THEN
         Z1 = XCOOR((POINT-1)*(IDIM+1) + 3) - ZCEN
      ELSE
         Z1 = 0.D0
      ENDIF
 
      RAYON = SQRT(X1*X1+Y1*Y1+Z1*Z1)
 
 
C     On calcule le vecteur definissant avec le vecteur precedent
C     le 2eme vecteur de la base du cercle
 
      X2 = Z1*Yn - Y1*Zn
      Y2 = X1*Zn - Z1*Xn
      Z2 = Y1*Xn - X1*Yn
 
 
C     Cas du cercle :
 
      IF (POINT.EQ.POINT2) THEN
         ANGLE =  (2.D0*XPI)
      ENDIF
 
 
      DLONG = ABS(ANGLE * RAYON)
      IF(DLONG.LT.XPETIT) DLONG = 1.D0
C      DEN1  = 1.D0/INBR
      DEN1  = DEN1/DLONG
      DEN2  = DEN2/DLONG
 
      CALL DECOUP(INBR,DEN1,DEN2,APROG,NBELEM,DENI,DECA,DLONG)
      DLONG = ANGLE * RAYON
 
C      NBELEM = INBR
      NX     = NBELEM - 1
 
      NBSOUS=0
      NBREF=0
      IF (ILCOUR.EQ.0) GOTO 666
      IF (KDEGRE(ILCOUR).EQ.0) GOTO 666
      NBNN=NBNNE(KDEGRE(ILCOUR))
      IF (NBNN.EQ.0) GOTO 666
 
 
      SEGINI MELEME
      SEGINI TABPAR
 
      ITYPEL = KDEGRE(ILCOUR)
 
      DIN      = DEN1
 
      IPOO     = nbpts
      NUM(1,1) = POINT
 
 
      IF (NX.EQ.0) GOTO 110
 
      DO 102 I=1,NX
         DIN = DIN * APROG
         IF (ITYPEL.EQ.2) GOTO 103
         IPOO=IPOO+1
         NUM(2,I)=IPOO
         NUM(3,I)=IPOO+1
 103     IPOO=IPOO+1
         TABPAR(I) = DIN
         NUM(1,I+1)=IPOO
         IF (ITYPEL.EQ.3) GOTO 102
         NUM(2,I)=IPOO
 
 102  CONTINUE
 
 
 110  DIN = DIN * APROG
      IF (ITYPEL.NE.3) GOTO 115
      IPOO=IPOO+1
      NUM(2,NBELEM)=IPOO
      TABPAR(NBELEM) = DIN
 115  CONTINUE
 
 
      NUM(ITYPEL,NBELEM)=POINT2
 
      SEGACT MCOORD*mod
      IADR=NBPTS
 
 
C     Creation des points
 
      DPAR=0
 
      IF (NX.EQ.0) GOTO 300
      NBPTS=IADR+NX*(ITYPEL-1)
 
 
      SEGADJ MCOORD
 
      DO 301 I=1,NX
 
         DIN = TABPAR(I)
         DPAR=DPAR+DIN
         IF (ITYPEL.EQ.2) GOTO 1103
         DPAR1=DPAR-DIN*0.5D0
 
         PARA=ANGLE*DPAR1
         COSA=COS(PARA)
         SINA=SIN(PARA)
 
         XCOOR(IADR*(IDIM+1)+1)=XCEN+COSA*X1+SINA*X2
         XCOOR(IADR*(IDIM+1)+2)=YCEN+COSA*Y1+SINA*Y2
         IF (IDIM.EQ.3) THEN
            XCOOR(IADR*(IDIM+1)+3)=ZCEN+COSA*Z1+SINA*Z2
         ENDIF
 
         XCOOR((IADR+1)*(IDIM+1))=DENI+DECA*DPAR1
         IADR=IADR+1
 
 1103    PARA=ANGLE*DPAR
         COSA=COS(PARA)
         SINA=SIN(PARA)
 
         XCOOR(IADR*(IDIM+1)+1)=XCEN+COSA*X1+SINA*X2
         XCOOR(IADR*(IDIM+1)+2)=YCEN+COSA*Y1+SINA*Y2
         IF (IDIM.EQ.3) THEN
            XCOOR(IADR*(IDIM+1)+3)=ZCEN+COSA*Z1+SINA*Z2
         ENDIF
 
         XCOOR((IADR+1)*(IDIM+1))=DENI+DECA*DPAR
         IADR=IADR+1
 
 301  CONTINUE
 
      NBPTS = IADR
 
 
      SEGADJ MCOORD
 
 
 300  DIN = TABPAR(NBELEM)
      DPAR=DPAR+DIN
 
      IF (ITYPEL.NE.3) GOTO 303
 
      NBPTS = IADR+1
      SEGADJ MCOORD
      DPAR1 = DPAR-DIN*0.5D0
 
      PARA = ANGLE * DPAR1
      COSA = COS(PARA)
      SINA = SIN(PARA)
 
      XCOOR(IADR*(IDIM+1)+1) = XCEN + X1*COSA + X2*SINA
      XCOOR(IADR*(IDIM+1)+2) = YCEN + Y1*COSA + Y2*SINA
      IF (IDIM.EQ.3) THEN
         XCOOR(IADR*(IDIM+1)+3) = ZCEN + Z1*COSA + Z2*SINA
      ENDIF
 
      XCOOR((IADR+1)*(IDIM+1)) = DENI+DECA*DPAR1
      IADR=IADR + 1
 
 
 303  DO 309 I=1,NUM(/2)
         ICOLOR(I)=IDCOUL
 309  CONTINUE
      SEGSUP TABPAR
      SEGDES MELEME
 
      CALL ECROBJ('MAILLAGE',MELEME)
      RETURN
* Erreur 16 Type d'element incorrect
 666  CONTINUE
      CALL ERREUR(16) 
      RETURN
      END