contri.f

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

        ENDIF
   40 CONTINUE
      END
************************************************************************
      SUBROUTINE DELAUN(NPTS,N,X,Y,LIST,STACK,V,T,NTRI)
      use mod_output
      implicit double precision (a-h,o-z)
************************************************************************
*
*     PURPOSE:
*     --------
*
*     Assemble Delaunay triangulation
*
*     INPUT:
*     ------
*
*     NPTS   - Total number of points in data set
*     N      - Total number of points to be triangulated (N le NPTS)
*     X      - X-coords of all points in data set
*            - X-coord of point I given by X(I)
*            - If point is in LIST, coordinate must be normalised
*              such that X=(X-XMIN)/DMAX
*            - Last three locations are used to store x-coords of
*              supertriangle vertices in subroutine delaun
*     Y      - Y-coords of all points in data set
*            - Y-coord of point I given by Y(I)
*            - If point is in LIST, coordinate must be normalised
*              such that Y=(Y-YMIN)/DMAX
*            - Last three locations are used to store y-coords of
*              supertriangle vertices in subroutine delaun
*     LIST   - List of points to be triangulated
*            - Coincident points are flagged by an error message
*            - Points are ordered such that consecutive points are 
*              close to one another in the x-y plane
*     STACK  - Not defined
*     V      - Not defined
*     T      - Not defined
*     NTRI   - Not defined
*
*     OUTPUT:
*     ------- 
*
*     NPTS   - Unchanged
*     N      - Unchanged
*     X      - Unchanged
*     Y      - Unchanged
*     LIST   - Unchanged
*     STACK  - Not defined
*     V      - Vertex array for triangulation
*            - Vertices listed in anticlockwise sequence
*            - Vertices for triangle J are found in V(I,J) for I=1,2,3
*              and J=1,2,...,NTRI
*            - First vertex is at point of contact of first and third
*              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
*     NTRI   - Number of triangles in triangulation, including those
*              formed with the supertriangle vertices
*            - NTRI = 2*N+1
*
*     NOTES:
*     ------
*
*     - This is a faster version of the code appearing in AES 1987 vol 9
*       (small subroutines and functions have been coded in-line)
*     - Also some changes in code to improve efficiency
*     - Saving in cpu-time about 60 percent for larger problems
*     - A test has been implemented to detect coincident points
*     - Triangles formed with supertriangle vertices have not been 
*       deleted
*
*     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,B,C,I,J,L,N,P,R
      INTEGER V1,V2,V3
      INTEGER NPTS,NTRI
      INTEGER TOPSTK
      INTEGER T(3,2*N+1),V(3,2*N+1)
      INTEGER LIST(N)
      INTEGER STACK(NPTS)
*
      DOUBLE PRECISION D
      DOUBLE PRECISION XP,YP
      DOUBLE PRECISION TOL,X13,X23,X1P,X2P,Y13,Y23,Y1P,Y2P
      DOUBLE PRECISION COSA,COSB
      DOUBLE PRECISION C00000,C00100
      DOUBLE PRECISION X(NPTS+3),Y(NPTS+3)
*
*     TOL is the tolerance used to detect coincident points
*     The square of the distance between any two points must be
*     greater then TOL to avoid an error message
*     This value of TOL is suitable for single precision on most
*     32-bit machines (which typically have a precision of 0.000001)
*
      PARAMETER (TOL=1.d-10)
      PARAMETER (C00000=0.d0, C00100=100.d0)
*---------------------------------------------------------------------
*     Define vertices for supertriangle
*
      V1=NPTS+1
      V2=NPTS+2
      V3=NPTS+3
*---------------------------------------------------------------------
*     Set coords of supertriangle
*
      X(V1)=-C00100
      X(V2)= C00100
      X(V3)= C00000
      Y(V1)=-C00100
      Y(V2)=-C00100
      Y(V3)= C00100
*---------------------------------------------------------------------
*     Introduce first point
*     Define vertex and adjacency lists for first 3 triangles
*
      P=LIST(1)
      V(1,1)=P
      V(2,1)=V1
      V(3,1)=V2
      T(1,1)=3
      T(2,1)=0
      T(3,1)=2
      V(1,2)=P
      V(2,2)=V2
      V(3,2)=V3
      T(1,2)=1
      T(2,2)=0
      T(3,2)=3
      V(1,3)=P
      V(2,3)=V3
      V(3,3)=V1
      T(1,3)=2
      T(2,3)=0
      T(3,3)=1
*---------------------------------------------------------------------
*     Loop over each point (except first) and construct triangles
*   
      NTRI=3
      TOPSTK=0
      DO 140 I=2,N
        P=LIST(I)
        XP=X(P)
        YP=Y(P)
*
*       Locate triangle J in which point lies
*
        J=NTRI
   10   CONTINUE
        V1=V(1,J)
        V2=V(2,J)
        IF((Y(V1)-YP)*(X(V2)-XP).GT.(X(V1)-XP)*(Y(V2)-YP))THEN
          J=T(1,J)
          GOTO 10
        ENDIF
        V3=V(3,J)
        IF((Y(V2)-YP)*(X(V3)-XP).GT.(X(V2)-XP)*(Y(V3)-YP))THEN
          J=T(2,J)
          GOTO 10
        ENDIF
        IF((Y(V3)-YP)*(X(V1)-XP).GT.(X(V3)-XP)*(Y(V1)-YP))THEN
          J=T(3,J)
          GOTO 10
        ENDIF
*
*       Check that new point is not coincident with any previous point
*
        D=(X(V1)-XP)**2.d0
        IF(D.LT.TOL)THEN
           D=D+(Y(V1)-YP)**2.d0
           IF(D.LT.TOL)THEN
             write(string,2000)V1,P
             call output
             write(string,2010)
             call output
           ENDIF
        ENDIF
        D=(X(V2)-XP)**2.d0
        IF(D.LT.TOL)THEN
           D=D+(Y(V2)-YP)**2.d0
           IF(D.LT.TOL)THEN
             write(string,2000)V2,P
             call output
             write(string,2010)
             call output
           ENDIF
        ENDIF
        D=(X(V3)-XP)**2.d0
        IF(D.LT.TOL)THEN
           D=D+(Y(V3)-YP)**2.d0
           IF(D.LT.TOL)THEN
             write(string,2000)V3,P
             call output
             write(string,2010)
             call output
           ENDIF
        ENDIF
*
*       Create new vertex and adjacency lists for triangle J
*
        A=T(1,J)
        B=T(2,J)
        C=T(3,J)
        V(1,J)=P
        V(2,J)=V1
        V(3,J)=V2
        T(1,J)=NTRI+2    
        T(2,J)=A
        T(3,J)=NTRI+1
*
*       Create new triangles
*
        NTRI=NTRI+1
        V(1,NTRI)=P
        V(2,NTRI)=V2
        V(3,NTRI)=V3
        T(1,NTRI)=J
        T(2,NTRI)=B
        T(3,NTRI)=NTRI+1
        NTRI=NTRI+1
        V(1,NTRI)=P
        V(2,NTRI)=V3
        V(3,NTRI)=V1
        T(1,NTRI)=NTRI-1
        T(2,NTRI)=C
        T(3,NTRI)=J
*
*       Put each edge of triangle J on STACK
*       Store triangles on left side of each edge
*       Update adjacency lists for triangles B and C
*
        TOPSTK=TOPSTK+2
        STACK(TOPSTK-1)=J
        IF(T(1,C).EQ.J)THEN
          T(1,C)=NTRI
        ELSE
          T(2,C)=NTRI
        ENDIF
        STACK(TOPSTK)=NTRI
        IF(B.NE.0)THEN
          IF(T(1,B).EQ.J)THEN
            T(1,B)=NTRI-1
          ELSEIF(T(2,B).EQ.J)THEN
            T(2,B)=NTRI-1
          ELSE
            T(3,B)=NTRI-1
          ENDIF
          TOPSTK=TOPSTK+1
          STACK(TOPSTK)=NTRI-1
        ENDIF
*
*       Loop while STACK is not empty
*
   60   IF(TOPSTK.GT.0)THEN
*
*         Find triangles L and R which are either side of stacked edge
*         triangle L is defined by V3-V1-V2 and is left of V1-V2
*         triangle R is defined by V4-V2-V1 and is right of V1-V2

          R=STACK(TOPSTK)
          TOPSTK=TOPSTK-1
          L=T(2,R)
*
*         Check if new point P is in circumcircle for triangle L
*
          IF(T(1,L).EQ.R)THEN
            V1=V(1,L)
            V2=V(2,L)
            V3=V(3,L)
            A=T(2,L)
            B=T(3,L)
          ELSEIF(T(2,L).EQ.R)THEN
            V1=V(2,L)
            V2=V(3,L)
            V3=V(1,L)
            A=T(3,L)
            B=T(1,L)
          ELSE
            V1=V(3,L)
            V2=V(1,L)
            V3=V(2,L)
            A=T(1,L)
            B=T(2,L)
          ENDIF
          X13=X(V1)-X(V3)
          Y13=Y(V1)-Y(V3)
          X23=X(V2)-X(V3)
          Y23=Y(V2)-Y(V3)
          X1P=X(V1)-XP
          Y1P=Y(V1)-YP
          X2P=X(V2)-XP
          Y2P=Y(V2)-YP
          COSA=X13*X23+Y13*Y23
          COSB=X2P*X1P+Y1P*Y2P
          IF(COSA.LT.C00000.OR.COSB.LT.C00000)THEN
            IF((COSA.LT.C00000.AND.COSB.LT.C00000).OR.
     +        ((X13*Y23-X23*Y13)*COSB.LT.(X1P*Y2P-X2P*Y1P)*COSA))THEN
*
*             New point is inside circumcircle for triangle L
*             Swap diagonal for convex quad formed by P-V2-V3-V1
*
              C=T(3,R)
*
*             Update vertex and adjacency list for triangle R
*
              V(3,R)=V3
              T(2,R)=A
              T(3,R)=L
* 
*             Update vertex and adjacency list for triangle L
*
              V(1,L)=P
              V(2,L)=V3
              V(3,L)=V1
              T(1,L)=R
              T(2,L)=B
              T(3,L)=C
*
*             Put edges R-A and L-B on STACK
*             Update adjacency lists for triangles A and C
*
              IF(A.NE.0)THEN
                IF(T(1,A).EQ.L)THEN
                  T(1,A)=R
                ELSEIF(T(2,A).EQ.L)THEN
                  T(2,A)=R
                ELSE
                  T(3,A)=R
                ENDIF
                TOPSTK=TOPSTK+1
                IF(TOPSTK.GT.NPTS)THEN
                  nerror = 1
                  write(error(1),1000)
                  call erroroutput
                  STOP
                ENDIF
                STACK(TOPSTK)=R
              ENDIF
              IF(B.NE.0)THEN
                TOPSTK=TOPSTK+1
                IF(TOPSTK.GT.NPTS)THEN
                  nerror = 1
                  write(error(1),1000)
                  call erroroutput
                  STOP
                ENDIF
                STACK(TOPSTK)=L
              ENDIF
              T(1,C)=L
            ENDIF
          ENDIF
          GOTO 60
        ENDIF
  140 CONTINUE
*---------------------------------------------------------------------
*     Check consistency of triangulation
*
      IF(NTRI.NE.2*N+1)THEN
        nerror = 1
        error(1) = 
     &    'Incorrect number of triangles formed in subroutine DELAUN'
        call erroroutput
        STOP
      ENDIF
*---------------------------------------------------------------------
 1000 FORMAT('Stack overflow in subroutine DELAUN')
 2000 FORMAT('* WARNING: Points',i5,' and',i5,
     &       ' are coincident in subroutine DELAUN:')
 2010 FORMAT(11x,'Delete either point from list vector.')
      END
************************************************************************
      SUBROUTINE EDGE(VI,VJ,NPTS,NTRI,NEF,JW,X,Y,V,T,TN,W)
      use mod_output
      implicit double precision (a-h,o-z)
************************************************************************
*
*     PURPOSE:
*     --------
*
*     Force edge VI-VJ to be present in Delaunay triangulation
*
*     INPUT:
*     ------
*
*     VI,VJ  - Vertices defining edge to be present in triangulation
*     NPTS   - Total number of points in data set
*     NTRI   - Number of triangles in triangulation
*     NEF    - Running total of edges that have been forced
*            - Set to zero before first call to EDGE
*     JW     - Number of cols in workspace vector W 
*            - JW must not be less than the number of edges in the
*              existing triangulation that intersect VI-VJ
*     X      - X-coords of all points in data set
*            - X-coord of point I given by X(I)
*            - Last three locations are used to store x-coords of
*              supertriangle vertices in subroutine delaun
*     Y      - Y-coords of all points in data set
*            - Y-coord of point I given by Y(I)
*            - Last three locations are used to store y-coords of
*              supertriangle vertices in subroutine delaun
*     V      - Vertex array for triangulation
*            - Vertices listed in anticlockwise sequence
*            - Vertices for triangle J are found in V(I,J) for I=1,2,3
*              and J=1,2,...,NTRI
*            - First vertex is at point of contact of first and third