function angle ( xa, ya, xb, yb, xc, yc ) !*****************************************************************************80 ! !! ANGLE computes the size of an angle in 2D. ! ! Discussion: ! ! This routine computes the interior angle in radians at vertex ! (XB,YB) of the chain formed by the directed edges from ! (XA,YA) to (XB,YB) to (XC,YC). The interior is to the ! left of the two directed edges. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) XA, YA, XB, YB, XC, YC, the vertex coordinates. ! ! Output, real ( kind = 8 ) ANGLE, the angle, between 0 and 2*PI. ! ANGLE is set to PI/2 in the degenerate case. ! implicit none real ( kind = 8 ) angle real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) t real ( kind = 8 ) tol real ( kind = 8 ) x1 real ( kind = 8 ) x2 real ( kind = 8 ) xa real ( kind = 8 ) xb real ( kind = 8 ) xc real ( kind = 8 ) y1 real ( kind = 8 ) y2 real ( kind = 8 ) ya real ( kind = 8 ) yb real ( kind = 8 ) yc tol = 100.0D+00 * epsilon ( tol ) x1 = xa - xb y1 = ya - yb x2 = xc - xb y2 = yc - yb t = sqrt ( ( x1 * x1 + y1 * y1 ) * ( x2 * x2 + y2 * y2 ) ) if ( t == 0.0D+00 ) then t = 1.0D+00 end if t = ( x1 * x2 + y1 * y2 ) / t if ( 1.0D+00 - tol < abs ( t ) ) then t = sign ( 1.0D+00, t ) end if angle = acos ( t ) if ( x2 * y1 - y2 * x1 < 0.0D+00 ) then angle = 2.0D+00 * pi - angle end if return end function angle3 ( u, v, rtolsq ) !*****************************************************************************80 ! !! ANGLE3 computes the size of a plane angle in 3D. ! ! Discussion: ! ! This routine computes the angle, in the range [0,PI], between two ! 3D vectors U and V. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) U(1:3), V(1:3), the vectors. ! ! Input, real ( kind = 8 ) RTOLSQ, the relative tolerance used to ! detect a zero vector, based on the square of the Euclidean length. ! ! Output, real ( kind = 8 ) ANGLE3, the angle between the two vectors, ! in the range [0,PI]. If U or V is the zero vector, ANGLE3 = PI ! is returned. ! implicit none real ( kind = 8 ) angle3 real ( kind = 8 ) dotp real ( kind = 8 ) lu real ( kind = 8 ) lv real ( kind = 8 ), parameter :: pi = 3.141592653589793D+00 real ( kind = 8 ) rtolsq real ( kind = 8 ) t real ( kind = 8 ) tol real ( kind = 8 ) u(3) real ( kind = 8 ) v(3) tol = 100.0D+00 * epsilon ( tol ) dotp = dot_product ( u(1:3), v(1:3) ) lu = dot_product ( u(1:3), u(1:3) ) lv = dot_product ( v(1:3), v(1:3) ) if ( rtolsq < lu .and. rtolsq < lv ) then t = dotp / sqrt ( lu * lv ) if ( 1.0D+00 - tol < abs ( t ) ) then t = sign ( 1.0D+00, t ) end if angle3 = acos ( t ) else angle3 = pi end if return end function areapg ( nvrt, xc, yc ) !*****************************************************************************80 ! !! AREAPG computes twice the signed area of a simple polygon. ! ! Discussion: ! ! This routine computes twice the signed area of a simple polygon, ! with vertices given in circular (counterclockwise or clockwise) order. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer ( kind = 4 ) NVRT, the number of vertices on the boundary of ! polygon. 3 <= NVRT. ! ! Input, real ( kind = 8 ) XC(1:NVRT), YC(1:NVRT), the vertex coordinates ! in counterclockwise or clockwise order. ! ! Output, real ( kind = 8 ) AREAPG, twice the signed area of the polygon, ! positive if counterclockwise. ! implicit none real ( kind = 8 ) areapg integer ( kind = 4 ) nvrt integer ( kind = 4 ) i real ( kind = 8 ) sum2 real ( kind = 8 ) xc(nvrt) real ( kind = 8 ) yc(nvrt) sum2 = xc(1) * ( yc(2) - yc(nvrt) ) + xc(nvrt) * ( yc(1) - yc(nvrt-1) ) do i = 2, nvrt-1 sum2 = sum2 + xc(i) * ( yc(i+1) - yc(i-1) ) end do areapg = sum2 return end function areatr ( xa, ya, xb, yb, xc, yc ) !*****************************************************************************80 ! !! AREATR computes twice the signed area of a triangle. ! ! Discussion: ! ! This routine computes twice the signed area of the triangle with ! vertices (XA,YA), (XB,YB), and (XC,YC) in counterclockwise or ! clockwise order. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) XA, YA, XB, YB, XC, YC, the vertex coordinates. ! ! Output, real ( kind = 8 ) AREATR, twice the signed area of triangle, ! positive if counterclockwise. ! implicit none real ( kind = 8 ) areatr real ( kind = 8 ) xa real ( kind = 8 ) xb real ( kind = 8 ) xc real ( kind = 8 ) ya real ( kind = 8 ) yb real ( kind = 8 ) yc areatr = ( xb - xa ) * ( yc - ya ) - ( xc - xa ) * ( yb - ya ) return end subroutine availf ( hdavfc, fc_num, fc_max, fc, ind, ierr ) !*****************************************************************************80 ! !! AVAILF returns the index of the next available record in the FC array. ! ! Discussion: ! ! This routine returns the index of the next available record in the ! FC array, either HDAVFC or FC_NUM+1. ! ! Modified: ! ! 07 September 2005 ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input/output, integer ( kind = 4 ) HDAVFC, head pointer of available ! records in FC. ! ! Input/output, integer ( kind = 4 ) FC_NUM, current number of records used ! in FC. ! ! Input, integer ( kind = 4 ) FC_MAX, the maximum number of records ! available in FC. ! ! Input, integer ( kind = 4 ) FC(1:7,1:*), array of face records; see ! routine DTRIS3. ! ! Output, integer ( kind = 4 ) IND, the index of available record (if FC ! not full). ! ! Output, integer ( kind = 4 ) IERR, error flag, which is zero unless an ! error occurred. ! implicit none integer ( kind = 4 ) fc(7,*) integer ( kind = 4 ) fc_num integer ( kind = 4 ) hdavfc integer ( kind = 4 ) ierr integer ( kind = 4 ) ind integer ( kind = 4 ) fc_max ierr = 0 if ( hdavfc /= 0 ) then ind = hdavfc hdavfc = -fc(1,hdavfc) else if ( fc_max <= fc_num ) then ierr = 11 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'AVAILF - Fatal error!' write ( *, '(a)' ) ' Memory requirements for array FC exceed the' write ( *, '(a,i12)' ) ' current limit of FC_MAX = ', fc_max else fc_num = fc_num + 1 ind = fc_num end if return end subroutine availk ( k, hdavfc, nfc, fc_max, fc, pos, ierr ) !*****************************************************************************80 ! !! AVAILK returns the position of the next available record in the FC array. ! ! Discussion: ! ! This routine returns the position of the next available record in ! the FC array, either HDAVFC or NFC+1. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer ( kind = 4 ) K, the number of vertices in a face. ! ! Input/output, integer ( kind = 4 ) HDAVFC, head pointer of available ! records in FC. ! ! Input/output, integer ( kind = 4 ) NFC, current number of records used ! in FC. ! ! Input, integer ( kind = 4 ) FC_MAX, maximum number of records available ! in FC. ! ! Input, integer ( kind = 4 ) FC(1:K+4,1:*), array of face records; see ! routine DTRISK. ! ! Output, integer ( kind = 4 ) POS, position of available record (if FC ! not full). ! ! Output, integer ( kind = 4 ) IERR, error flag, which is zero unless an ! error occurred. ! implicit none integer ( kind = 4 ) k integer ( kind = 4 ) fc(k+4,*) integer ( kind = 4 ) hdavfc integer ( kind = 4 ) ierr integer ( kind = 4 ) fc_max integer ( kind = 4 ) nfc integer ( kind = 4 ) pos ierr = 0 if ( hdavfc /= 0 ) then pos = hdavfc hdavfc = -fc(1,hdavfc) else if ( fc_max <= nfc ) then ierr = 22 else nfc = nfc + 1 pos = nfc end if return end subroutine baryck ( k, ind, vcl, pt, alpha, degen, mat, ipvt ) !*****************************************************************************80 ! !! BARYCK computes the barycentric coordinates of a point in KD. ! ! Discussion: ! ! This routine computes the barycentric coordinates of a point with ! respect to a simplex of K+1 vertices in K-dimensional space. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, integer ( kind = 4 ) K, dimension of points and simplex, ! ! Input, integer ( kind = 4 ) IND(1:K+1), indices in VCL of K-D vertices ! of simplex. ! ! Input, real ( kind = 8 ) VCL(1:K,1:*), K-D vertex coordinate list. ! ! Input, real ( kind = 8 ) PT(1:K), K-D point for which barycentric ! coordinates are to be computed. ! ! Output, real ( kind = 8 ) ALPHA(1:K+1), barycentric coordinates ! (if DEGEN = .FALSE.) such that ! PT = ALPHA(1)*V[IND(1)] + ... + ALPHA(K+1)*V[IND(K+1)]. ! ! Output, logical DEGEN, is TRUE if the K+1 vertices form a ! degenerate simplex. ! ! Workspace, real ( kind = 8 ) MAT(1:K,1:K), matrix used for solving ! system of linear equations. ! ! Workspace, integer IPVT(1:K-1), pivot indices. ! implicit none integer ( kind = 4 ) k real ( kind = 8 ) alpha(k+1) logical degen integer ( kind = 4 ) i integer ( kind = 4 ) ind(k+1) integer ( kind = 4 ) ipvt(k-1) integer ( kind = 4 ) j integer ( kind = 4 ) l integer ( kind = 4 ) m real ( kind = 8 ) mat(k,k) real ( kind = 8 ) pt(k) real ( kind = 8 ) tol real ( kind = 8 ) vcl(k,*) tol = 100.0D+00 * epsilon ( tol ) m = ind(k+1) do j = 1, k l = ind(j) do i = 1, k mat(i,j) = vcl(i,l) - vcl(i,m) end do end do alpha(1:k) = pt(1:k) - vcl(1:k,m) call lufac ( mat, k, k, tol, ipvt, degen ) if ( .not. degen ) then call lusol ( mat, k, k, ipvt, alpha ) alpha(k+1) = 1.0D+00 - sum ( alpha(1:k) ) end if return end subroutine baryth ( a, b, c, d, e, alpha, degen ) !*****************************************************************************80 ! !! BARYTH computes barycentric coordinates of a point in 3D. ! ! Discussion: ! ! This routine computes the barycentric coordinates of a 3D point with ! respect to the four vertices of a tetrahedron. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, real ( kind = 8 ) A(1:3), B(1:3), C(1:3), D(1:3), 4 vertices ! of tetrahedron. ! ! Input, real ( kind = 8 ) E(1:3), fifth point for which ! barycentric coordinates found ! ! Output, real ( kind = 8 ) ALPHA(1:4), the scaled barycentric coordinates ! (if DEGEN = .FALSE.) such that ! E = ( ALPHA(1) * A + ALPHA(2) * B + ALPHA(3) * C +ALPHA(4) * D ) / DET ! where DET = 6 * (volume of tetra ABCD); an ALPHA(I) may be set to 0 ! after tolerance test to indicate that E is coplanar with a face, so ! sum of ALPHA(I)/DET may not be 1; if the actual barycentric ! coordinates rather than just their signs are needed, ! modify this routine to divide ALPHA(I) by DET. ! ! Output, logical DEGEN, TRUE iff A, B, C, D are coplanar. ! implicit none real ( kind = 8 ) a(3) real ( kind = 8 ) alpha(4) real ( kind = 8 ) amax real ( kind = 8 ) b(3) real ( kind = 8 ) bmax real ( kind = 8 ) c(3) real ( kind = 8 ) cmax real ( kind = 8 ) cp1 real ( kind = 8 ) cp2 real ( kind = 8 ) cp3 real ( kind = 8 ) d(3) real ( kind = 8 ) da(3) real ( kind = 8 ) db(3) real ( kind = 8 ) dc(3) real ( kind = 8 ) de(3) logical degen real ( kind = 8 ) det real ( kind = 8 ) dmax real ( kind = 8 ) e(3) real ( kind = 8 ) ea(3) real ( kind = 8 ) eb(3) real ( kind = 8 ) ec(3) real ( kind = 8 ) emax real ( kind = 8 ) tol tol = 100.0D+00 * epsilon ( tol ) degen = .false. da(1:3) = a(1:3) - d(1:3) db(1:3) = b(1:3) - d(1:3) dc(1:3) = c(1:3) - d(1:3) amax = max ( abs ( a(1) ), abs ( a(2) ), abs ( a(3) ) ) bmax = max ( abs ( b(1) ), abs ( b(2) ), abs ( b(3) ) ) cmax = max ( abs ( c(1) ), abs ( c(2) ), abs ( c(3) ) ) dmax = max ( abs ( d(1) ), abs ( d(2) ), abs ( d(3) ) ) cp1 = db(2) * dc(3) - db(3) * dc(2) cp2 = db(3) * dc(1) - db(1) * dc(3) cp3 = db(1) * dc(2) - db(2) * dc(1) det = da(1) * cp1 + da(2) * cp2 + da(3) * cp3 if ( abs ( det ) <= 0.01D+00 * tol * max ( amax, bmax, cmax, dmax ) ) then degen = .true. return end if de(1:3) = e(1:3) - d(1:3) ea(1:3) = a(1:3) - e(1:3) eb(1:3) = b(1:3) - e(1:3) ec(1:3) = c(1:3) - e(1:3) alpha(1) = de(1) * cp1 + de(2) * cp2 + de(3) * cp3 cp1 = da(2) * de(3) - da(3) * de(2) cp2 = da(3) * de(1) - da(1) * de(3) cp3 = da(1) * de(2) - da(2) * de(1) alpha(2) = dc(1) * cp1 + dc(2) * cp2 + dc(3) * cp3 alpha(3) = db(1) * cp1 + db(2) * cp2 + db(3) * cp3 alpha(4) = ea(1) * ( eb(2) * ec(3) - eb(3) * ec(2) ) & + ea(2) * ( eb(3) * ec(1) - eb(1) * ec(3) ) & + ea(3) * ( eb(1) * ec(2) - eb(2) * ec(1) ) if ( det < 0.0D+00 ) then alpha(1) = -alpha(1) alpha(2) = -alpha(2) alpha(4) = -alpha(4) else alpha(3) = -alpha(3) end if emax = max ( abs ( e(1) ), abs ( e(2) ), abs ( e(3) ) ) if ( abs ( alpha(1) ) <= tol * max ( bmax, cmax, dmax, emax ) ) then alpha(1) = 0.0D+00 end if if ( abs ( alpha(2) ) <= tol * max ( amax, cmax, dmax, emax ) ) then alpha(2) = 0.0D+00 end if if ( abs ( alpha(3) ) <= tol * max ( amax, bmax, dmax, emax ) ) then alpha(3) = 0.0D+00 end if if ( abs ( alpha(4) ) <= tol * max ( amax, bmax, cmax, emax ) ) then alpha(4) = 0.0D+00 end if return end subroutine bcdtri ( rflag, bf_num, nbpt, nipt, sizht, fc_max, vcl, vm, nfc, & nface, ntetra, fc, ht, ierr ) !*****************************************************************************80 ! !! BCDTRI constructs a boundary-constrained Delaunay triangulation in 3D. ! ! Discussion: ! ! This routine constructs the boundary-constrained Delaunay triangulation ! of a set of 3D vertices, based on the local empty circumsphere criterion, ! by using incremental approach and local transformations. ! ! Vertices in interior of the convex hull are inserted one at a time in ! the order given by the end of the VM array. The initial tetrahedra ! created due to a new vertex are obtained by a walk through the ! triangulation until the location of the vertex is known. ! ! If there are no interior vertices specified, then one may be added if ! needed to produce a boundary-constrained triangulation. ! ! Author: ! ! Barry Joe, ! Department of Computing Science, ! University of Alberta, ! Edmonton, Alberta, Canada T6G 2H1 ! Email: barry@cs.ualberta.ca ! ! Reference: ! ! Barry Joe, ! GEOMPACK - a software package for the generation of meshes ! using geometric algorithms, ! Advances in Engineering Software, ! Volume 13, pages 325-331, 1991. ! ! Parameters: ! ! Input, logical RFLAG, is TRUE iff return immediately (with input unchanged) ! when NIPT = 0 and extra mesh vertex not removed. ! ! Input, integer ( kind = 4 ) BF_NUM, the number of boundary faces or ! triangles. ! ! Input, integer ( kind = 4 ) NBPT, the number of vertices (points) on ! boundary of convex hull. ! ! Input, integer ( kind = 4 ) NIPT, the number of vertices (points) in ! interior of convex hull. ! ! Input, integer ( kind = 4 ) SIZHT, the size of hash table HT; a good ! choice is a prime number which is about 1/8 * NFACE (or 3/2 * NPT for ! random points from the uniform distribution). ! ! Input, integer ( kind = 4 ) FC_MAX, the maximum size available for FC array. ! ! Input, real ( kind = 8 ) VCL(1:3,1:*), the vertex coordinate list. ! ! Input/output, integer ( kind = 4 ) VM(1:NPT). On input, indices of vertices of VCL ! being triangulated where NPT = NBPT + MAX ( NIPT, 1 ); VM(1:NBPT) ! are boundary points, rest are interior points; if NIPT = 0 then ! VM(NPT) is an interior point which may be added to triangulation; ! interior points are inserted in order VM(NBPT+1:NPT). ! On output, VM(NPT) is set to 0 if NIPT = 0 and extra point not needed. ! ! Input, integer ( kind = 4 ) FC(1:3,1:BF_NUM) - boundary triangles desired in triangulation; ! entries are local vertex indices 1:NBPT (indices of VM) ! ! Output, integer ( kind = 4 ) NFC, the number of positions used in FC array; ! NFC <= FC_MAX. ! ! Output, integer ( kind = 4 ) NFACE, the number of faces in triangulation; NFACE <= NFC. ! ! Output, integer ( kind = 4 ) NTETRA, the number of tetrahedra in triangulation. ! ! Output, integer ( kind = 4 ) FC(1:7,1:NFC), the array of face records which are in ! linked lists in hash table with direct chaining. Fields are: ! FC(1:3,*) - A,B,C with 1<=A