Sources Esope | Cast3M

sens1
C SENS1     SOURCE    CB215821  23/01/25    21:15:34     11573          
      SUBROUTINE SENS1
c
c    sous routine castem 2000
c
c    utilisée pour l'opérateur SENS
c    determine dans quel sens tourne un contour en 2 D
c
c   synthaxe de l'opérateur   TAB2 = SENS TAB1 ;
c
c   tab1 est issue de CCON et contient des maillages indicé par
c        des entiers. ces maillages doivent etre des contours fermes
c        orientes.
c
c   tab2 contient des entiers (+/-1) indice par les meme entiers
c       +1 indique que le contour est dans le sens trigonometrique
c       -1 indique que le contour est dans le sens inverse
c
c  langage Fortran 77 + esope
c  auteur A de Gayffier
c
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      CHARACTER CHARI,CHARR
      LOGICAL LOGII,LOGIR
c
-INC CCREEL
-INC SMTABLE
-INC SMELEME
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
c
c   declaration du segment INCR
c
      SEGMENT INCR
        INTEGER NUM0(NBELEM)
      ENDSEGMENT
c
      IOBIN=0
c   operateur disponible seulement den dimension 2
      IF  (IDIM .NE. 2 ) THEN
         CALL ERREUR(19)
         RETURN
      ENDIF
c
      CALL LIROBJ ('TABLE',MTABLE,1,IRET)
      IF (IRET .EQ. 0) RETURN
      SEGACT ,MTABLE
c
c     compte combien d'indice entier
      NBENT = 0
      DO 5 I=1,MTABTI(/2)
         IF (MTABTI(I) .EQ. 'ENTIER' ) NBENT = NBENT +1
 5    CONTINUE
c
      M = 0
      SEGINI  ,MTAB1
c
c    boucle sur tous les contours
c
      SEGACT,MCOORD
      DO 10 INDTAB=1,NBENT
c
         CALL ACCTAB(MTABLE,'ENTIER',INDTAB,XVALI,CHARI,LOGII,0,
     &        'MAILLAGE',IVALR,XVALR,CHARR,LOGIR,IRETR)
         IF (IERR  .NE. 0 ) RETURN
c
c
         MELEME = IRETR
         SEGACT MELEME
c
c         quelques controles sur le maillage
c
         IF  (LISOUS(/1) .NE. 0 .OR. (ITYPEL .NE. 2
     &        .AND. ITYPEL .NE. 3)) THEN
            CALL ERREUR(26)
            SEGDES MELEME
            SEGSUP MTAB1
            RETURN
         ENDIF
c
c   initialisation de la boucle sur les segments2
c
         NBELEM = NUM(/2)
         SEGINI INCR
         I = 1
         IF (ITYPEL .EQ. 2 ) THEN
            I1 = NUM(1,1)
            I2 = NUM(2,1)
         ELSE IF  ( ITYPEL .EQ. 3 ) THEN
            I1 = NUM(1,1)
            I2 = NUM(3,1)
         ENDIF
         STHETA = 0.D0
         ISEG2  = 0
c
c   Initialisation du tableau NUM0
c
         DO 9 J=1,NBELEM
            NUM0(J) = NUM(1,J)
 9       CONTINUE
c
c    debut de la boucle sur les segments seg2 du maillage
c
 15      CONTINUE
c
         DO 20 J=1,NBELEM
            IF ( NUM0(J) .EQ. I2 ) THEN
c            on a trouve le segment suivant
               ISEG2 = J
               NUM0(J) = 0 ;
               IF (ITYPEL .EQ. 2 ) THEN
                  I3 = NUM(2,J)
               ELSE IF  ( ITYPEL .EQ. 3 ) THEN
                  I3 = NUM(3,J)
               ENDIF
               GOTO 30
            ENDIF
 20      CONTINUE
c
c      on n'est pas parvenu a boucler le contour
c
c
c    dans le cas axisymmetrique le contour n'est pas necessairement fermé
c
         IF (IFOUR .EQ. 0) THEN
c
c           on verifie que le point est sur l'axe OZ
c
            XI1 = ABS(XCOOR((I2-1)*(IDIM+1) +1 ))
            DENS = MAX(XCOOR((I2-1)*(IDIM+1)),1.D-10)
c       print *,'i2=',i2,'r=',xi1,'dens=',dens
            IF ( XI1 .GT. (DENS/100.d0)) GOTO 25
c
c           on cherche l'autre point sur l'axe
c
            DO 23 J=1,NBELEM
c
               XI2 =ABS(XCOOR((NUM(1,J)-1)*(IDIM+1) +1 ))
               DENS =MAX(XCOOR((NUM(1,J)-1)*(IDIM+1)),1.D-10)
c
               IF (XI2.LT.(DENS/100.d0)) THEN
                  I3 = NUM(1,J)
*       print *, 'on a trouvé le second point',i3,'r=',xi2
                  GOTO 30
               ENDIF
 23         CONTINUE
*        print *, 'pas trouvé le second point sur l axe'
c
         ENDIF
C
 25      CONTINUE
c
c    gestion de l'erreur contour non ferme
c
         CALL ERREUR(28)
         SEGSUP MTAB1
         SEGDES MELEME,MTABLE
         RETURN
c
 30      CONTINUE
c
c       I sert a compter le nombre de segments pris en compte
c
         I = I + 1
c
c    calcul de l'angle entre les deux segments
c
         XI1 = XCOOR((I1-1)*(IDIM+1) +1 )
         YI1 = XCOOR((I1-1)*(IDIM+1) +2 )
         XI2 = XCOOR((I2-1)*(IDIM+1) +1 )
         YI2 = XCOOR((I2-1)*(IDIM+1) +2 )
         XI3 = XCOOR((I3-1)*(IDIM+1) +1 )
         YI3 = XCOOR((I3-1)*(IDIM+1) +2 )
         X12 = XI2-XI1
         Y12 = YI2-YI1
         X23 = XI3-XI2
         Y23 = YI3-YI2
         D12 = SQRT( X12*X12+Y12*Y12)
         D23 = SQRT ( X23*X23+Y23*Y23)
         SITHET = ( X12*Y23 - X23*Y12 ) / D12 / D23
         COTHET = ( X12*X23 + Y12*Y23 ) / D12 / D23
         IF ( (ABS(SITHET)) .LT. 1.D-10 ) THEN
            THETA = 0.D0
         ELSE
            IF ( ABS(SITHET) .GT. 1.D0 ) THEN
               THETA = SIGN((XPI/2.D0),SITHET)
            ELSE
               THETA = ASIN( SITHET )
            ENDIF
         ENDIF
         IF (COTHET.LT.0.D0 .AND. SITHET.GT.0.D0) THETA=XPI-THETA
         IF (COTHET.LT.0.D0 .AND. SITHET.LT.0.D0) THETA=-XPI-THETA
         STHETA = STHETA + THETA
         I1 = I2
         I2 = I3
c
c       recherche du point suivant si on n'est pas parvenu
c       au point initial
c
         IF ( ISEG2 .NE. 1 ) GOTO 15
c
c        on controle que tous les segments ont ete examinés
c
         IF (I.NE.(NBELEM+1) .AND.  I.NE.(NBELEM+2)) THEN
c       nbelem+2 correspond au cas axis avec un contour non fermé
            CALL ERREUR(28)
            SEGSUP MTAB1
            SEGDES MELEME,MTABLE
            RETURN
         ENDIF
c
c      determination du sens du contour la somme vaut +/- XPI
c
         IF ( ABS(STHETA - XPI*2.D0 ) .LT. 1.D-3 ) THEN
            IVALR = 1
         ELSE IF ( ABS(STHETA + XPI*2.D0 ) .LT. 1.D-3 ) THEN
            IVALR = -1
         ELSE
            IVALR = 0
c          il y un problème: c'est sur!
c           print *, 'somme des angles en radian',stheta
            INTERR(1)=INDTAB
            CALL ERREUR(718)
         ENDIF
c
c    ecriture dans la table du resultat
c
         CALL ECCTAB(MTAB1,'ENTIER',INDTAB,XVALI,CHARI,LOGII,IOBIN,
     &        'ENTIER',IVALR,XVALR,CHARR,LOGIR,IRETR)
         SEGDES MELEME
         SEGSUP INCR
c
 10   CONTINUE
      SEGDES,MCOORD
c
c        il n'y a plus rien dans la table
      CALL ECROBJ('TABLE',MTAB1)
      SEGDES MTABLE,MTAB1
      RETURN
c
      END