contri.f

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

* 
*             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
                  WRITE(66,1000)
                  RETURN !STOP
                ENDIF
                STACK(TOPSTK)=R
              ENDIF
              IF(B.NE.0)THEN
                TOPSTK=TOPSTK+1
                IF(TOPSTK.GT.NPTS)THEN
                  WRITE(66,1000)
                  RETURN !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
        WRITE(66,'(//,''***ERROR IN DELAUN***'',
     +             /,''INCORRECT NUMBER OF TRIANGLES FORMED'')')
        RETURN !STOP
      ENDIF
	DELAUN = .TRUE.
*---------------------------------------------------------------------
 1000 FORMAT(//,'***ERROR IN SUBROUTINE DELAUN***',
     +        /,'STACK OVERFLOW')
 2000 FORMAT(//,'***WARNING IN DELAUN***',
     +        /,'POINTS',i8,' AND',i8,' ARE COINCIDENT',
     +        /,'DELETE EITHER POINT FROM LIST VECTOR' )
      END
c----------
      FUNCTION EDGE(VI,VJ,NPTS,NTRI,NEF,JW,X,Y,V,T,TN,W,
     +      XMIN,YMIN,DMAX)
	LOGICAL EDGE
************************************************************************
*
*     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
*              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
*     XMIN   - Min x-coord of points in LIST
*     YMIN   - Min y-coord of points in LIST
*     DMAX   - Max of range in x-coord or y-coord of points in LIST
*
*     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
*     XMIN   - Unchanged
*     YMIN   - Unchanged
*     DMAX   - Unchanged
*
*     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
*
************************************************************************
      REAL XMIN, YMIN, DMAX
      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)
*
      REAL X1,X2,X3,X4,XI,XJ,Y1,Y2,Y3,Y4,YI,YJ
      REAL TOL,X13,X14,X23,X24,Y13,Y14,Y24,Y23
      REAL COSA,COSB
      REAL C00000,DETIJ3
      REAL X(NPTS+3),Y(NPTS+3)
* 
      LOGICAL SWAP
*
c      PARAMETER (TOL=1.E-6)  !32 bit
      PARAMETER (TOL=1.E-16)   !64 bit
      PARAMETER (C00000=0.0)

	EDGE = .FALSE.
*---------------------------------------------------------------------
*     Check data
*
      IF(VI.LE.0.OR.VI.GT.NPTS)THEN
        call PigPutMessage('Triangulation error. See gridit.log.')
        call PigUWait(3.0)
        WRITE(66,'(//,''***ERROR IN EDGE***'',
     +             /,''ILLEGAL VALUE FOR VI'',
     +             /,''VI='',i8)')VI
        RETURN !STOP
      ENDIF
      IF(VJ.LE.0.OR.VJ.GT.NPTS)THEN
        call PigPutMessage('Triangulation error. See gridit.log.')
        call PigUWait(3.0)
        WRITE(66,'(//,''***ERROR IN EDGE***'',
     +             /,''ILLEGAL VALUE FOR VJ'',
     +             /,''VJ='',i8)')VJ
        RETURN !STOP
      ENDIF
*---------------------------------------------------------------------- 
*     Find any triangle which has VI as a vertex
*
      S=TN(VI)
      IF(S.LE.0)THEN
        call PigPutMessage('Triangulation error. See gridit.log.')
        call PigUWait(3.0)
        WRITE(66,1000)VI
        RETURN !STOP
      ENDIF           
      IF(TN(VJ).LE.0)THEN
        call PigPutMessage('Triangulation error. See gridit.log.')
        call PigUWait(3.0)
        WRITE(66,1000)VJ
        RETURN !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)THEN
		EDGE = .TRUE.
		RETURN
	END IF
*
*     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
      call addmarker(x(V1)*DMAX+XMIN,y(V1)*DMAX+YMIN) !experimental
c      call addmarker(x(VJ)*DMAX+XMIN,y(VJ)*DMAX+YMIN) !experimental
c next message wrong - changed below agd 1998/may/25
c      WRITE(66,'(//,''***ERROR IN EDGE***'',
c     +  /,''VERTEX ADJACENT TO '',i8,'' IS ON ARC'',
c     +  /,''BETWEEN VERTICES'',i8,'' AND'',i8)')VI,VI,VJ
      WRITE(66,'(//,''***ERROR IN EDGE***'',
     +  /,''VERTEX '',i8,'' ADJACENT TO '',i8,'' IS ON ARC'',
     +  /,''BETWEEN VERTICES'',i8,'' AND'',i8)')V1,V2,VI,VJ
      RETURN !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
        WRITE(66,'(//,''***ERROR IN EDGE***'',
     +             /,''NOT ENOUGH WORKSPACE'',
     +             /,''INCREASE JW'')')
        RETURN !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 
        call addmarker(x(V3)*DMAX+XMIN,y(V3)*DMAX+YMIN) !experimental
        WRITE(66,'(//,''***ERROR IN EDGE***'',
     +             /,''VERTEX'',i8,'' IS ON ARC'',
     +             /,''BETWEEN VERTICES'',i8,'' AND'',i8)')V3,VI,VJ
        RETURN !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
              WRITE(66,'(//,''***ERROR IN EDGE***'',
     +                   /,''ARC BETWEEN VERTICES'',i8,'' AND'',i8,
     +                   /,''CROSSES SUPERTRIANGLE BOUNDARY DEFINED '',
     +                     ''BY VERTICES'',i8,''AND'',i8)')VI,VJ,V1,V2
               RETURN !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)