extrac.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE extrac
!                       *****************
!
     &(x,y,som,ikle,indic,nelem,nelmax,npoin,nsom,projec)
!
!***********************************************************************
! PROGICIEL : STBTEL V5.2    07/12/88    J-M HERVOUET (LNH) 30 87 80 18
!                            19/02/93    J-M JANIN    (LNH) 30 87 72 84
!                                        A   WATRIN
!***********************************************************************
!
!  FONCTION  :  PREPARATION DE DONNEES AVANT L'APPEL DE FMTSEL
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________
! !      NOM       !MODE!                   ROLE
! !________________!____!______________________________________________
! !________________!____!______________________________________________
! MODE : -->(DONNEE NON MODIFIEE), <--(RESULTAT), <-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
!
! APPELE PAR : PREDON
! APPEL DE : -
!
!***********************************************************************
!
      USE declarations_special
      IMPLICIT NONE
!
      INTEGER, INTENT(INOUT) :: NELEM
      INTEGER, INTENT(IN) :: NELMAX,NPOIN,NSOM
      INTEGER, INTENT(INOUT) :: IKLE(nelmax,3),INDIC(npoin)
      DOUBLE PRECISION, INTENT(INOUT) :: X(npoin),Y(npoin)
      DOUBLE PRECISION, INTENT(IN) :: SOM(10,2)
      LOGICAL, INTENT(IN) :: PROJEC
!
      LOGICAL FLAG,F1,F2,F3
!
      INTEGER IELEM,IPOIN,ISOM,IDP,I1,I2,I3
      DOUBLE PRECISION DX,DY,A1,A2,A3
!
!=======================================================================
! BOUCLE SUR TOUS LES PLANS DE COUPE
!=======================================================================
!
      DO isom = 1,nsom
!
        dx = som(isom+1,1) - som(isom,1)
        dy = som(isom+1,2) - som(isom,2)
!
!=======================================================================
! POUR UN DEMI PLAN DE COUPE DONNE :
!      RECHERCHE DES POINTS EXT.(=0) , INT.(=1) , SUR LE BORD (=2)
!=======================================================================
!
        DO ipoin = 1,npoin
          indic(ipoin) = 0
          IF (dx*(y(ipoin)-som(isom,2)).GE.dy*(x(ipoin)-som(isom,1)))
     &        indic(ipoin) = 1
        ENDDO
!
        ielem = 1
20      CONTINUE
        i1 = indic(ikle(ielem,1))
        i2 = indic(ikle(ielem,2))
        i3 = indic(ikle(ielem,3))
        IF (i1.EQ.0.OR.i2.EQ.0.OR.i3.EQ.0) THEN
          IF (i1.EQ.1) indic(ikle(ielem,1)) = 2
          IF (i2.EQ.1) indic(ikle(ielem,2)) = 2
          IF (i3.EQ.1) indic(ikle(ielem,3)) = 2
          ikle(ielem,1) = ikle(nelem,1)
          ikle(ielem,2) = ikle(nelem,2)
          ikle(ielem,3) = ikle(nelem,3)
          nelem = nelem - 1
        ELSE
          ielem = ielem + 1
        ENDIF
        IF (ielem.NE.nelem+1) GOTO 20
!
!=======================================================================
! POUR UN DEMI PLAN DE COUPE DONNE :
!      ELIMINATION DES ELEMENTS DEGENERES
!=======================================================================
!
30      CONTINUE
        ielem = 1
        flag = .false.
35      CONTINUE
        i1 = ikle(ielem,1)
        i2 = ikle(ielem,2)
        i3 = ikle(ielem,3)
        f1 = indic(i1).EQ.2
        f2 = indic(i2).EQ.2
        f3 = indic(i3).EQ.2
        IF (f1.AND.f2.AND.f3) THEN
          ikle(ielem,1) = ikle(nelem,1)
          ikle(ielem,2) = ikle(nelem,2)
          ikle(ielem,3) = ikle(nelem,3)
          nelem = nelem - 1
        ELSE
          IF (f1.AND.f2) THEN
            IF (dx*(x(i2)-x(i1))+dy*(y(i2)-y(i1)).LE.0.d0) THEN
              flag = .true.
              indic(i3) = 2
            ENDIF
          ENDIF
          IF (f2.AND.f3) THEN
            IF (dx*(x(i3)-x(i2))+dy*(y(i3)-y(i2)).LE.0.d0) THEN
              flag = .true.
              indic(i1) = 2
            ENDIF
          ENDIF
          IF (f3.AND.f1) THEN
            IF (dx*(x(i1)-x(i3))+dy*(y(i1)-y(i3)).LE.0.d0) THEN
              flag = .true.
              indic(i2) = 2
            ENDIF
          ENDIF
          ielem = ielem + 1
        ENDIF
        IF (ielem.NE.nelem+1) GOTO 35
        IF (flag) GOTO 30
!
!=======================================================================
! POUR UN DEMI PLAN DE COUPE DONNE :
!      PROJECTION DES NOUVEAUX POINTS DE BORD
!=======================================================================
!
        IF (projec) THEN
          a1 = 1.d0 / (dx*dx + dy*dy)
          a2 = a1 * (som(isom,1)*som(isom+1,2) -
     &               som(isom,2)*som(isom+1,1) )
          DO idp = 1,3
            DO ielem = 1,nelem
              ipoin = ikle(ielem,idp)
              IF (indic(ipoin).EQ.2) THEN
                a3 = a1*(x(ipoin)*dx+y(ipoin)*dy)
                x(ipoin) = dx*a3 + dy*a2
                y(ipoin) = dy*a3 - dx*a2
              ENDIF
            ENDDO
          ENDDO
        ENDIF
!
      ENDDO !ISOM
!
!=======================================================================
!
      RETURN
      END