contri.f

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

      SUBROUTINE BCLIP(NPTS,NCB,BLIST,TN,V,T,NTRI,NB,maxelist)
      use mod_output
      implicit double precision (a-h,o-z)
***********************************************************************
*
*     PURPOSE:
*     --------
*
*     Clip constrained Delaunay triangulation to boundaries in BLIST
*     If BLIST is empty, then the triangulation is clipped to a convex
*     hull by removing all triangles that are formed with the 
*     supertriangle vertices
*     The triangulation MUST initially include the supertriangle 
*     vertices
*
*     INPUT:
*     ------
*
*     NPTS   - Total number of points in data set (NPTS ge N)
*     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
*     BLIST  - List of edges which must be present in triangulation
*              and define boundaries
*            - Edge I defined by vertices in BLIST(1,I) and BLIST(2,I)
*              where I ranges from 1,...,NCB
*            - 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
*     TN     - List of triangle numbers such that vertex I can be 
*              found in triangle TN(I)
*     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   - Number of triangles, including those formed with
*              vertices of supertriangle
*     NB     - Not defined
*
*     OUTPUT:
*     -------
*
*     NPTS   - Unchanged
*     NCB    - Unchanged
*     BLIST  - Unchanged
*     TN     - Not defined 
*     V      - Updated 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      - Updated 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   - Updated number of triangles in final triangulation
*     NB     - Number of boundaries defining the triangulation
*            - NB=1 for a simply connected domain with no holes 
*            - NB=H+1 for a domain with H holes
*
*     PROGRAMMER:           Scott Sloan
*     -----------
*
*     LAST MODIFIED:        3 march 1991        Scott Sloan
*     --------------
*
*     COPYRIGHT 1990:       Scott Sloan
*     ---------------       Department of Civil Engineering
*                           University of Newcastle
*                           NSW 2308
*
***********************************************************************
      integer maxelist
      INTEGER A,C,E,I,J,L,R,S
      INTEGER NB,TS,V1,V2
      INTEGER NCB,NEV
      INTEGER NPTS,NTRI
      INTEGER NNTRI,NPTS1,NTRI1
      INTEGER T(3,NTRI),V(3,NTRI)
      INTEGER TN(NPTS+3)
      INTEGER BLIST(2,maxelist)
*---------------------------------------------------------------------
*     Skip to triangle deletion step if triangulation has no
*     boundary constraints
*
      IF(NCB.EQ.0)THEN
        NB=1
        GOTO 55
      ENDIF
*---------------------------------------------------------------------
*     Mark all edges which define the boundaries
*     S=triangle in which search starts
*     R=triangle to right of edge V1-V2
*     
      DO 20 I=1,NCB
        V1=BLIST(1,I)
        V2=BLIST(2,I)
        S=TN(V1)
        R=S
*
*       Circle anticlockwise round V1 until edge V1-V2 is found
*
   10   IF(V(1,R).EQ.V1)THEN
          IF(V(3,R).EQ.V2)THEN
            T(3,R)=-T(3,R)
            GOTO 20
          ELSE
            R=ABS(T(3,R))
          ENDIF
        ELSEIF(V(2,R).EQ.V1)THEN
          IF(V(1,R).EQ.V2)THEN
            T(1,R)=-T(1,R)
            GOTO 20
          ELSE
            R=ABS(T(1,R))
          ENDIF
        ELSEIF(V(2,R).EQ.V2)THEN
            T(2,R)=-T(2,R)
            GOTO 20
        ELSE
            R=ABS(T(2,R))
        ENDIF
        IF(R.EQ.S)THEN
          nerror = 4
          error(1) = 
     &      'Constrained boundary edge not found in subroutine BCLIP'
          write(error(2),'(''V1='',i5,'' V2='',i5)')V1,V2
          error(3) = 'check that this edge is not crossed by another'
          error(4) = 'boundary edge.'
          call erroroutput
          STOP
        ENDIF
        GOTO 10
   20 CONTINUE
*--------------------------------------------------------------------
*     Mark all triangles which are to right of boundaries by
*     circling anticlockwise around the outside of each boundary
*     Loop while some boundary edges have not been visited
*     S = triangle from which search starts
*     NEV = number of edges visited
*
      S  =0
      NB =0
      NEV=0
      NTRI1=NTRI+1
      NPTS1=NPTS+1
   30 IF(NEV.LT.NCB)THEN
*
*       Find an edge on a boundary
*
   40   S=S+1
        IF(T(1,S).LT.0)THEN
          E=1
        ELSEIF(T(2,S).LT.0)THEN
          E=2
        ELSEIF(T(3,S).LT.0)THEN
          E=3
        ELSE
          GOTO 40
        ENDIF
*
*       Store and mark starting edge
*       Mark starting triangle for deletion 
*       Increment count of edges visited
*       C = current triangle
*
        TS =T(E,S)
        T(E,S)=NTRI1
        V(1,S)=NPTS1
        NEV=NEV+1
        C=S
*
*       Increment to next edge
*
        E=MOD(E+1,3)+1
*
*       Loop until trace of current boundary is complete
*
   50   IF(T(E,C).NE.NTRI1)THEN
          IF(T(E,C).LT.0)THEN
*
*           Have found next boundary edge
*           Increment count of boundary edges visited, unmark the edge
*           and move to next edge 
*
            NEV=NEV+1
            T(E,C)=-T(E,C)
            E=MOD(E+1,3)+1
          ELSE
*
*           Non-boundary edge, circle anticlockwise around boundary 
*           vertex to move to next triangle, and mark next
*           triangle for deletion
*
            L=C
            C=T(E,L)
            IF(T(1,C).EQ.L)THEN
              E=3
            ELSEIF(T(2,C).EQ.L)THEN
              E=1
            ELSE
              E=2
            ENDIF
            V(1,C)=NPTS1
          ENDIF
          GOTO 50
        ENDIF
*
*       Trace of current boundary is complete
*       Reset marked edge to original value and check for any more
*       boundaries
*
        T(E,C)=-TS
        NB=NB+1
        GOTO 30
      ENDIF
*---------------------------------------------------------------------
*     Remove all triangles which have been marked for deletion
*     Any triangle with a vertex number greater than NPTS is deleted
*     
   55 CONTINUE
      NNTRI=NTRI
      NTRI =0
      DO 80 J=1,NNTRI
        IF(MAX(V(1,J),V(2,J),V(3,J)).GT.NPTS)THEN
*
*         Triangle J is to be deleted
*         Update adjacency lists for triangles adjacent to J
* 
          DO 60 I=1,3
            A=T(I,J)
            IF(A.NE.0)THEN
              IF(T(1,A).EQ.J)THEN
                T(1,A)=0
              ELSEIF(T(2,A).EQ.J)THEN
                T(2,A)=0
              ELSE
                T(3,A)=0
              ENDIF
            ENDIF
   60     CONTINUE
        ELSE
*
*         Triangle J is not to be deleted
*         Update count of triangles
*
          NTRI=NTRI+1
          IF(NTRI.LT.J)THEN
*
*           At least one triangle has already been deleted
*           Relabel triangle J as triangle NTRI
*
            DO 70 I=1,3
              A=T(I,J)
              T(I,NTRI)=A
              V(I,NTRI)=V(I,J)
              IF(A.NE.0)THEN
                IF(T(1,A).EQ.J)THEN
                  T(1,A)=NTRI
                ELSEIF(T(2,A).EQ.J)THEN
                  T(2,A)=NTRI
                ELSE
                  T(3,A)=NTRI
                ENDIF
              ENDIF
  70        CONTINUE
          ENDIF
        ENDIF
  80  CONTINUE 
      END
************************************************************************
      SUBROUTINE BSORT(N,NPTS,X,Y,XMIN,XMAX,YMIN,YMAX,DMAX,BIN,LIST,W)
      implicit double precision (a-h,o-z)
************************************************************************
*
*     PURPOSE:
*     --------
*
*     Sort points such that consecutive points are close to one another 
*     in the x-y plane using a bin sort
*
*     INPUT:
*     ------
*
*     N      - Total number of points to be triangulated (N le NPTS)
*     NPTS   - Total number of points in data set
*     X      - X-coords of all points in data set
*            - If point is in list,the coordinate must be normalised
*              according to X=(X-XMIN)/DMAX
*            - 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
*            - If point is in list,the coordinate must be normalised
*              according to Y=(Y-YMIN)/DMAX
*            - Y-coord of point I given by Y(I)
*            - Last three locations are used to store y-coords of
*              supertriangle vertices in subroutine delaun
*     XMIN   - Min x-coord of points in LIST
*     XMAX   - Max x-coord of points in LIST
*     YMIN   - Min y-coord of points in LIST
*     YMAX   - Max y-coord of points in LIST
*     DMAX   - DMAX=MAX(XMAX-XMIN,YMAX-YMIN)
*     BIN    - Not defined
*     LIST   - List of points to be triangulated
*     W      - Undefined workspace
*
*     OUTPUT:
*     -------
*
*     N      - Unchanged
*     NPTS   - Unchanged
*     X      - Unchanged
*     Y      - Unchanged
*     XMIN   - Unchanged
*     XMAX   - Unchanged
*     YMIN   - Unchanged
*     YMAX   - Unchanged
*     DMAX   - Unchanged
*     BIN    - Not used
*     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 used
*
*     PROGRAMMER:             Scott Sloan
*     -----------
*
*     LAST MODIFIED:          3 march 1991        Scott Sloan
*     --------------
*
*     COPYRIGHT 1990:         Scott Sloan
*     ---------------         Department of Civil Engineering
*                             University of Newcastle
*                             NSW 2308
*
************************************************************************
      INTEGER B,I,J,K,N,P
      INTEGER LI,LT,NB
      INTEGER NDIV,NPTS
      INTEGER W(N)
      INTEGER BIN(NPTS)
      INTEGER LIST(N)
*
      DOUBLE PRECISION DMAX,XMAX,XMIN,YMAX,YMIN
      DOUBLE PRECISION FACTX,FACTY
      DOUBLE PRECISION X(NPTS+3),Y(NPTS+3)
*---------------------------------------------------------------------
*     Compute number of bins in x-y directions
*     Compute inverse of bin size in x-y directions
*
      NDIV=NINT(REAL(N)**0.25d0)
      FACTX=REAL(NDIV)/((XMAX-XMIN)*1.01d0/DMAX)
      FACTY=REAL(NDIV)/((YMAX-YMIN)*1.01d0/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)