contri.f

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

!-------------------------------------------------------------------------------
!
!  Quasicontinuum (QC) Method: Mixed Continuum and Atomistic Simulation Package
!  QC Package distribution version 1.2 (February 2005)    
!
!  Copyright (C) 2003 R. Miller, M. Ortiz, R. Phillips, D. Rodney, E. B. Tadmor
!  Copyright (C) 2004, 2005 R. Miller, E. B. Tadmor
!
!  This program is free software; you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation; either version 2 of the License, or
!  (at your option) any later version.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with this program; if not, write to the Free Software
!  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
!
!  For more information please visit the QC website at:
!
!    www.qcmethod.com
!
!  or contact Ron Miller or Ellad Tadmor with contact information below:
!
!  NOTE HOWEVER THAT NEITHER RON MILLER NOR ELLAD TADMOR WILL BE ABLE TO
!  PROVIDE SUPPORT OR DEBUGGING SERVICES FOR THE CODE. THE CODE IS SUPPLIED 
!  AS IS FOR YOUR USE.
!
!    Ron Miller
!    Department of Mechanical and Aerospace Engineering
!    Carleton University
!    1125 Colonel By Drive
!    Ottawa, ON, K1S 5B6
!    CANADA
!   
!    email: [e-mail]
!    URL: http://www.mae.carleton.ca/~rmiller
!
!  -or-
!
!    Ellad B. Tadmor
!    Faculty of Mechanical Engineering
!    Technion - Israel Institute of Technology
!    32000 Haifa
!    ISRAEL
!    
!    email: [e-mail]
!    URL: http://techunix.technion.ac.il/~tadmor/
!   
!-------------------------------------------------------------------------------

      SUBROUTINE CONTRI(NPTS,N,NCE,NCB,ELIST,XX,NXDM,LIST,W,V,T,NTRI,
     &    MAXELIST,MAXXX)
      use mod_output
      implicit double precision (a-h,o-z)
***********************************************************************
*
*     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 maxelist,maxxx
      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,maxelist)
*
      DOUBLE PRECISION DMAX,XMIN,XMAX,YMIN,YMAX
      DOUBLE PRECISION C00001
      DOUBLE PRECISION FACT
      integer nxdm
      double precision xx(nxdm,maxxx)
      double precision, pointer:: x(:),y(:)
*
      PARAMETER(C00001=1.d0)
c
c Ron's addition: just putting xx(nxdm,1) into X and Y
c
      allocate(x(npts+3),y(npts+3))
      do i=1,npts
         x(i)=xx(1,i)
         y(i)=xx(2,i)
      enddo
*---------------------------------------------------------------------
*     Check input for obvious errors
*
      IF(NPTS.LT.3)THEN
        nerror = 2
        error(1) = 'Less than 3 points in dataset in subroutine CONTRI'
        error(2) = 'check value of NPTS'
        call erroroutput
        STOP
      ENDIF
      IF(N.LT.3)THEN
        nerror = 2
        error(1) = 
     &   'Less than 3 points to be triangulated in subroutine CONTRI'
        error(2) = 'check value of N'
        call erroroutput
        STOP
      ELSEIF(N.GT.NPTS)THEN
        nerror = 2
        error(1) = 
     &   'Too many points to be triangulated in subroutine CONTRI'
        error(2) = 'N must be less than or equal to NPTS'
        call erroroutput
        STOP
      ENDIF
      IF(NCB.GT.NCE)THEN
        nerror = 2
        error(1) = 'NCB greater than NCE in subroutine CONTRI'
        error(2) = 'check both values'
        call erroroutput
        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
          nerror = 2
          write(error(1),1000)
          write(error(2),1010)I,P
          call erroroutput
          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
          nerror = 2
          write(error(1),
     &       '(''Point'',i5,'' occurs'',i5,'' times in list'')')
     &       P,W(P,1)
          error(2) = 'in subroutine CONTRI.'
     &               // ' Point should appear only once.'
          call erroroutput
          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
            nerror = 2
            write(error(1),2000)
            write(error(2),2010)1,I,VI
            call erroroutput
            STOP
          ELSEIF(W(VI,1).NE.1)THEN
            nerror = 2
            write(error(1),3000)
            write(error(2),3010)1,I,VI
            call erroroutput
            STOP
          ENDIF
          VJ=ELIST(2,I)
          IF(VJ.LT.1.OR.VJ.GT.NPTS)THEN
            nerror = 2
            write(error(1),2000)
            write(error(2),2010)2,I,VJ
            call erroroutput
            STOP
          ELSEIF(W(VJ,1).NE.1)THEN
            nerror = 2
            write(error(1),3000)
            write(error(2),3010)2,I,VJ
            call erroroutput
            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
            nerror = 3
            error(1) = 'Boundary segments do not form a closed loop'
            error(2) = 'in subroutine CONTRI. Check entries in cols'
            write(error(3),'(''1...NCB of ELIST for node'',i5)') VI
            call erroroutput
            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
        nerror = 2
        error(1) = 'All points have same x-coord in subroutine CONTRI.'
        error(2) = 'All points are collinear.'
        call erroroutput
        STOP
      ENDIF
      IF(YMIN.EQ.YMAX)THEN
        nerror = 2
        error(1) = 'All points have same y-coord in subroutine CONTRI.'
        error(2) = 'All points are collinear.'
        call erroroutput
        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
*
      CALL BSORT(N,NPTS,X,Y,XMIN,XMAX,YMIN,YMAX,DMAX,W,LIST,W(1,2))
*---------------------------------------------------------------------
*     Compute Delaunay triangulation
*
      CALL DELAUN(NPTS,N,X,Y,LIST,W,V,T,NTRI)
*---------------------------------------------------------------------
*     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)
        VJ=ELIST(2,I)
        CALL EDGE(VI,VJ,NPTS,NTRI,NEF,JW,X,Y,V,T,W,W(1,2))
  140 CONTINUE
*---------------------------------------------------------------------
*     Clip triangulation and check it
*
      CALL BCLIP(NPTS,NCB,ELIST,W,V,T,NTRI,NB,maxelist)
      CALL TCHECK(NPTS,N,X,Y,LIST,V,T,NTRI,NEF,NB,NCE,NCB,ELIST,W,
     &            maxelist)
*---------------------------------------------------------------------
*     Reset x-y coords to original values
*
      DO 150 I=1,N
        P=LIST(I)
        X(P)=X(P)*DMAX+XMIN
        Y(P)=Y(P)*DMAX+YMIN
  150 CONTINUE
*---------------------------------------------------------------------
 1000 FORMAT('Illegal value in LIST in subroutine CONTRI:')
 1010 FORMAT('LIST(',I5,' )=',I5)
 2000 FORMAT('Illegal value in ELIST in subroutine CONTRI:')
 2010 FORMAT('ELIST(',I5,',',I5,' )=',I5)
 3000 FORMAT('Illegal value in ELIST in subroutine CONTRI:')
 3010 FORMAT('ELIST(',I5,',',I5,' )=',I5,' -- point is not in LIST')
      deallocate(x,y)
      RETURN
      END 
************************************************************************