surcon.f

Go to the documentation of this file.

!                       *****************
                        SUBROUTINE surcon
!                       *****************
!
     &(x,y,ikle,ipo,nbor,nptfr,ncolor,ifabor,color)
!
!***********************************************************************
! PROGICIEL : STBTEL V5.2                   19/04/91  J-C GALLAND  (LNH)
!                                           19/02/93  J-M JANIN    (LNH)
!***********************************************************************
!
! FONCTION : DECOUPAGE DES TRIANGLES SURCONTRAINTS
!
!-----------------------------------------------------------------------
!                             ARGUMENTS
! .________________.____.______________________________________________.
! |      NOM       |MODE|                   ROLE                       |
! |________________|____|______________________________________________|
! |   X,Y          |<-->| COORDONNEES DU MAILLAGE .
! |   IKLE         |<-->| NUMEROS GLOBAUX DES NOEUDS DE CHAQUE ELEMENT
! |   TRAV1        |<-->| TABLEAU DE TRAVAIL
! |   NBOR         | -->| TABLEAU DES POINTS DE BORD
! |   NPTFR        | -->| NOMBRE DE POINT FRONTIERE
! |   NCOLOR       |<-->| TABLEAU DES COULEURS DES POINTS
! |   IFABOR       |<-->| TABLEAU DES VOISINS DES ELEMENTS
! |    COLOR       |<-->| STOCKAGE COULEURS DES NOEUDS SUR FICHIER GEO
! |________________|____|______________________________________________
! | COMMON:        |    |
! |  GEO:          |    |
! |    MESH        | -->| TYPE DES ELEMENTS DU MAILLAGE
! |    NDP         | -->| NOMBRE DE NOEUDS PAR ELEMENTS
! |    NPOIN       |<-->| NOMBRE TOTAL DE NOEUDS DU MAILLAGE
! |    NELEM       |<-->| NOMBRE TOTAL D'ELEMENTS DU MAILLAGE
! |    NPMAX       | -->| DIMENSION EFFECTIVE DES TABLEAUX X ET Y
! |                |    | (NPMAX = NPOIN + 0.1*NELEM)
! |    NELMAX      | -->| DIMENSION EFFECTIVE DES TABLEAUX CONCERNANT
! |                |    | LES ELEMENTS (NELMAX = NELEM + 0.2*NELEM)
! |                |    |
! |________________|____|______________________________________________|
! MODE : -->(DONNEE NON MODIFIEE), <--(RESULTAT), <-->(DONNEE MODIFIEE)
!-----------------------------------------------------------------------
! APPELE PAR : STBTEL
! APPEL DE : DECOUP
!***********************************************************************
!
      USE declarations_special
      USE declarations_stbtel, ONLY: nelem,npoin,nelmax,npmax
      USE interface_stbtel, ex_surcon => surcon
      IMPLICIT NONE
!
      DOUBLE PRECISION, INTENT(INOUT) :: X(*) , Y(*)
      INTEGER, INTENT(INOUT) :: NBOR(*) , IKLE(nelmax,4) , NCOLOR(*)
      INTEGER, INTENT(INOUT) :: IFABOR(nelmax,*) , IPO(*)
      LOGICAL, INTENT(INOUT) :: COLOR
      INTEGER, INTENT(IN) :: NPTFR
!
      INTEGER NPOIN2 , NELEM2
      INTEGER ITEST , IELEM , KELEM , ISWAP , KSWAP
      INTEGER I , J , K
      INTEGER :: IP(3) , KP , ISUI(3)
!
      parameter( isui = (/ 2 , 3 , 1 /) )
!
!=======================================================================
!
      WRITE(lu,4050)
!
      npoin2 = npoin
      nelem2 = nelem
      itest  = 0
!
!=======================================================================
! RECHERCHE DES TRIANGLES SURCONTRAINTS
!=======================================================================
!
      DO k = 1 , npmax
        ipo(k) = 0
      ENDDO
!
      DO k = 1 , nptfr
        ipo(nbor(k)) = 1
      ENDDO

      DO k = 1 , nelem
        IF (ipo(ikle(k,1))+ipo(ikle(k,2))+ipo(ikle(k,3)).EQ.3) THEN
!
! LE TRIANGLE EST SURCONTRAINT, ON LE DECOUPE EN TROIS
!
          itest  = itest + 1
          IF (itest.GT.int(0.1*nelmax)) THEN
            WRITE(lu,4000) int(0.1*nelmax)
            CALL plante(1)
            stop
          ENDIF
!
          WRITE(lu,4070) x(ikle(k,1)),y(ikle(k,1)),
     &       x(ikle(k,2)),y(ikle(k,2)),x(ikle(k,3)),y(ikle(k,3))
!
          CALL decoup (k,x,y,ikle,ncolor,ifabor,
     &                 nelem2,npoin2,color)
!
        ENDIF
      ENDDO !K
!
!=======================================================================
! REMISE A JOUR DE NPOIN ET NELEM
!=======================================================================
!
      npoin = npoin2
      nelem = nelem2
!
!=======================================================================
! MAINTENANT ON SWAPE LES ARETES CONTENANT 2 POINTS DE BORD ET QUI NE
! SONT PAS DES ARETES DE BORD
!
! REMARQUE : LES IFABOR NE SONT PAS MODIFIES DANS CE QUI SUIT CAR
!            ILS NE SERVENT PLUS APRES POUR LE LOGICIEL
!=======================================================================
!
      iswap = 0
      kswap = 0
      DO ielem = 1,nelem
        ip(1) = ikle(ielem,1)
        ip(2) = ikle(ielem,2)
        ip(3) = ikle(ielem,3)
        DO i = 1,3
          j = isui(i)
          IF(ifabor(ielem,i).GT.0.AND.ipo(ip(i))+ipo(ip(j)).EQ.2) THEN
            kelem = ifabor(ielem,i)
!
! COMPTE TENU DU PREMIER DECOUPAGE, ON EST SUR QUE TOUT ELEMENT CONTIENT
! AU MOINS UN POINT INTERIEUR, ET DANS CE CAS UN SEUL.
! K EST SUR D'ETRE REMPLI.
!
            IF (ipo(ikle(kelem,1)).EQ.0) k=1
            IF (ipo(ikle(kelem,2)).EQ.0) k=2
            IF (ipo(ikle(kelem,3)).EQ.0) k=3
            kp = ikle(kelem,k)
!
            IF ((x(kp)-x(ip(i)))*(y(ip(isui(j)))-y(ip(i))).GT.
     &          (y(kp)-y(ip(i)))*(x(ip(isui(j)))-x(ip(i))).AND.
     &          (x(kp)-x(ip(j)))*(y(ip(isui(j)))-y(ip(j))).LT.
     &          (y(kp)-y(ip(j)))*(x(ip(isui(j)))-x(ip(j)))) THEN
!
              iswap = iswap + 1
              ikle(ielem,i) = kp
              ikle(kelem,isui(k)) = ip(isui(j))
!
              WRITE(lu,4080) x(ip(i)),y(ip(i)),
     &                       x(ip(j)),y(ip(j))
!
            ELSE
!
              kswap = kswap + 1
!
            ENDIF
          ENDIF
        ENDDO
      ENDDO
!
!=======================================================================
!
      WRITE(lu,4100) itest,npoin,nelem,iswap
!
      IF (kswap.NE.0) THEN
        WRITE(lu,4200) kswap
      ENDIF
!
 4000 FORMAT(//,1x,'**************************************************',
     &        /,1x,'THE MAXIMUM NUMBER OF OVERSTRESSED TRIANGLES   ',
     &        /,1x,'IS',i5,' : IT IS GREATER IN YOUR CASE ||',
     &        /,1x,'**************************************************')
 4050 FORMAT(//,1x,'OVERSTRESSED ELEMENTS ARE CANCELLED',/,
     &          1x,'-----------------------------------',/)
 4070 FORMAT   (1x,'ADDITIONAL NODE AT CENTRE OF TRIANGLE :',/,
     &          1x,'(',d9.3,',',d9.3,'),(',d9.3,',',d9.3,'),',
     &             '(',d9.3,',',d9.3,')')
 4080 FORMAT   (1x,'SWAP OF FACE :',/,
     &          1x,'(',d9.3,',',d9.3,'),(',d9.3,',',d9.3,')')
 4100 FORMAT (/,1x,'NUMBER OF CANCELLED ELEMENTS : ',i5,/,
     &          1x,'AFTER BEING CANCELLED :',/,
     &          1x,'          NUMBER OF POINTS      : ',i5,/,
     &          1x,'          NUMBER OF ELEMENTS    : ',i5,//,
     &          1x,'MOREOVER,',i4,' TRIANGLES HAVE BEEN SWAPPED',//)
 4200 FORMAT   (1x,'      BUT',i4,' COULD NOT BE',//)
!
      RETURN
      END