contri.f

Programs in the irregular grid design package described in this manual are used to carry out five ma

*              to one another in the x-y plane
*     W      - Not used
*
*     PROGRAMMER:             Scott Sloan
*     -----------
*
*     LAST MODIFIED:          3 march 1991        Scott Sloan
*     --------------
*
*     COPYRIGHT 1990:         Scott Sloan
*     ---------------         Department of Civil Engineering
*                             University of Newcastle
*                             NSW 2308
c
c Note: AG Dolling, April 1998. This routine is an implementation of
c the bucket sort of Cormen et al, Introduction to Algorithms,
c 1990. MIT Press. Sec 9.4, p 180.
c
*
************************************************************************
      INTEGER B,I,J,K,N,P
      INTEGER LI,LT,NB
      INTEGER NDIV,NPTS
      INTEGER W(N)
      INTEGER BIN(NPTS)
      INTEGER LIST(N)
*
      REAL DMAX,XMAX,XMIN,YMAX,YMIN
      REAL FACTX,FACTY
      REAL X(NPTS+3),Y(NPTS+3)
	BSORT = .FALSE.
*---------------------------------------------------------------------
*     Compute number of bins in x-y directions
*     Compute inverse of bin size in x-y directions
*
      NDIV=NINT(REAL(N)**0.25)
      FACTX=REAL(NDIV)/((XMAX-XMIN)*1.01/DMAX)
      FACTY=REAL(NDIV)/((YMAX-YMIN)*1.01/DMAX)
*---------------------------------------------------------------------
*     Zero count of points in each bin
*
      NB=NDIV*NDIV
      DO 5 I=1,NB
        W(I)=0
    5 CONTINUE
*---------------------------------------------------------------------
*     Assign bin numbers to each point
*     Count entries in each bin and store in W
*
      DO 10 K=1,N
        P=LIST(K)
        I=INT(Y(P)*FACTY)
        J=INT(X(P)*FACTX)
        IF(MOD(I,2).EQ.0)THEN
          B=I*NDIV+J+1
        ELSE
          B=(I+1)*NDIV-J
        ENDIF
        BIN(P)=B
        W(B)=W(B)+1
   10 CONTINUE
*---------------------------------------------------------------------
*     Form pointers to end of each bin in final sorted list
*     Note that some bins may contain no entries
*
      DO 20 I=2,NB
        W(I)=W(I-1)+W(I)
   20 CONTINUE
*---------------------------------------------------------------------
*     Now perform linear sort
*
      DO 40 I=1,N
        IF(LIST(I).GT.0)THEN
          LI=LIST(I)
          B=BIN(LI)
          P=W(B)
   30     IF(P.NE.I)THEN
            W(B)=W(B)-1
            LT =LIST(P)
            LIST(P)=LI
            LI=LT
            B=BIN(LI)
            LIST(P)=-LIST(P)
            P=W(B)
            GOTO 30
          ENDIF
          W(B)=W(B)-1
          LIST(I)=LI
        ELSE
          LIST(I)=-LIST(I)
        ENDIF
   40 CONTINUE
c next line added by ag dolling 30 oct 1997
	bin(1) = nb
	BSORT = .TRUE.
      END
c----------
      FUNCTION DELAUN(NPTS,N,X,Y,LIST,STACK,V,T,NTRI)
	LOGICAL DELAUN
************************************************************************
*
*     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)
*
      REAL D
      REAL XP,YP
      REAL TOL,X13,X23,X1P,X2P,Y13,Y23,Y1P,Y2P
      REAL COSA,COSB
      REAL C00000,C00100
      REAL 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)
*
C      PARAMETER (TOL=1.E-10)   !32 bit reals
      PARAMETER (TOL=1.E-16)    !64 bit reals
      PARAMETER (C00000=0.0, C00100=100.0)

	DELAUN = .FALSE.
*---------------------------------------------------------------------
*     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
        IF(D.LT.TOL)THEN
           D=D+(Y(V1)-YP)**2
           IF(D.LT.TOL)THEN
             WRITE(66,2000)V1,P
           ENDIF
        ENDIF
        D=(X(V2)-XP)**2
        IF(D.LT.TOL)THEN
           D=D+(Y(V2)-YP)**2
           IF(D.LT.TOL)THEN
             WRITE(66,2000)V2,P
           ENDIF
        ENDIF
        D=(X(V3)-XP)**2
        IF(D.LT.TOL)THEN
           D=D+(Y(V3)-YP)**2
           IF(D.LT.TOL)THEN
             WRITE(66,2000)V3,P
           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