# include # include # include # include # include # include # include # include "pwl_interp_2d_scattered.h" # include "r8lib.h" /******************************************************************************/ int diaedg ( double x0, double y0, double x1, double y1, double x2, double y2, double x3, double y3 ) /******************************************************************************/ /* Purpose: DIAEDG chooses a diagonal edge. Discussion: The routine determines whether 0--2 or 1--3 is the diagonal edge that should be chosen, based on the circumcircle criterion, where (X0,Y0), (X1,Y1), (X2,Y2), (X3,Y3) are the vertices of a simple quadrilateral in counterclockwise order. Licensing: This code is distributed under the MIT license. Modified: 28 August 2003 Author: Original FORTRAN77 version by Barry Joe. C version by John Burkardt. 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. Input: double X0, Y0, X1, Y1, X2, Y2, X3, Y3, the coordinates of the vertices of a quadrilateral, given in counter clockwise order. Output: int DIAEDG, chooses a diagonal: +1, if diagonal edge 02 is chosen; -1, if diagonal edge 13 is chosen; 0, if the four vertices are cocircular. */ { double ca; double cb; double dx10; double dx12; double dx30; double dx32; double dy10; double dy12; double dy30; double dy32; double s; double tol; double tola; double tolb; int value; tol = 100.0 * DBL_EPSILON; dx10 = x1 - x0; dy10 = y1 - y0; dx12 = x1 - x2; dy12 = y1 - y2; dx30 = x3 - x0; dy30 = y3 - y0; dx32 = x3 - x2; dy32 = y3 - y2; tola = tol * fmax ( fabs ( dx10 ), fmax ( fabs ( dy10 ), fmax ( fabs ( dx30 ), fabs ( dy30 ) ) ) ); tolb = tol * fmax ( fabs ( dx12 ), fmax ( fabs ( dy12 ), fmax ( fabs ( dx32 ), fabs ( dy32 ) ) ) ); ca = dx10 * dx30 + dy10 * dy30; cb = dx12 * dx32 + dy12 * dy32; if ( tola < ca && tolb < cb ) { value = -1; } else if ( ca < -tola && cb < -tolb ) { value = 1; } else { tola = fmax ( tola, tolb ); s = ( dx10 * dy30 - dx30 * dy10 ) * cb + ( dx32 * dy12 - dx12 * dy32 ) * ca; if ( tola < s ) { value = -1; } else if ( s < -tola ) { value = 1; } else { value = 0; } } return value; } /******************************************************************************/ void i4mat_transpose_print ( int m, int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the MIT license. Modified: 31 January 2005 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, int A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { i4mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void i4mat_transpose_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: I4MAT_TRANSPOSE_PRINT_SOME prints some of an I4MAT, transposed. Discussion: An I4MAT is an MxN array of I4's, stored by (I,J) -> [I+J*M]. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 Author: John Burkardt Parameters: Input, int M, the number of rows of the matrix. M must be positive. Input, int N, the number of columns of the matrix. N must be positive. Input, int A[M*N], the matrix. Input, int ILO, JLO, IHI, JHI, designate the first row and column, and the last row and column to be printed. Input, char *TITLE, a title. */ { # define INCX 10 int i; int i2hi; int i2lo; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); if ( m <= 0 || n <= 0 ) { fprintf ( stdout, "\n" ); fprintf ( stdout, " (None)\n" ); return; } /* Print the columns of the matrix, in strips of INCX. */ for ( i2lo = ilo; i2lo <= ihi; i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( i2hi, ihi ); fprintf ( stdout, "\n" ); /* For each row I in the current range... Write the header. */ fprintf ( stdout, " Row: " ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, "%6d ", i - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { /* Print out (up to INCX) entries in column J, that lie in the current strip. */ fprintf ( stdout, "%5d: ", j - 1 ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, "%6d ", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void i4vec_heap_d ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_HEAP_D reorders an I4VEC into a descending heap. Discussion: An I4VEC is a vector of I4's. A heap is an array A with the property that, for every index J, A[J] >= A[2*J+1] and A[J] >= A[2*J+2], (as long as the indices 2*J+1 and 2*J+2 are legal). Diagram: A(0) / \ A(1) A(2) / \ / \ A(3) A(4) A(5) A(6) / \ / \ A(7) A(8) A(9) A(10) Licensing: This code is distributed under the MIT license. Modified: 30 April 1999 Author: John Burkardt Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the size of the input array. Input/output, int A[N]. On input, an unsorted array. On output, the array has been reordered into a heap. */ { int i; int ifree; int key; int m; /* Only nodes (N/2)-1 down to 0 can be "parent" nodes. */ for ( i = ( n / 2 ) - 1; 0 <= i; i-- ) { /* Copy the value out of the parent node. Position IFREE is now "open". */ key = a[i]; ifree = i; for ( ; ; ) { /* Positions 2*IFREE + 1 and 2*IFREE + 2 are the descendants of position IFREE. (One or both may not exist because they equal or exceed N.) */ m = 2 * ifree + 1; /* Does the first position exist? */ if ( n <= m ) { break; } else { /* Does the second position exist? */ if ( m + 1 < n ) { /* If both positions exist, take the larger of the two values, and update M if necessary. */ if ( a[m] < a[m+1] ) { m = m + 1; } } /* If the large descendant is larger than KEY, move it up, and update IFREE, the location of the free position, and consider the descendants of THIS position. */ if ( key < a[m] ) { a[ifree] = a[m]; ifree = m; } else { break; } } } /* When you have stopped shifting items up, return the item you pulled out back to the heap. */ a[ifree] = key; } return; } /******************************************************************************/ int i4vec_min ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_MIN returns the minimum element in an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 17 May 2003 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, int A[N], the array to be checked. Output, int I4VEC_MIN, the value of the minimum element. This is set to 0 if N <= 0. */ { int i; int value; if ( n <= 0 ) { return 0; } value = a[0]; for ( i = 1; i < n; i++ ) { if ( a[i] < value ) { value = a[i]; } } return value; } /******************************************************************************/ void i4vec_sort_heap_a ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SORT_HEAP_A ascending sorts an I4VEC using heap sort. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 30 April 1999 Author: John Burkardt Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the number of entries in the array. Input/output, int A[N]. On input, the array to be sorted; On output, the array has been sorted. */ { int n1; int temp; if ( n <= 1 ) { return; } /* 1: Put A into descending heap form. */ i4vec_heap_d ( n, a ); /* 2: Sort A. The largest object in the heap is in A[0]. Move it to position A[N-1]. */ temp = a[0]; a[0] = a[n-1]; a[n-1] = temp; /* Consider the diminished heap of size N1. */ for ( n1 = n-1; 2 <= n1; n1-- ) { /* Restore the heap structure of the initial N1 entries of A. */ i4vec_heap_d ( n1, a ); /* Take the largest object from A[0] and move it to A[N1-1]. */ temp = a[0]; a[0] = a[n1-1]; a[n1-1] = temp; } return; } /******************************************************************************/ int i4vec_sorted_unique ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SORTED_UNIQUE finds the unique elements in a sorted I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 30 June 2009 Author: John Burkardt Parameters: Input, int N, the number of elements in A. Input/output, int A[N]. On input, the sorted integer array. On output, the unique elements in A. Output, int I4VEC_SORTED_UNIQUE, the number of unique elements in A. */ { int i; int unique_num; unique_num = 0; if ( n <= 0 ) { return unique_num; } unique_num = 1; for ( i = 1; i < n; i++ ) { if ( a[i] != a[unique_num-1] ) { unique_num = unique_num + 1; a[unique_num-1] = a[i]; } } return unique_num; } /******************************************************************************/ int lrline ( double xu, double yu, double xv1, double yv1, double xv2, double yv2, double dv ) /******************************************************************************/ /* Purpose: LRLINE determines where a point lies in relation to a directed line. Discussion: LRLINE determines whether a point is to the left of, right of, or on a directed line parallel to a line through given points. Licensing: This code is distributed under the MIT license. Modified: 28 August 2003 Author: Original FORTRAN77 version by Barry Joe. C version by John Burkardt. 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, double XU, YU, XV1, YV1, XV2, YV2, are vertex coordinates; the directed line is parallel to and at signed distance DV to the left of the directed line from (XV1,YV1) to (XV2,YV2); (XU,YU) is the vertex for which the position relative to the directed line is to be determined. Input, double DV, the signed distance, positive for left. Output, int LRLINE, is +1, 0, or -1 depending on whether (XU,YU) is to the right of, on, or left of the directed line. LRLINE is 0 if the line degenerates to a point. */ { double dx; double dxu; double dy; double dyu; double t; double tol = 0.0000001; double tolabs; int value; dx = xv2 - xv1; dy = yv2 - yv1; dxu = xu - xv1; dyu = yu - yv1; tolabs = tol * fmax ( fabs ( dx ), fmax ( fabs ( dy ), fmax ( fabs ( dxu ), fmax ( fabs ( dyu ), fabs ( dv ) ) ) ) ); t = dy * dxu - dx * dyu + dv * sqrt ( dx * dx + dy * dy ); if ( tolabs < t ) { value = 1; } else if ( -tolabs <= t ) { value = 0; } else if ( t < -tolabs ) { value = -1; } return value; } /******************************************************************************/ int perm_check2 ( int n, int p[], int base ) /******************************************************************************/ /* Purpose: PERM_CHECK2 checks that a vector represents a permutation. Discussion: The routine verifies that each of the integers from 1 to N occurs among the N entries of the permutation. Set the input quantity BASE to 0, if P is a 0-based permutation, or to 1 if P is a 1-based permutation. Licensing: This code is distributed under the MIT license. Modified: 18 October 2012 Author: John Burkardt Parameters: Input, int N, the number of entries. Input, int P[N], the permutation, in standard index form. Input, int BASE, the index base. Output, int PERM_CHECK, is 1 if the array is NOT a permutation. */ { int error; int ifind; int iseek; error = 0; for ( iseek = base; iseek < base + n; iseek++ ) { error = 1; for ( ifind = 1; ifind <= n; ifind++ ) { if ( p[ifind-1] == iseek ) { error = 0; break; } } if ( error ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "PERM_CHECK2 - Fatal error!\n" ); fprintf ( stderr, " Could not find occurrence of value %d\n", iseek ); return 1; } } return 0; } /******************************************************************************/ void perm_inverse ( int n, int p[] ) /******************************************************************************/ /* Purpose: PERM_INVERSE inverts a permutation "in place". Licensing: This code is distributed under the MIT license. Modified: 28 March 2009 Author: John Burkardt Parameters: Input, int N, the number of objects being permuted. Input/output, int P[N], the permutation, in standard index form. On output, P describes the inverse permutation */ { int base; int error; int i; int i0; int i1; int i2; int is; int p_min; if ( n <= 0 ) { printf ( "\n" ); printf ( "PERM_INVERSE - Fatal error!\n" ); printf ( " Input value of N = %d\n", n ); exit ( 1 ); } /* Find the least value, and shift data so it begins at 1. */ p_min = i4vec_min ( n, p ); base = 1; for ( i = 0; i < n; i++ ) { p[i] = p[i] - p_min + base; } /* Check the permutation. */ error = perm_check2 ( n, p, base ); if ( error ) { printf ( "\n" ); printf ( "PERM_INVERSE - Fatal error!\n" ); printf ( " The input array does not represent\n" ); printf ( " a proper permutation.\n" ); exit ( 1 ); } is = 1; for ( i = 1; i <= n; i++ ) { i1 = p[i-1]; while ( i < i1 ) { i2 = p[i1-1]; p[i1-1] = -i2; i1 = i2; } is = - i4_sign ( p[i-1] ); p[i-1] = abs ( p[i-1] ) * i4_sign ( is ); } for ( i = 1; i <= n; i++ ) { i1 = -p[i-1]; if ( 0 <= i1 ) { i0 = i; for ( ; ; ) { i2 = p[i1-1]; p[i1-1] = i0; if ( i2 < 0 ) { break; } i0 = i1; i1 = i2; } } } /* Now we can restore the permutation. */ for ( i = 0; i < n; i++ ) { p[i] = p[i] + p_min - base; } return; } /******************************************************************************/ double *pwl_interp_2d_scattered_value ( int nd, double xyd[], double zd[], int t_num, int t[], int t_neighbor[], int ni, double xyi[] ) /******************************************************************************/ /* Purpose: PWL_INTERP_2D_SCATTERED_VALUE evaluates a 2d interpolant of scattered data Licensing: This code is distributed under the MIT license. Modified: 25 October 2012 Author: John Burkardt Parameters: Input, int ND, the number of data points. Input, double XYD[2*ND], the data point coordinates. Input, double ZD[ND], the data values. Input, int T_NUM, the number of triangles. Input, int T[3*T_NUM], the triangle information. Input, int T_NEIGHBOR[3*T_NUM], the triangle neighbors. Input, int NI, the number of interpolation points. Input, double XYI[2*NI], the interpolation point coordinates. Output, double PWL_INTERP_2D_SCATTERED_VALUE[NI], the interpolated values. */ { double alpha; double beta; double gamma; int edge; int i; int j; int step_num; double *zi; zi = ( double * ) malloc ( ni * sizeof ( double ) ); for ( i = 0; i < ni; i++ ) { triangulation_search_delaunay ( nd, xyd, 3, t_num, t, t_neighbor, xyi+2*i, &j, &alpha, &beta, &gamma, &edge, &step_num ); if ( j == -1 ) { zi[i] = -1.0; } zi[i] = alpha * zd[t[0+j*3]] + beta * zd[t[1+j*3]] + gamma * zd[t[2+j*3]]; } return zi; } /******************************************************************************/ int r8tris2 ( int node_num, double node_xy[], int *triangle_num, int triangle_node[], int triangle_neighbor[] ) /******************************************************************************/ /* Purpose: R8TRIS2 constructs a Delaunay triangulation of 2D vertices. Discussion: The routine constructs the Delaunay triangulation of a set of 2D vertices using an incremental approach and diagonal edge swaps. Vertices are first sorted in lexicographically increasing (X,Y) order, and then are inserted one at a time from outside the convex hull. Licensing: This code is distributed under the MIT license. Modified: 15 January 2004 Author: Original FORTRAN77 version by Barry Joe. C version by John Burkardt. 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, int NODE_NUM, the number of nodes. Input, double NODE_XY[2*NODE_NUM], the coordinates of the nodes. Output, int *TRIANGLE_NUM, the number of triangles in the triangulation; TRIANGLE_NUM is equal to 2*node_num - NB - 2, where NB is the number of boundary vertices. Output, int TRIANGLE_NODE[3*TRIANGLE_NUM], the nodes that make up each triangle. The elements are indices of NODE_XY. The vertices of the triangles are in counterclockwise order. Output, int TRIANGLE_NEIGHBOR[3*TRIANGLE_NUM], the triangle neighbor list. Positive elements are indices of TIL; negative elements are used for links of a counter clockwise linked list of boundary edges; LINK = -(3*I + J-1) where I, J = triangle, edge index; TRIANGLE_NEIGHBOR[I,J] refers to the neighbor along edge from vertex J to J+1 (mod 3). Output, int R8TRIS2, is 0 for no error. */ { double cmax; int e; int error; int i; int *indx; int j; int k; int l; int ledg; int lr; int ltri; int m; int m1; int m2; int n; int redg; int rtri; int *stack; int t; double tol; int top; stack = ( int * ) malloc ( node_num * sizeof ( int ) ); tol = 100.0 * DBL_EPSILON; /* Sort the vertices by increasing (x,y). */ indx = r82row_sort_heap_index_a ( node_num, node_xy ); r82row_permute ( node_num, indx, node_xy ); /* Make sure that the nodes are "reasonably" distinct. */ m1 = 1; for ( i = 2; i <= node_num; i++ ) { m = m1; m1 = i; k = -1; for ( j = 0; j <= 1; j++ ) { cmax = fmax ( fabs ( node_xy[2*(m-1)+j] ), fabs ( node_xy[2*(m1-1)+j] ) ); if ( tol * ( cmax + 1.0 ) < fabs ( node_xy[2*(m-1)+j] - node_xy[2*(m1-1)+j] ) ) { k = j; break; } } if ( k == -1 ) { printf ( "\n" ); printf ( "R8TRIS2 - Fatal error!\n" ); printf ( " Fails for point number I = %d\n", i ); printf ( " M = %d\n", m ); printf ( " M1 = %d\n", m1 ); printf ( " X,Y(M) = %g %g\n", node_xy[2*(m-1)+0], node_xy[2*(m-1)+1] ); printf ( " X,Y(M1) = %g %g\n", node_xy[2*(m1-1)+0], node_xy[2*(m1-1)+1] ); exit ( 1 ); } } /* Starting from nodes M1 and M2, search for a third point M that makes a "healthy" triangle (M1,M2,M) */ m1 = 1; m2 = 2; j = 3; for ( ; ; ) { if ( node_num < j ) { printf ( "\n" ); printf ( "R8TRIS2 - Fatal error!\n" ); free ( stack ); return 225; } m = j; lr = lrline ( node_xy[2*(m-1)+0], node_xy[2*(m-1)+1], node_xy[2*(m1-1)+0], node_xy[2*(m1-1)+1], node_xy[2*(m2-1)+0], node_xy[2*(m2-1)+1], 0.0 ); if ( lr != 0 ) { break; } j = j + 1; } /* Set up the triangle information for (M1,M2,M), and for any other triangles you created because nodes were collinear with M1, M2. */ *triangle_num = j - 2; if ( lr == -1 ) { triangle_node[3*0+0] = m1; triangle_node[3*0+1] = m2; triangle_node[3*0+2] = m; triangle_neighbor[3*0+2] = -3; for ( i = 2; i <= *triangle_num; i++ ) { m1 = m2; m2 = i+1; triangle_node[3*(i-1)+0] = m1; triangle_node[3*(i-1)+1] = m2; triangle_node[3*(i-1)+2] = m; triangle_neighbor[3*(i-1)+0] = -3 * i; triangle_neighbor[3*(i-1)+1] = i; triangle_neighbor[3*(i-1)+2] = i - 1; } triangle_neighbor[3*(*triangle_num-1)+0] = -3 * (*triangle_num) - 1; triangle_neighbor[3*(*triangle_num-1)+1] = -5; ledg = 2; ltri = *triangle_num; } else { triangle_node[3*0+0] = m2; triangle_node[3*0+1] = m1; triangle_node[3*0+2] = m; triangle_neighbor[3*0+0] = -4; for ( i = 2; i <= *triangle_num; i++ ) { m1 = m2; m2 = i+1; triangle_node[3*(i-1)+0] = m2; triangle_node[3*(i-1)+1] = m1; triangle_node[3*(i-1)+2] = m; triangle_neighbor[3*(i-2)+2] = i; triangle_neighbor[3*(i-1)+0] = -3 * i - 3; triangle_neighbor[3*(i-1)+1] = i - 1; } triangle_neighbor[3*(*triangle_num-1)+2] = -3 * (*triangle_num); triangle_neighbor[3*0+1] = -3 * (*triangle_num) - 2; ledg = 2; ltri = 1; } /* Insert the vertices one at a time from outside the convex hull, determine visible boundary edges, and apply diagonal edge swaps until Delaunay triangulation of vertices (so far) is obtained. */ top = 0; for ( i = j+1; i <= node_num; i++ ) { m = i; m1 = triangle_node[3*(ltri-1)+ledg-1]; if ( ledg <= 2 ) { m2 = triangle_node[3*(ltri-1)+ledg]; } else { m2 = triangle_node[3*(ltri-1)+0]; } lr = lrline ( node_xy[2*(m-1)+0], node_xy[2*(m-1)+1], node_xy[2*(m1-1)+0], node_xy[2*(m1-1)+1], node_xy[2*(m2-1)+0], node_xy[2*(m2-1)+1], 0.0 ); if ( 0 < lr ) { rtri = ltri; redg = ledg; ltri = 0; } else { l = -triangle_neighbor[3*(ltri-1)+ledg-1]; rtri = l / 3; redg = (l % 3) + 1; } vbedg ( node_xy[2*(m-1)+0], node_xy[2*(m-1)+1], node_num, node_xy, *triangle_num, triangle_node, triangle_neighbor, <ri, &ledg, &rtri, &redg ); n = *triangle_num + 1; l = -triangle_neighbor[3*(ltri-1)+ledg-1]; for ( ; ; ) { t = l / 3; e = ( l % 3 ) + 1; l = -triangle_neighbor[3*(t-1)+e-1]; m2 = triangle_node[3*(t-1)+e-1]; if ( e <= 2 ) { m1 = triangle_node[3*(t-1)+e]; } else { m1 = triangle_node[3*(t-1)+0]; } *triangle_num = *triangle_num + 1; triangle_neighbor[3*(t-1)+e-1] = *triangle_num; triangle_node[3*(*triangle_num-1)+0] = m1; triangle_node[3*(*triangle_num-1)+1] = m2; triangle_node[3*(*triangle_num-1)+2] = m; triangle_neighbor[3*(*triangle_num-1)+0] = t; triangle_neighbor[3*(*triangle_num-1)+1] = *triangle_num - 1; triangle_neighbor[3*(*triangle_num-1)+2] = *triangle_num + 1; top = top + 1; if ( node_num < top ) { printf ( "\n" ); printf ( "R8TRIS2 - Fatal error!\n" ); printf ( " Stack overflow.\n" ); free ( stack ); return 8; } stack[top-1] = *triangle_num; if ( t == rtri && e == redg ) { break; } } triangle_neighbor[3*(ltri-1)+ledg-1] = -3 * n - 1; triangle_neighbor[3*(n-1)+1] = -3 * (*triangle_num) - 2; triangle_neighbor[3*(*triangle_num-1)+2] = -l; ltri = n; ledg = 2; error = swapec ( m, &top, <ri, &ledg, node_num, node_xy, *triangle_num, triangle_node, triangle_neighbor, stack ); if ( error != 0 ) { printf ( "\n" ); printf ( "R8TRIS2 - Fatal error!\n" ); printf ( " Error return from SWAPEC.\n" ); free ( stack ); return error; } } /* Undo the sorting. */ for ( i = 0; i < 3; i++ ) { for ( j = 0; j < *triangle_num; j++ ) { triangle_node[i+j*3] = indx [ triangle_node[i+j*3] - 1 ]; } } perm_inverse ( node_num, indx ); r82row_permute ( node_num, indx, node_xy ); free ( indx ); free ( stack ); return 0; } /******************************************************************************/ int swapec ( int i, int *top, int *btri, int *bedg, int node_num, double node_xy[], int triangle_num, int triangle_node[], int triangle_neighbor[], int stack[] ) /******************************************************************************/ /* Purpose: SWAPEC swaps diagonal edges until all triangles are Delaunay. Discussion: The routine swaps diagonal edges in a 2D triangulation, based on the empty circumcircle criterion, until all triangles are Delaunay, given that I is the index of the new vertex added to the triangulation. Licensing: This code is distributed under the MIT license. Modified: 03 September 2003 Author: Original FORTRAN77 version by Barry Joe. C version by John Burkardt. 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, int I, the index of the new vertex. Input/output, int *TOP, the index of the top of the stack. On output, TOP is zero. Input/output, int *BTRI, *BEDG; on input, if positive, are the triangle and edge indices of a boundary edge whose updated indices must be recorded. On output, these may be updated because of swaps. Input, int NODE_NUM, the number of nodes. Input, double NODE_XY[2*NODE_NUM], the coordinates of the nodes. Input, int TRIANGLE_NUM, the number of triangles. Input/output, int TRIANGLE_NODE[3*TRIANGLE_NUM], the triangle incidence list. May be updated on output because of swaps. Input/output, int TRIANGLE_NEIGHBOR[3*TRIANGLE_NUM], the triangle neighbor list; negative values are used for links of the counter-clockwise linked list of boundary edges; May be updated on output because of swaps. LINK = -(3*I + J-1) where I, J = triangle, edge index. Workspace, int STACK[MAXST]; on input, entries 1 through TOP contain the indices of initial triangles (involving vertex I) put in stack; the edges opposite I should be in interior; entries TOP+1 through MAXST are used as a stack. Output, int SWAPEC, is set to 8 for abnormal return. */ { int a; int b; int c; int e; int ee; int em1; int ep1; int f; int fm1; int fp1; int l; int r; int s; int swap; int t; int tt; int u; double x; double y; /* Determine whether triangles in stack are Delaunay, and swap diagonal edge of convex quadrilateral if not. */ x = node_xy[2*(i-1)+0]; y = node_xy[2*(i-1)+1]; for ( ; ; ) { if ( *top <= 0 ) { break; } t = stack[(*top)-1]; *top = *top - 1; if ( triangle_node[3*(t-1)+0] == i ) { e = 2; b = triangle_node[3*(t-1)+2]; } else if ( triangle_node[3*(t-1)+1] == i ) { e = 3; b = triangle_node[3*(t-1)+0]; } else { e = 1; b = triangle_node[3*(t-1)+1]; } a = triangle_node[3*(t-1)+e-1]; u = triangle_neighbor[3*(t-1)+e-1]; if ( triangle_neighbor[3*(u-1)+0] == t ) { f = 1; c = triangle_node[3*(u-1)+2]; } else if ( triangle_neighbor[3*(u-1)+1] == t ) { f = 2; c = triangle_node[3*(u-1)+0]; } else { f = 3; c = triangle_node[3*(u-1)+1]; } swap = diaedg ( x, y, node_xy[2*(a-1)+0], node_xy[2*(a-1)+1], node_xy[2*(c-1)+0], node_xy[2*(c-1)+1], node_xy[2*(b-1)+0], node_xy[2*(b-1)+1] ); if ( swap == 1 ) { em1 = i4_wrap ( e - 1, 1, 3 ); ep1 = i4_wrap ( e + 1, 1, 3 ); fm1 = i4_wrap ( f - 1, 1, 3 ); fp1 = i4_wrap ( f + 1, 1, 3 ); triangle_node[3*(t-1)+ep1-1] = c; triangle_node[3*(u-1)+fp1-1] = i; r = triangle_neighbor[3*(t-1)+ep1-1]; s = triangle_neighbor[3*(u-1)+fp1-1]; triangle_neighbor[3*(t-1)+ep1-1] = u; triangle_neighbor[3*(u-1)+fp1-1] = t; triangle_neighbor[3*(t-1)+e-1] = s; triangle_neighbor[3*(u-1)+f-1] = r; if ( 0 < triangle_neighbor[3*(u-1)+fm1-1] ) { *top = *top + 1; stack[(*top)-1] = u; } if ( 0 < s ) { if ( triangle_neighbor[3*(s-1)+0] == u ) { triangle_neighbor[3*(s-1)+0] = t; } else if ( triangle_neighbor[3*(s-1)+1] == u ) { triangle_neighbor[3*(s-1)+1] = t; } else { triangle_neighbor[3*(s-1)+2] = t; } *top = *top + 1; if ( node_num < *top ) { return 8; } stack[(*top)-1] = t; } else { if ( u == *btri && fp1 == *bedg ) { *btri = t; *bedg = e; } l = - ( 3 * t + e - 1 ); tt = t; ee = em1; while ( 0 < triangle_neighbor[3*(tt-1)+ee-1] ) { tt = triangle_neighbor[3*(tt-1)+ee-1]; if ( triangle_node[3*(tt-1)+0] == a ) { ee = 3; } else if ( triangle_node[3*(tt-1)+1] == a ) { ee = 1; } else { ee = 2; } } triangle_neighbor[3*(tt-1)+ee-1] = l; } if ( 0 < r ) { if ( triangle_neighbor[3*(r-1)+0] == t ) { triangle_neighbor[3*(r-1)+0] = u; } else if ( triangle_neighbor[3*(r-1)+1] == t ) { triangle_neighbor[3*(r-1)+1] = u; } else { triangle_neighbor[3*(r-1)+2] = u; } } else { if ( t == *btri && ep1 == *bedg ) { *btri = u; *bedg = f; } l = - ( 3 * u + f - 1 ); tt = u; ee = fm1; while ( 0 < triangle_neighbor[3*(tt-1)+ee-1] ) { tt = triangle_neighbor[3*(tt-1)+ee-1]; if ( triangle_node[3*(tt-1)+0] == b ) { ee = 3; } else if ( triangle_node[3*(tt-1)+1] == b ) { ee = 1; } else { ee = 2; } } triangle_neighbor[3*(tt-1)+ee-1] = l; } } } return 0; } /******************************************************************************/ void triangulation_order3_print ( int node_num, int triangle_num, double node_xy[], int triangle_node[], int triangle_neighbor[] ) /******************************************************************************/ /* Purpose: TRIANGULATION_ORDER3_PRINT prints information defining a triangulation. Discussion: Triangulations created by R8TRIS2 include extra information encoded in the negative values of TRIANGLE_NEIGHBOR. Because some of the nodes counted in NODE_NUM may not actually be used in the triangulation, I needed to compute the true number of vertices. I added this calculation on 13 October 2001. Ernest Fasse pointed out an error in the indexing of VERTEX_LIST, which was corrected on 19 February 2004. Licensing: This code is distributed under the MIT license. Modified: 11 June 2005 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int TRIANGLE_NUM, the number of triangles. Input, double NODE_XY[2*NODE_NUM], the coordinates of the nodes. Input, int TRIANGLE_NODE[3*TRIANGLE_NUM], the nodes that make up the triangles. Input, int TRIANGLE_NEIGHBOR[3*TRIANGLE_NUM], the triangle neighbors on each side. If there is no triangle neighbor on a particular side, the value of TRIANGLE_NEIGHBOR should be negative. If the triangulation data was created by R8TRIS2, then there is more information encoded in the negative values. */ { # define DIM_NUM 2 int boundary_num; int i; int j; int k; int n1; int n2; int s; int s1; int s2; int skip; int t; int *vertex_list; int vertex_num; printf ( "\n" ); printf ( "TRIANGULATION_ORDER3_PRINT\n" ); printf ( " Information defining a triangulation.\n" ); printf ( "\n" ); printf ( " The number of nodes is %d\n", node_num ); r8mat_transpose_print ( DIM_NUM, node_num, node_xy, " Node coordinates" ); printf ( "\n" ); printf ( " The number of triangles is %d\n", triangle_num ); printf ( "\n" ); printf ( " Sets of three nodes are used as vertices of\n" ); printf ( " the triangles. For each triangle, the nodes\n" ); printf ( " are listed in counterclockwise order.\n" ); i4mat_transpose_print ( 3, triangle_num, triangle_node, " Triangle nodes" ); printf ( "\n" ); printf ( " On each side of a given triangle, there is either\n" ); printf ( " another triangle, or a piece of the convex hull.\n" ); printf ( " For each triangle, we list the indices of the three\n" ); printf ( " neighbors, or (if negative) the codes of the\n" ); printf ( " segments of the convex hull.\n" ); i4mat_transpose_print ( 3, triangle_num, triangle_neighbor, " Triangle neighbors" ); /* Determine VERTEX_NUM, the number of vertices. */ vertex_list = ( int * ) malloc ( 3 * triangle_num * sizeof ( int ) ); k = 0; for ( t = 0; t < triangle_num; t++ ) { for ( s = 0; s < 3; s++ ) { vertex_list[k] = triangle_node[s+t*3]; k = k + 1; } } i4vec_sort_heap_a ( 3*triangle_num, vertex_list ); vertex_num = i4vec_sorted_unique ( 3*triangle_num, vertex_list ); free ( vertex_list ); /* Determine the number of boundary points. */ boundary_num = 2 * vertex_num - triangle_num - 2; printf ( "\n" ); printf ( " The number of boundary points is %d\n", boundary_num ); printf ( "\n" ); printf ( " The segments that make up the convex hull can be\n" ); printf ( " determined from the negative entries of the triangle\n" ); printf ( " neighbor list.\n" ); printf ( "\n" ); printf ( " # Tri Side N1 N2\n" ); printf ( "\n" ); skip = 0; k = 0; for ( i = 0; i < triangle_num; i++ ) { for ( j = 0; j < 3; j++ ) { if ( triangle_neighbor[j+i*3] < 0 ) { s = -triangle_neighbor[j+i*3]; t = s / 3; if ( t < 1 || triangle_num < t ) { printf ( "\n" ); printf ( " Sorry, this data does not use the R8TRIS2\n" ); printf ( " convention for convex hull segments.\n" ); skip = 1; break; } s1 = ( s % 3 ) + 1; s2 = i4_wrap ( s1+1, 1, 3 ); k = k + 1; n1 = triangle_node[s1-1+(t-1)*3]; n2 = triangle_node[s2-1+(t-1)*3]; printf ( " %4d %4d %4d %4d %4d\n", k, t, s1, n1, n2 ); } } if ( skip ) { break; } } return; # undef DIM_NUM } /******************************************************************************/ void triangulation_search_delaunay ( int node_num, double node_xy[], int triangle_order, int triangle_num, int triangle_node[], int triangle_neighbor[], double p[2], int *triangle_index, double *alpha, double *beta, double *gamma, int *edge, int *step_num ) /******************************************************************************/ /* Purpose: TRIANGULATION_SEARCH_DELAUNAY searches a triangulation for a point. Discussion: The algorithm "walks" from one triangle to its neighboring triangle, and so on, until a triangle is found containing point P, or P is found to be outside the convex hull. The algorithm computes the barycentric coordinates of the point with respect to the current triangle. If all three quantities are positive, the point is contained in the triangle. If the I-th coordinate is negative, then (X,Y) lies on the far side of edge I, which is opposite from vertex I. This gives a hint as to where to search next. For a Delaunay triangulation, the search is guaranteed to terminate. For other triangulations, a cycle may occur. Note the surprising fact that, even for a Delaunay triangulation of a set of nodes, the nearest point to (X,Y) need not be one of the vertices of the triangle containing (X,Y). The code can be called for triangulations of any order, but only the first three nodes in each triangle are considered. Thus, if higher order triangles are used, and the extra nodes are intended to give the triangle a polygonal shape, these will have no effect, and the results obtained here might be misleading. Licensing: This code is distributed under the MIT license. Modified: 18 August 2009 Author: Original FORTRAN77 version by Barry Joe. C version by John Burkardt. 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, int NODE_NUM, the number of nodes. Input, double NODE_XY[2*NODE_NUM], the coordinates of the nodes. Input, int TRIANGLE_ORDER, the order of the triangles. Input, int TRIANGLE_NUM, the number of triangles in the triangulation. Input, int TRIANGLE_NODE[TRIANGLE_ORDER*TRIANGLE_NUM], the nodes of each triangle. Input, int TRIANGLE_NEIGHBOR[3*TRIANGLE_NUM], the triangle neighbor list. Input, double P[2], the coordinates of a point. Output, int *TRIANGLE_INDEX, the index of the triangle where the search ended. If a cycle occurred, then TRIANGLE_INDEX = -1. Output, double *ALPHA, *BETA, *GAMMA, the barycentric coordinates of the point relative to triangle TRIANGLE_INDEX. Output, int *EDGE, indicates the position of the point (X,Y) in triangle TRIANGLE: 0, the interior or boundary of the triangle; -1, outside the convex hull of the triangulation, past edge 1; -2, outside the convex hull of the triangulation, past edge 2; -3, outside the convex hull of the triangulation, past edge 3. Output, int *STEP_NUM, the number of steps. */ { int a; int b; int c; double det; double dxp; double dxa; double dxb; double dyp; double dya; double dyb; static int triangle_index_save = -1; *step_num = - 1; *edge = 0; if ( triangle_index_save < 0 || triangle_num <= triangle_index_save ) { *triangle_index = ( triangle_num + 1 ) / 2 - 1; } else { *triangle_index = triangle_index_save; } for ( ; ; ) { *step_num = *step_num + 1; if ( triangle_num < *step_num ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRIANGULATION_SEARCH_DELAUNAY - Fatal error!\n" ); fprintf ( stderr, " The algorithm seems to be cycling.\n" ); fprintf ( stderr, " Current triangle is %d\n", *triangle_index ); exit ( 1 ); } /* Get the vertices of triangle TRIANGLE. */ a = triangle_node[0+(*triangle_index)*triangle_order]; b = triangle_node[1+(*triangle_index)*triangle_order]; c = triangle_node[2+(*triangle_index)*triangle_order]; /* Using vertex C as a base, compute the distances to vertices A and B, and the point (X,Y). */ dxa = node_xy[0+a*2] - node_xy[0+c*2]; dya = node_xy[1+a*2] - node_xy[1+c*2]; dxb = node_xy[0+b*2] - node_xy[0+c*2]; dyb = node_xy[1+b*2] - node_xy[1+c*2]; dxp = p[0] - node_xy[0+c*2]; dyp = p[1] - node_xy[1+c*2]; det = dxa * dyb - dya * dxb; /* Compute the barycentric coordinates of the point (X,Y) with respect to this triangle. */ *alpha = ( dxp * dyb - dyp * dxb ) / det; *beta = ( dxa * dyp - dya * dxp ) / det; *gamma = 1.0 - *alpha - *beta; /* If the barycentric coordinates are all positive, then the point is inside the triangle and we're done. */ if ( 0.0 <= *alpha && 0.0 <= *beta && 0.0 <= *gamma ) { break; } /* At least one barycentric coordinate is negative. If there is a negative barycentric coordinate for which there exists an opposing triangle neighbor closer to the point, move to that triangle. (Two coordinates could be negative, in which case we could go for the most negative one, or the most negative one normalized by the actual distance it represents). */ if ( *alpha < 0.0 && 0 <= triangle_neighbor[1+(*triangle_index)*3] ) { *triangle_index = triangle_neighbor[1+(*triangle_index)*3]; continue; } else if ( *beta < 0.0 && 0 <= triangle_neighbor[2+(*triangle_index)*3] ) { *triangle_index = triangle_neighbor[2+(*triangle_index)*3]; continue; } else if ( *gamma < 0.0 && 0 <= triangle_neighbor[0+(*triangle_index)*3] ) { *triangle_index = triangle_neighbor[0+(*triangle_index)*3]; continue; } /* All negative barycentric coordinates correspond to vertices opposite sides on the convex hull. Note the edge and exit. */ if ( *alpha < 0.0 ) { *edge = -2; break; } else if ( *beta < 0.0 ) { *edge = -3; break; } else if ( *gamma < 0.0 ) { *edge = -1; break; } else { printf ( "\n" ); printf ( "TRIANGULATION_ORDER3_SEARCH - Fatal error!\n" ); printf ( " The algorithm seems to have reached a dead end\n" ); printf ( " after %d steps.\n", *step_num ); printf ( " Currently, TRIANGLE_INDEX = %d\n", *triangle_index ); *triangle_index = -1; *edge = -1; return; } } triangle_index_save = *triangle_index; return; } /******************************************************************************/ void vbedg ( double x, double y, int node_num, double node_xy[], int triangle_num, int triangle_node[], int triangle_neighbor[], int *ltri, int *ledg, int *rtri, int *redg ) /******************************************************************************/ /* Purpose: VBEDG determines which boundary edges are visible to a point. Discussion: The point (X,Y) is assumed to be outside the convex hull of the region covered by the 2D triangulation. Licensing: This code is distributed under the MIT license. Modified: 30 September 2008 Author: Original FORTRAN77 version by Barry Joe. C version by John Burkardt. 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, double X, Y, the coordinates of a point outside the convex hull of the current triangulation. Input, int NODE_NUM, the number of nodes. Input, double NODE_XY[2*NODE_NUM], the coordinates of the nodes. Input, int TRIANGLE_NUM, the number of triangles. Input, int TRIANGLE_NODE[3*TRIANGLE_NUM], the triangle incidence list. Input, int TRIANGLE_NEIGHBOR[3*TRIANGLE_NUM], the triangle neighbor list; negative values are used for links of a counter clockwise linked list of boundary edges; LINK = -(3*I + J-1) where I, J = triangle, edge index. Input/output, int *LTRI, *LEDG. If LTRI != 0 then these values are assumed to be already computed and are not changed, else they are updated. On output, LTRI is the index of boundary triangle to the left of the leftmost boundary triangle visible from (X,Y), and LEDG is the boundary edge of triangle LTRI to the left of the leftmost boundary edge visible from (X,Y). 1 <= LEDG <= 3. Input/output, int *RTRI. On input, the index of the boundary triangle to begin the search at. On output, the index of the rightmost boundary triangle visible from (X,Y). Input/output, int *REDG, the edge of triangle RTRI that is visible from (X,Y). 1 <= REDG <= 3. */ { int a; double ax; double ay; int b; double bx; double by; int done; int e; int l; int lr; int t; /* Find the rightmost visible boundary edge using links, then possibly leftmost visible boundary edge using triangle neighbor information. */ if ( *ltri == 0 ) { done = 0; *ltri = *rtri; *ledg = *redg; } else { done = 1; } for ( ; ; ) { l = -triangle_neighbor[3*((*rtri)-1)+(*redg)-1]; t = l / 3; e = 1 + l % 3; a = triangle_node[3*(t-1)+e-1]; if ( e <= 2 ) { b = triangle_node[3*(t-1)+e]; } else { b = triangle_node[3*(t-1)+0]; } ax = node_xy[2*(a-1)+0]; ay = node_xy[2*(a-1)+1]; bx = node_xy[2*(b-1)+0]; by = node_xy[2*(b-1)+1]; lr = lrline ( x, y, ax, ay, bx, by, 0.0 ); if ( lr <= 0 ) { break; } *rtri = t; *redg = e; } if ( done ) { return; } t = *ltri; e = *ledg; for ( ; ; ) { b = triangle_node[3*(t-1)+e-1]; e = i4_wrap ( e-1, 1, 3 ); while ( 0 < triangle_neighbor[3*(t-1)+e-1] ) { t = triangle_neighbor[3*(t-1)+e-1]; if ( triangle_node[3*(t-1)+0] == b ) { e = 3; } else if ( triangle_node[3*(t-1)+1] == b ) { e = 1; } else { e = 2; } } a = triangle_node[3*(t-1)+e-1]; ax = node_xy[2*(a-1)+0]; ay = node_xy[2*(a-1)+1]; bx = node_xy[2*(b-1)+0]; by = node_xy[2*(b-1)+1]; lr = lrline ( x, y, ax, ay, bx, by, 0.0 ); if ( lr <= 0 ) { break; } } *ltri = t; *ledg = e; return; }