contri.f

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

      FUNCTION CONTRI(NPTS,N,NCE,NCB,ELIST,X,Y,LIST,W,V,T,NTRI,
     +    nerrors, maxerrors, err_elist)
	LOGICAL CONTRI
***********************************************************************
*
*     PURPOSE:
*     --------
*
*     Assemble constrained Delaunay triangulation for collection of 
*     points in the plane. This code has a total memory requirement
*     equal to 4*NPTS + 13*N + 2*NCE + 18
*
*     INPUT:
*     ------
*
*     NPTS   - Total number of points in data set (NPTS ge N)
*     N      - Total number of points to be triangulated (N ge 3)
*     NCE    - Total number of constrained edges which must be 
*              present, including those which define boundaries
*            - NCE=0 indicates triangulation is unconstrained so that
*              ELIST is empty and the code will produce a 
*              triangulation of the convex hull
*     NCB    - Total number of constrained edges which define one
*              external boundary and any internal boundaries (holes)
*            - NCB=0 indicates there are no boundary edge constraints
*              and the code will produce a triangulation of the convex
*              hull
*            - NCB must not be greater than NCE
*            - If NCB gt 0, then at least one of the boundaries
*              specified must be external
*     ELIST  - List of edges which must be present in triangulation
*            - These may be internal edges or edges which define 
*              boundaries
*            - Edge I defined by vertices in ELIST(1,I) and ELIST(2,I)
*            - Edges which define boundaries must come first in ELIST
*              and thus occupy the first NCB columns
*            - Edges which define an external boundary must be listed
*              anticlockwise but may be presented in any order
*            - Edges which define an internal boundary (hole) must be 
*              listed clockwise but may be presented in any order
*            - An internal boundary (hole) cannot be specified unless
*              an external boundary is also specified
*            - All boundaries must form closed loops
*            - An edge may not appear more than once in ELIST
*            - An external or internal boundary may not cross itself
*              and may not share a common edge with any other boundary
*            - Internal edges, which are not meant to define boundaries 
*              but must be present in the final triangulation, occupy
*              columns NCB+1,... ,NCE of ELIST
*     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
*     LIST   - List of points to be triangulated
*            - If N eq NPTS, set LIST(I)=I for I=1,2,...,NPTS
*              prior to calling this routine
*            - No point in LIST may lie outside any external boundary
*              or inside any internal boundary 
*     W      - Not defined, workspace of length ge 2*(NPTS+3)
*     V      - Not defined
*     T      - Not defined
*     NTRI   - Not defined
*
*     OUTPUT:
*     -------
*
*     NPTS   - Unchanged
*     N      - Unchanged
*     NCE    - Unchanged
*     NCB    - Unchanged
*     ELIST  - Unchanged
*     X      - Unchanged, except that last three locations contain
*              normalised x-coords of supertriangle vertices
*     Y      - Unchanged, except that last three locations contain
*              normalised y-coords of supertriangle vertices
*     LIST   - List of points to be triangulated
*            - Points ordered such that consecutive points are close
*              to one another in the x-y plane
*     W      - 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
*              and J=1,2,...,NTRI
*            - Adjacent triangles listed in anticlockwise sequence
*            - Zero denotes no adjacent triangle
*     NTRI   - Total number of triangles in final triangulation
*
*     SUBROUTINES  CALLED:  BSORT, DELAUN, EDGE, TCHECK, BCLIP
*     --------------------
*
*     PROGRAMMER:             Scott Sloan
*     -----------
*
*     LAST MODIFIED:          3 march 1991        Scott Sloan
*     --------------
*
*     COPYRIGHT 1990:         Scott Sloan
*     ---------------         Department of Civil Engineering
*                             University of Newcastle
*                             NSW 2308
*
***********************************************************************
      INTEGER I,J,N,P
      INTEGER JW,NB,VI,VJ
      INTEGER NCB,NCE,NEF
      INTEGER NPTS,NTRI
      INTEGER T(3,2*N+1),V(3,2*N+1),W(NPTS+3,2)
      INTEGER LIST(N)
      INTEGER ELIST(2,NCE)
c      INTEGER ELIST(2,*)

	integer nerrors, maxerrors
	integer err_elist(maxerrors)
*
      REAL DMAX,XMIN,XMAX,YMIN,YMAX
      REAL C00001
      REAL FACT
      REAL X(NPTS+3),Y(NPTS+3)
*
      PARAMETER(C00001=1.0)

	integer start, end, nbins, JJ, L1, L2
	LOGICAL BCLIP
	LOGICAL BSORT
	LOGICAL DELAUN
	LOGICAL EDGE
	LOGICAL TCHECK

	CONTRI = .FALSE.
*---------------------------------------------------------------------
*     Check input for obvious errors
*
      IF(NPTS.LT.3)THEN
        nerrors = -1
        call PigPutMessage('Less than 3 points.')
        WRITE(66,'(//,''***ERROR IN CONTRI***'',
     +             /,''LESS THAN 3 POINTS IN DATASET'',
     +             /,''CHECK VALUE OF NPTS'')')
        RETURN !STOP
      ENDIF
      IF(N.LT.3)THEN
        nerrors = -1
        call PigPutMessage('Less than 3 points.')
        WRITE(66,'(//,''***ERROR IN CONTRI***'',
     +             /,''LESS THAN 3 POINTS TO BE TRIANGULATED'',
     +             /,''CHECK VALUE OF N'')')
        RETURN !STOP
      ELSEIF(N.GT.NPTS)THEN
        nerrors = -1
        call PigPutMessage('Inconsistent database. Check gridit.log.')
        WRITE(66,'(//,''***ERROR IN CONTRI***'',
     +             /,''TOO MANY POINTS TO BE TRIANGULATED'',
     +             /,''N MUST BE LESS THAN OR EQUAL TO NPTS'',
     +             /,''CHECK VALUES OF N AND NPTS'')')
        RETURN !STOP
      ENDIF
      IF(NCB.GT.NCE)THEN
        nerrors = -1
        call PigPutMessage('Inconsistent database. Check gridit.log.')
        WRITE(66,'(//,''***ERROR IN CONTRI***'',
     +             /,''NCB GREATER THAN NCE'',
     +             /,''CHECK BOTH VALUES'')')
        RETURN !STOP
      ENDIF
*---------------------------------------------------------------------
*     Check for invalid entries in LIST
*
      DO 10 I=1,N
        P=LIST(I)
        IF(P.LT.1.OR.P.GT.NPTS)THEN
          nerrors = -1
          call PigPutMessage('Inconsistent database. Check gridit.log.')
          WRITE(66,1000)I,P
          RETURN !STOP
        ENDIF
        W(P,1)=0
   10 CONTINUE  
      DO 20 I=1,N
        P=LIST(I)
        W(P,1)=W(P,1)+1
   20 CONTINUE
      DO 30 I=1,N
        P=LIST(I)
        IF(W(P,1).GT.1)THEN
          nerrors = -1
          call PigPutMessage('Inconsistent database. Check gridit.log.')
          WRITE(66,'(//,''***ERROR IN CONTRI***'',
     +               /,''POINT'',i8,'' OCCURS'',i8,'' TIMES IN LIST'',
     +               /,''POINT SHOULD APPEAR ONLY ONCE'')')P,W(P,1)
          RETURN !STOP
        ENDIF
   30 CONTINUE
*---------------------------------------------------------------------
*     Check for invalid entries in ELIST
*
      IF(NCE.GT.0)THEN
        DO 40 I=1,NPTS
          W(I,1)=0
   40   CONTINUE
        DO 50 I=1,N
          W(LIST(I),1)=1
   50   CONTINUE
        DO 60 I=1,NCE
          VI=ELIST(1,I)
          IF(VI.LT.1.OR.VI.GT.NPTS)THEN
            nerrors = -1
            call PigPutMessage(
     +        'Inconsistent database. Check gridit.log.')
            WRITE(66,2000)1,I,VI
            RETURN !STOP
          ELSEIF(W(VI,1).NE.1)THEN
            nerrors = -1
            call PigPutMessage(
     +        'Inconsistent database. Check gridit.log.')
            WRITE(66,3000)1,I,VI
            RETURN !STOP
          ENDIF
          VJ=ELIST(2,I)
          IF(VJ.LT.1.OR.VJ.GT.NPTS)THEN
            nerrors = -1
            call PigPutMessage(
     +        'Inconsistent database. Check gridit.log.')
            WRITE(66,2000)2,I,VJ
            RETURN !STOP
          ELSEIF(W(VJ,1).NE.1)THEN
            nerrors = -1
            call PigPutMessage(
     +        'Inconsistent database. Check gridit.log.')
            WRITE(66,3000)2,I,VJ
	      RETURN !STOP
          ENDIF
   60   CONTINUE
      ENDIF
*---------------------------------------------------------------------
*     Check that all boundaries (if there are any) form closed loops
*     Count appearances in ELIST(1,.) and ELIST(2,.) of each node
*     These must match if each boundary forms a closed loop
*
      IF(NCB.GT.0)THEN
        DO 70 I=1,NPTS
          W(I,1)=0
          W(I,2)=0
   70   CONTINUE
        DO 80 I=1,NCB
          VI=ELIST(1,I)
          VJ=ELIST(2,I)
          W(VI,1)=W(VI,1)+1
          W(VJ,2)=W(VJ,2)+1
   80   CONTINUE
        DO 90 I=1,NCB
          VI=ELIST(1,I)
          IF(W(VI,1).NE.W(VI,2))THEN
            nerrors = -1
            call PigPutMessage(
     +        'Inconsistent database. Check gridit.log.')
            WRITE(66,'(//,''***ERROR IN CONTRI***'',
     +                 /,''BOUNDARY SEGMENTS DO NOT FORM A '',
     +                   ''CLOSED LOOP'',
     +                 /,''CHECK ENTRIES IN COLS 1...NCB OF ELIST '',
     +                   ''FOR NODE'',i8)')VI
			RETURN !STOP
          ENDIF
   90   CONTINUE        
      ENDIF
*---------------------------------------------------------------------
*     Compute min and max coords for x and y
*     Compute max overall dimension
*
      P=LIST(1)
      XMIN=X(P)
      XMAX=XMIN
      YMIN=Y(P)
      YMAX=YMIN
      DO 100 I=2,N
        P=LIST(I)
        XMIN=MIN(XMIN,X(P))
        XMAX=MAX(XMAX,X(P))
        YMIN=MIN(YMIN,Y(P))
        YMAX=MAX(YMAX,Y(P))
  100 CONTINUE
      IF(XMIN.EQ.XMAX)THEN
        nerrors = -1
        call PigPutMessage('All points have same X-Coordinate.')
        WRITE(66,'(//,''***ERROR IN CONTRI***'',
     +             /,''ALL POINTS HAVE SAME X-COORD'',
     +             /,''ALL POINTS ARE COLLINEAR'')')
        RETURN !STOP
      ENDIF
      IF(YMIN.EQ.YMAX)THEN
        nerrors = -1
        call PigPutMessage('All points have same Y-Coordinate.')
        WRITE(66,'(//,''***ERROR IN CONTRI***'',
     +             /,''ALL POINTS HAVE SAME Y-COORD'',
     +             /,''ALL POINTS ARE COLLINEAR'')')
        RETURN !STOP
      ENDIF
      DMAX=MAX(XMAX-XMIN,YMAX-YMIN)
*---------------------------------------------------------------------
*     Normalise x-y coords of points
*
      FACT=C00001/DMAX
      DO 110 I=1,N
        P=LIST(I)
        X(P)=(X(P)-XMIN)*FACT
        Y(P)=(Y(P)-YMIN)*FACT
  110 CONTINUE
*---------------------------------------------------------------------
*     Sort points into bins using a linear bin sort
*     This call is optional
*
      IF(.NOT. 
     +	BSORT(N,NPTS,X,Y,XMIN,XMAX,YMIN,YMAX,DMAX,W,LIST,W(1,2))
     +	) THEN
		RETURN
	END IF
c
c Add check for coincident points - adjacent points in list must not be
c at same location
c This test added by ag Dolling 30 Oct 1997
c depends on the results of BSORT, including its use of work areas
c W(1,1) as the number of bins used for the sort, and
c W(*,2) as the count of points in each of the W(1,1) bins.
c
      CONTRI = .true.
	nbins = W(1,1)
	start = 1
	end = 0
	do jj=1,nbins
		start = end + 1
		end = w(jj,2)
		do i = start,end
			do j = i+1, end
				l1 = list(i)
				l2 = list(j)
c	    call PigDrawSymbols(1, x(l1)*DMAX+XMIN, y(l1)*DMAX+YMIN)
c          write(66,'(2(i8,g24.15))') l1,x(l1),y(l1),l2,x(l2),y(l2)
			  if  (        (x(l1).eq.x(l2))
     +			    .and.    (y(l1).eq.y(l2))
     +		        ) then
	        if(nerrors.eq.0)then
                call PigPutMessage(
     +'Coincident points detected. Markers being placed...')
              end if
              nerrors = nerrors - 1 !note: negative count treated specially...
					write(66,'(a,i8,a,i8)') 'Points ',l1,' and ',l2,
     +		            ' are coincident.'
					call addmarker(x(l1)*DMAX+XMIN,y(l1)*DMAX+YMIN)
					CONTRI = .false.
				end if
			end do
		end do
	end do
      if(.not.CONTRI)then
*---------------------------------------------------------------------
*     Reset x-y coords to original values
*
          call PigPutMessage(
     +'Coincident points detected. Delete duplicates at markers.')
          DO I=1,N
            P=LIST(I)
            X(P)=X(P)*DMAX+XMIN
            Y(P)=Y(P)*DMAX+YMIN
	    end do
c loop below added by a g dolling 17 October 1997
c to return supertriangle in original coordinates
          DO P=NPTS+1,NPTS+3
            X(P)=X(P)*DMAX+XMIN
            Y(P)=Y(P)*DMAX+YMIN
          end do
c	    call PigUWait(2.0)
	    return
	endif
	CONTRI = .FALSE.
*---------------------------------------------------------------------
*     Compute Delaunay triangulation
*
      IF(.NOT. DELAUN(NPTS,N,X,Y,LIST,W,V,T,NTRI)) THEN
        call PigPutMessage('Triangulation error. See gridit.log.')
        call PigUWait(3.0)
		RETURN
	END IF
*---------------------------------------------------------------------
*     For each node, store any triangle in which it is a vertex
*     Include supertriangle vertices
*
      DO 130 J=1,NTRI
        DO 120 I=1,3
          VI=V(I,J)
          W(VI,1)=J
  120   CONTINUE
  130 CONTINUE
*---------------------------------------------------------------------
*     Constrain triangulation by forcing edges to be present
*
      NEF=0
      JW=(NPTS+3)/2
      DO 140 I=1,NCE
        VI=ELIST(1,I)