contri.f

作者：Scott Sloan 相关文献参考： Adv Eng Software 1987 V9No1,34-

*              adjacent triangles
*     T      - Adjacency array for triangulation
*            - Triangles adjacent to J are found in T(I,J) for I=1,2,3
*              J=1,2,...,NTRI
*            - Adjacent triangles listed in anticlockwise sequence
*            - Zero denotes no adjacent triangle
*     TN     - List of triangle numbers such that vertex I can be 
*              found in triangle TN(I)
*     W      - Not defined, used as workspace
*
*     OUTPUT:
*     ------- 
*
*     VI,VJ  - Unchanged
*     NPTS   - Unchanged
*     NTRI   - Unchanged
*     NEF    - If VI-VJ needs to be forced, NEF is incremented by unity
*            - Else NEF is unchanged
*     JW     - Unchanged
*     X      - Unchanged
*     Y      - Unchanged
*     V      - Vertex array for triangulation updated so that edge 
*              V1-V2 is present
*     T      - Adjacency array for triangulation updated so that edge 
*              V1-V2 is present
*     TN     - List of triangle numbers updated so that edge 
*              V1-V2 is present
*     W      - List of new edges that replace old edges in
*              triangulation
*            - Vertices in W(1,I) and W(2,I) define each new edge I
*
*     NOTES:
*     ------
*
*     - This routine assumes that the edge defined by VI-VJ does not
*       lie on an outer boundary of the triangulation and, thus, the
*       triangulation must include the triangles that are formed with
*       the supertriangle vertices
*
*     PROGRAMMER:             Scott Sloan
*     -----------
*
*     LAST MODIFIED:          3 march 1991          Scott Sloan
*     --------------
*
*     COPYRIGHT 1990:         Scott Sloan
*     ---------------         Department of Civil Engineering
*                             University of Newcastle
*                             NSW 2308
*
************************************************************************
      INTEGER A,C,E,I,L,R,S
      INTEGER JW,NC,V1,V2,V3,V4,VI,VJ
      INTEGER ELR,ERL,NEF
      INTEGER LAST,NTRI,NPTS
      INTEGER FIRST
      INTEGER T(3,NTRI),V(3,NTRI),W(2,JW)
      INTEGER TN(NPTS+3)
*
      DOUBLE PRECISION X1,X2,X3,X4,XI,XJ,Y1,Y2,Y3,Y4,YI,YJ
      DOUBLE PRECISION TOL,X13,X14,X23,X24,Y13,Y14,Y24,Y23
      DOUBLE PRECISION COSA,COSB
      DOUBLE PRECISION C00000,DETIJ3
      DOUBLE PRECISION X(NPTS+3),Y(NPTS+3)
* 
      LOGICAL SWAP
*
      PARAMETER (TOL=1.d-6)
      PARAMETER (C00000=0.d0)
*---------------------------------------------------------------------
*     Check data
*
      IF(VI.LE.0.OR.VI.GT.NPTS)THEN
        nerror = 2
        error(1) = 'Illegal value for VI in subroutine EDGE:'
        write(error(2),'(''VI='',i5)')VI
        call erroroutput
        STOP
      ENDIF
      IF(VJ.LE.0.OR.VJ.GT.NPTS)THEN
        nerror = 2
        error(1) = 'Illegal value for VJ in subroutine EDGE:'
        write(error(2),'(''VJ='',i5)')VJ
        call erroroutput
        STOP
      ENDIF
*---------------------------------------------------------------------- 
*     Find any triangle which has VI as a vertex
*
      S=TN(VI)
      IF(S.LE.0)THEN
        nerror = 1
        write(error(1),1000)VI
        call erroroutput
        STOP
      ENDIF           
      IF(TN(VJ).LE.0)THEN
        nerror = 1
        write(error(1),1000)VJ
        call erroroutput
        STOP
      ENDIF
      XI=X(VI)
      YI=Y(VI)
      XJ=X(VJ)
      YJ=Y(VJ)
*----------------------------------------------------------------------
*     Find an arc that crosses VI-VJ
*     C=current triangle
*     S=triangle in which search is started
*
      C=S
*
*     Vertices V1 and V2 are such that V1-V2 is opposite VI
*     Circle anticlockwise round VI until V1-V2 crosses VI-VJ
*
   10 CONTINUE
      IF(V(1,C).EQ.VI)THEN
        V2=V(3,C)
        E =2
      ELSEIF(V(2,C).EQ.VI)THEN
        V2=V(1,C)
        E =3
      ELSE   
        V2=V(2,C)
        E =1
      ENDIF
*
*     Test if arc VI-VJ already exists
*
      IF(V2.EQ.VJ)RETURN
*
*     Test if V1-V2 crosses VI-VJ
*
      X2=X(V2)
      Y2=Y(V2)
      IF((XI-X2)*(YJ-Y2).GT.(XJ-X2)*(YI-Y2))THEN
*
*       V2 is left of VI-VJ
*       Check if V1 is right of VI-VJ
*
        V1=V(E,C)
        X1=X(V1)
        Y1=Y(V1)
        IF((XI-X1)*(YJ-Y1).LT.(XJ-X1)*(YI-Y1))THEN
*
*         V1-V2 crosses VI-VJ , so edge needs to be forced
*
          NEF=NEF+1
          GOTO 15
        ENDIF
      ENDIF
*
*     No crossing, move anticlockwise around VI to the next triangle
*
      C=T(MOD(E,3)+1,C)
      IF(C.NE.S)GOTO 10
      nerror = 3
      write(error(1),'(''Vertex adjacent to'',i5,'' is on arc'')')VI
      write(error(2),'(''between vertices'',i5,'' and'',i5)')VI,VJ
      error(3) = 'in subroutine EDGE'
      call erroroutput
      STOP
   15 CONTINUE
*-------------------------------------------------------------------
*     Loop to store all arcs which cross VI-VJ
*     Vertices V1/V2 are right/left of   VI-VJ 
*
      NC=0
   20 CONTINUE
      NC=NC+1
      IF(NC.GT.JW)THEN
        nerror = 1
        error(1) ='Not enough workspace in subroutine EDGE, increase JW'
        call erroroutput
        STOP
      ENDIF
      W(1,NC)=V1
      W(2,NC)=V2
      C=T(E,C)
      IF(V(1,C).EQ.V2)THEN
        V3=V(3,C)
        E =2
      ELSEIF(V(2,C).EQ.V2)THEN
        V3=V(1,C)
        E =3
      ELSE
        V3=V(2,C)
        E =1
      ENDIF
*
*     Termination test, all arcs crossing VI-VJ have been stored
*
      IF(V3.EQ.VJ)GOTO 30
      X3=X(V3)
      Y3=Y(V3)
      DETIJ3=(XI-X3)*(YJ-Y3)-(XJ-X3)*(YI-Y3)
      IF(DETIJ3.LT.C00000)THEN
         E =MOD(E,3)+1
         V1=V3
      ELSEIF(DETIJ3.GT.C00000)THEN
         V2=V3
      ELSE 
        nerror = 3
        write(error(1),'(''Vertex'',i5,'' is on arc'')')V3
        write(error(2),'(''between vertices'',i5,'' and'',i5)')VI,VJ
        error(3) = 'in subroutine EDGE'
        call erroroutput
        STOP
      ENDIF
      GOTO 20
   30 CONTINUE
*-------------------------------------------------------------------
*     Swap each arc that crosses VI-VJ if it is a diagonal of a 
*     convex quadrilateral
*     Execute all possible swaps, even if newly formed arc also 
*     crosses VI-VJ, and iterate until no arcs cross
*
      LAST=NC
   35 IF(LAST.GT.0)THEN
        FIRST=1
   40   IF(FIRST.LE.LAST)THEN
*
*         Find triangle L which is left of V1-V2 
*         Find triangle R which is right of V1-V2
*
          V1=W(1,FIRST)
          V2=W(2,FIRST)
*
*         Exchange V1 and V2 if V1 is a supertriangle vertex
*
          IF(V1.GT.NPTS)THEN
            IF(V2.GT.NPTS)THEN
               nerror = 4
               write(error(1),
     &            '(''Arc between vertices'',i5,'' and'',i5)')VI,VJ
               error(2) = 'crosses supertriangle boundary defined by'
               write(error(3),'(''vertices'',i5,''and'',i5)')V1,V2
               error(4) = 'in subroutine EDGE'
               call erroroutput
               STOP
            ENDIF
            W(1,FIRST)=V2
            W(2,FIRST)=V1
            V2=V1
            V1=W(1,FIRST)
          ENDIF
          L=TN(V1)
   45     IF(V(1,L).EQ.V1)THEN
            IF(V(2,L).EQ.V2)THEN
              V3 =V(3,L)
              ELR=1
              R  =T(1,L)
            ELSE
              L=T(3,L)
              GOTO 45
            ENDIF
          ELSEIF(V(2,L).EQ.V1)THEN
            IF(V(3,L).EQ.V2)THEN
              V3 =V(1,L)
              ELR=2
              R  =T(2,L)
            ELSE
              L=T(1,L)
              GOTO 45
            ENDIF
          ELSEIF(V(1,L).EQ.V2)THEN
              V3 =V(2,L)
              ELR=3
              R  =T(3,L)
          ELSE
              L=T(2,L)
              GOTO 45
          ENDIF
*
*         Find vertices V3 and V4 where:
*         triangle L is defined by V3-V1-V2 
*         triangle R is defined by V4-V2-V1
*
          IF(T(1,R).EQ.L)THEN
            V4=V(3,R)
            A =T(2,R)
            ERL=1
          ELSEIF(T(2,R).EQ.L)THEN
            V4=V(1,R)
            A =T(3,R)
            ERL=2
          ELSE
            V4=V(2,R)
            A =T(1,R)
            ERL=3
          ENDIF
*
*         Test if quad formed by V3-V1-V4-V2 is convex
*
          X3=X(V3)
          Y3=Y(V3)
          X4=X(V4)
          Y4=Y(V4)
          X1=X(V1)
          Y1=Y(V1)
          X2=X(V2)
          Y2=Y(V2)
          IF((X3-X1)*(Y4-Y1).LT.(X4-X1)*(Y3-Y1))THEN
            IF((X3-X2)*(Y4-Y2).GT.(X4-X2)*(Y3-Y2))THEN
*
*             Quad is convex so swap diagonal arcs
*             Update vertex and adjacency lists for triangle L
*
              IF(ELR.EQ.1)THEN
                V(2,L)=V4
                C     =T(2,L)
                T(1,L)=A
                T(2,L)=R
              ELSEIF(ELR.EQ.2)THEN
                V(3,L)=V4
                C     =T(3,L)
                T(2,L)=A
                T(3,L)=R             
              ELSE
                V(1,L)=V4
                C     =T(1,L)
                T(3,L)=A
                T(1,L)=R             
              ENDIF                 
*
*             Update vertex and adjacency lists for triangle R
*
              IF(ERL.EQ.1)THEN
                V(2,R)=V3
                T(1,R)=C
                T(2,R)=L
              ELSEIF(ERL.EQ.2)THEN
                V(3,R)=V3
                T(2,R)=C
                T(3,R)=L             
              ELSE
                V(1,R)=V3
                T(3,R)=C
                T(1,R)=L             
              ENDIF
*
*             Update adjacency lists for triangles A and C
*
              IF(T(1,C).EQ.L)THEN
                T(1,C)=R
              ELSEIF(T(2,C).EQ.L)THEN
                T(2,C)=R
              ELSE
                T(3,C)=R
              ENDIF
              IF(T(1,A).EQ.R)THEN
                T(1,A)=L
              ELSEIF(T(2,A).EQ.R)THEN
                T(2,A)=L
              ELSE
                T(3,A)=L
              ENDIF
*
*             Update vertex-triangle list
*
              TN(V1)=L
              TN(V2)=R
*
*             Test if new diagonal arc crosses VI-VJ and store it if it
*             does
*
              IF(((XI-X3)*(YJ-Y3)-(XJ-X3)*(YI-Y3))*
     +          ((XI-X4)*(YJ-Y4)-(XJ-X4)*(YI-Y4)).LT.C00000)THEN
                W(1,FIRST)=V4
                W(2,FIRST)=V3
                FIRST=FIRST+1
              ELSE
                W(1,FIRST)=W(1,LAST)
                W(1,LAST) =V3
                W(2,FIRST)=W(2,LAST)
                W(2,LAST) =V4        
                LAST=LAST-1                    
              ENDIF
              GOTO 40
            ENDIF
          ENDIF
*
*         Arc cannot be swapped, so move to next intersecting arc
*
          FIRST=FIRST+1
          GOTO 40
        ENDIF
        GOTO 35
      ENDIF
*----------------------------------------------------------------------
*     Optimise all new arcs (except VI-VJ)
*        
      SWAP=.TRUE.
   50 IF(SWAP)THEN
        SWAP=.FALSE.
        DO 70 I=2,NC
*
*         Find triangle L which is left of V1-V2 
*         Find triangle R which is right of V1-V2
*
          V1=W(1,I)
          V2=W(2,I)
*
*         Exchange V1 and V2 if V1 is a supertriangle vertex
*
          IF(V1.GT.NPTS)THEN
            IF(V2.GT.NPTS)THEN
               nerror = 3
               write(error(1),
     &            '(''Arc between vertices'',i5,'' and'',i5)')V1,V2
               error(2) = 
     &            'cannot be optimised since it is a supertriangle'
               error(3) = 'boundary in subroutine EDGE'
               call erroroutput
               STOP
            ENDIF
            W(1,I)=V2