# include # include # include # include # include # include "tet_mesh.h" /******************************************************************************/ int i4_max ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MAX returns the maximum of two I4's. Licensing: This code is distributed under the MIT license. Modified: 13 October 1998 Author: John Burkardt Parameters: Input, int I1, I2, are two integers to be compared. Output, int I4_MAX, the larger of I1 and I2. */ { int value; if ( i2 < i1 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ int i4_min ( int i1, int i2 ) /******************************************************************************/ /* Purpose: I4_MIN returns the minimum of two I4's. Licensing: This code is distributed under the MIT license. Modified: 13 October 1998 Author: John Burkardt Parameters: Input, int I1, I2, two integers to be compared. Output, int I4_MIN, the smaller of I1 and I2. */ { int value; if ( i1 < i2 ) { value = i1; } else { value = i2; } return value; } /******************************************************************************/ void i4_swap ( int *i, int *j ) /******************************************************************************/ /* Purpose: I4_SWAP switches two I4's. Licensing: This code is distributed under the MIT license. Modified: 07 January 2002 Author: John Burkardt Parameters: Input/output, int *I, *J. On output, the values of I and J have been interchanged. */ { int k; k = *i; *i = *j; *j = k; return; } /******************************************************************************/ int i4_uniform_ab ( int a, int b, int *seed ) /******************************************************************************/ /* Purpose: I4_UNIFORM_AB returns a scaled pseudorandom I4 between A and B. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the MIT license. Modified: 24 May 2012 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Parameters: Input, int A, B, the limits of the interval. Input/output, int *SEED, the "seed" value, which should NOT be 0. On output, SEED has been updated. Output, int I4_UNIFORM_AB, a number between A and B. */ { int c; const int i4_huge = 2147483647; int k; float r; int value; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_UNIFORM_AB - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } /* Guaranteee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( float ) ( *seed ) * 4.656612875E-10; /* Scale R to lie between A-0.5 and B+0.5. */ r = ( 1.0 - r ) * ( ( float ) ( a ) - 0.5 ) + r * ( ( float ) ( b ) + 0.5 ); /* Round R to the nearest integer. */ value = round ( r ); /* Guarantee that A <= VALUE <= B. */ if ( value < a ) { value = a; } if ( b < value ) { value = b; } return value; } /******************************************************************************/ int i4col_compare ( int m, int n, int a[], int i, int j ) /******************************************************************************/ /* Purpose: I4COL_COMPARE compares columns I and J of an I4COL. Example: Input: M = 3, N = 4, I = 2, J = 4 A = ( 1 2 3 4 5 6 7 8 9 10 11 12 ) Output: I4COL_COMPARE = -1 Licensing: This code is distributed under the MIT license. Modified: 21 March 2009 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, int A[M*N], an array of N columns of vectors of length M. Input, int I, J, the columns to be compared. I and J must be between 1 and N. Output, int I4COL_COMPARE, the results of the comparison: -1, column I < column J, 0, column I = column J, +1, column J < column I. */ { int k; /* Check. */ if ( i < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4COL_COMPARE - Fatal error!\n" ); fprintf ( stderr, " Column index I = %d is less than 1.\n", i ); exit ( 1 ); } if ( n < i ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4COL_COMPARE - Fatal error!\n" ); fprintf ( stderr, " N = %d is less than column index I = %d.\n", n, i ); exit ( 1 ); } if ( j < 1 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4COL_COMPARE - Fatal error!\n" ); fprintf ( stderr, " Column index J = %d is less than 1.\n", j ); exit ( 1 ); } if ( n < j ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4COL_COMPARE - Fatal error!\n" ); fprintf ( stderr, " N = %d is less than column index J = %d.\n", n, j ); exit ( 1 ); } if ( i == j ) { return 0; } k = 1; while ( k <= m ) { if ( a[k-1+(i-1)*m] < a[k-1+(j-1)*m] ) { return (-1); } else if ( a[k-1+(j-1)*m] < a[k-1+(i-1)*m] ) { return 1; } k = k + 1; } return 0; }/******************************************************************************/ void i4col_sort_a ( int m, int n, int a[] ) /******************************************************************************/ /* Purpose: I4COL_SORT_A ascending sorts the columns of an I4COL. Discussion: In lexicographic order, the statement "X < Y", applied to two vectors X and Y of length M, means that there is some index I, with 1 <= I <= M, with the property that X(J) = Y(J) for J < I, and X(I) < Y(I). In other words, X is less than Y if, at the first index where they differ, the X value is less than the Y value. Licensing: This code is distributed under the MIT license. Modified: 30 June 2009 Author: John Burkardt Parameters: Input, int M, the number of rows of A. Input, int N, the number of columns of A. Input/output, int A[M*N]. On input, the array of N columns of M vectors; On output, the columns of A have been sorted in ascending lexicographic order. */ { int i; int indx; int isgn; int j; /* Initialize. */ i = 0; indx = 0; isgn = 0; j = 0; /* Call the external heap sorter. */ for ( ; ; ) { sort_heap_external ( n, &indx, &i, &j, isgn ); /* Interchange the I and J objects. */ if ( 0 < indx ) { i4col_swap ( m, n, a, i, j ); } /* Compare the I and J objects. */ else if ( indx < 0 ) { isgn = i4col_compare ( m, n, a, i, j ); } else if ( indx == 0 ) { break; } } return; } /******************************************************************************/ void i4col_sort2_a ( int m, int n, int a[] ) /******************************************************************************/ /* Purpose: I4COL_SORT2_A ascending sorts the elements of each column of an I4COL. Licensing: This code is distributed under the MIT license. Modified: 05 June 2010 Author: John Burkardt Parameters: Input, int M, the number of rows of A. Input, int N, the number of columns of A, and the length of a vector of data. Input/output, int A[M*N]. On input, the array of N columns of M vectors. On output, the elements of each column of A have been sorted in ascending order. */ { int col; int i; int indx; int isgn; int j; int temp; if ( m <= 1 ) { return; } if ( n <= 0 ) { return; } /* Initialize. */ for ( col = 0; col < n; col++ ) { i = 0; indx = 0; isgn = 0; j = 0; /* Call the external heap sorter. */ for ( ; ; ) { sort_heap_external ( m, &indx, &i, &j, isgn ); /* Interchange the I and J objects. */ if ( 0 < indx ) { temp = a[i-1+col*m]; a[i-1+col*m] = a[j-1+col*m]; a[j-1+col*m] = temp; } /* Compare the I and J objects. */ else if ( indx < 0 ) { if ( a[j-1+col*m] < a[i-1+col*m] ) { isgn = +1; } else { isgn = -1; } } else if ( indx == 0 ) { break; } } } return; } /******************************************************************************/ int i4col_sorted_unique_count ( int m, int n, int a[] ) /******************************************************************************/ /* Purpose: I4COL_SORTED_UNIQUE_COUNT counts unique elements in an I4COL. Discussion: The columns of the array may be ascending or descending sorted. Licensing: This code is distributed under the MIT license. Modified: 30 June 2009 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, int A[M*N], a sorted array, containing N columns of data. Output, int I4COL_SORTED_UNIQUE_COUNT, the number of unique columns. */ { int i; int j1; int j2; int unique_num; if ( n <= 0 ) { unique_num = 0; return unique_num; } unique_num = 1; j1 = 0; for ( j2 = 1; j2 < n; j2++ ) { for ( i = 0; i < m; i++ ) { if ( a[i+j1*m] != a[i+j2*m] ) { unique_num = unique_num + 1; j1 = j2; break; } } } return unique_num; } /******************************************************************************/ void i4col_swap ( int m, int n, int a[], int icol1, int icol2 ) /******************************************************************************/ /* Purpose: I4COL_SWAP swaps two columns of an I4COL. Discussion: The two dimensional information is stored as a one dimensional array, by columns. The row indices are 1 based, NOT 0 based. However, a preprocessor variable, called OFFSET, can be reset from 1 to 0 if you wish to use 0-based indices. Licensing: This code is distributed under the MIT license. Modified: 30 June 2009 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input/output, int A[M*N], an array of data. Input, int ICOL1, ICOL2, the two columns to swap. These indices should be between 1 and N. */ { # define OFFSET 1 int i; int t; /* Check. */ if ( icol1 - OFFSET < 0 || n-1 < icol1 - OFFSET ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4COL_SWAP - Fatal error!\n" ); fprintf ( stderr, " ICOL1 is out of range.\n" ); exit ( 1 ); } if ( icol2 - OFFSET < 0 || n-1 < icol2 - OFFSET ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4COL_SWAP - Fatal error!\n" ); fprintf ( stderr, " ICOL2 is out of range.\n" ); exit ( 1 ); } if ( icol1 == icol2 ) { return; } for ( i = 0; i < m; i++ ) { t = a[i+(icol1-OFFSET)*m]; a[i+(icol1-OFFSET)*m] = a[i+(icol2-OFFSET)*m]; a[i+(icol2-OFFSET)*m] = t; } return; # undef OFFSET } /******************************************************************************/ void i4i4_sort_a ( int i1, int i2, int *j1, int *j2 ) /******************************************************************************/ /* Purpose: I4I4_SORT_A ascending sorts a pair of I4's. Discussion: The program allows the reasonable call: i4i4_sort_a ( i1, i2, &i1, &i2 ); and this will return the reasonable result. Licensing: This code is distributed under the MIT license. Modified: 11 October 2005 Author: John Burkardt Parameters: Input, int I1, I2, the values to sort. Output, int J1, J2, the sorted values. */ { int k1; int k2; /* Copy arguments, so that the user can make "reasonable" calls like: i4i4_sort_a ( i1, i2, &i1, &i2 ); */ k1 = i1; k2 = i2; *j1 = i4_min ( k1, k2 ); *j2 = i4_max ( k1, k2 ); return; } /******************************************************************************/ void i4i4i4_sort_a ( int i1, int i2, int i3, int *j1, int *j2, int *j3 ) /******************************************************************************/ /* Purpose: I4I4I4_SORT_A ascending sorts a triple of I4's. Discussion: The program allows the reasonable call: i4i4i4_sort_a ( i1, i2, i3, &i1, &i2, &i3 ); and this will return the reasonable result. Licensing: This code is distributed under the MIT license. Modified: 12 October 2005 Author: John Burkardt Parameters: Input, int I1, I2, I3, the values to sort. Output, int *J1, *J2, *J3, the sorted values. */ { int k1; int k2; int k3; /* Copy arguments, so that the user can make "reasonable" calls like: i4i4i4_sort_a ( i1, i2, i3, &i1, &i2, &i3 ); */ k1 = i1; k2 = i2; k3 = i3; *j1 = i4_min ( i4_min ( k1, k2 ), i4_min ( k2, k3 ) ); *j2 = i4_min ( i4_max ( k1, k2 ), i4_min ( i4_max ( k2, k3 ), i4_max ( k3, k1 ) ) ); *j3 = i4_max ( i4_max ( k1, k2 ), i4_max ( k2, k3 ) ); return; } /******************************************************************************/ 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; if ( m < i2hi ) { i2hi = m; } if ( ihi < i2hi ) { 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 = jlo; if ( j2lo < 1 ) { j2lo = 1; } j2hi = jhi; if ( n < jhi ) { j2hi = 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_print ( int n, int a[], char *title ) /******************************************************************************/ /* Purpose: I4VEC_PRINT prints an I4VEC. Discussion: An I4VEC is a vector of I4's. Licensing: This code is distributed under the MIT license. Modified: 14 November 2003 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, int A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %6d: %8d\n", i, a[i] ); } return; } /******************************************************************************/ int i4vec_sum ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_SUM sums the entries of an I4VEC. Discussion: An I4VEC is a vector of I4's. Example: Input: A = ( 1, 2, 3, 4 ) Output: I4VEC_SUM = 10 Licensing: This code is distributed under the MIT license. Modified: 26 May 1999 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, int A[N], the vector to be summed. Output, int I4VEC_SUM, the sum of the entries of A. */ { int i; int sum; sum = 0; for ( i = 0; i < n; i++ ) { sum = sum + a[i]; } return sum; } /******************************************************************************/ void i4vec_zero ( int n, int a[] ) /******************************************************************************/ /* Purpose: I4VEC_ZERO zeroes an I4VEC. Licensing: This code is distributed under the MIT license. Modified: 01 August 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, int A[N], a vector of zeroes. */ { int i; for ( i = 0; i < n; i++ ) { a[i] = 0; } return; } /******************************************************************************/ void mesh_base_one ( int node_num, int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: MESH_BASE_ONE ensures that the element definition is one-based. Discussion: The ELEMENT_NODE array contains nodes indices that form elements. The convention for node indexing might start at 0 or at 1. If this function detects 0-based indexing, it converts to 1-based. Licensing: This code is distributed under the MIT license. Modified: 18 October 2014 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int ELEMENT_ORDER, the order of the elements. Input, int ELEMENT_NUM, the number of elements. Input/output, int ELEMENT_NODE[ELEMENT_ORDER*ELEMENT_NUM], the element definitions. */ { int element; const int i4_huge = 2147483647; int node; int node_max; int node_min; int order; node_min = + i4_huge; node_max = - i4_huge; for ( element = 0; element < element_num; element++ ) { for ( order = 0; order < element_order; order++ ) { node = element_node[order+element*element_order]; if ( node < node_min ) { node_min = node; } if ( node_max < node ) { node_max = node; } } } if ( node_min == 0 && node_max == node_num - 1 ) { printf ( "\n" ); printf ( "MESH_BASE_ONE:\n" ); printf ( " The element indexing appears to be 0-based!\n" ); printf ( " This will be converted to 1-based.\n" ); for ( element = 0; element < element_num; element++ ) { for ( order = 0; order < element_order; order++ ) { element_node[order+element*element_order] = element_node[order+element*element_order] + 1; } } } else if ( node_min == 1 && node_max == node_num ) { printf ( "\n" ); printf ( "MESH_BASE_ONE:\n" ); printf ( " The element indexing appears to be 1-based!\n" ); printf ( " No conversion is necessary.\n" ); } else { printf ( "\n" ); printf ( "MESH_BASE_ONE - Warning!\n" ); printf ( " The element indexing is not of a recognized type.\n" ); printf ( " NODE_MIN = %d\n", node_min ); printf ( " NODE_MAX = %d\n", node_max ); printf ( " NODE_NUM = %d\n", node_num ); } return; } /******************************************************************************/ void mesh_base_zero ( int node_num, int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: MESH_BASE_ZERO ensures that the element definition is zero-based. Discussion: The ELEMENT_NODE array contains nodes indices that form elements. The convention for node indexing might start at 0 or at 1. Since a C program will naturally assume a 0-based indexing, it is necessary to check a given element definition and, if it is actually 1-based, to convert it. This function attempts to detect 1-based node indexing and correct it. Licensing: This code is distributed under the MIT license. Modified: 18 October 2014 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int ELEMENT_ORDER, the order of the elements. Input, int ELEMENT_NUM, the number of elements. Input/output, int ELEMENT_NODE[ELEMENT_ORDER*ELEMENT_NUM], the element definitions. */ { int element; const int i4_huge = 2147483647; int node; int node_max; int node_min; int order; node_min = + i4_huge; node_max = - i4_huge; for ( element = 0; element < element_num; element++ ) { for ( order = 0; order < element_order; order++ ) { node = element_node[order+element*element_order]; if ( node < node_min ) { node_min = node; } if ( node_max < node ) { node_max = node; } } } if ( node_min == 1 && node_max == node_num ) { printf ( "\n" ); printf ( "MESH_BASE_ZERO:\n" ); printf ( " The element indexing appears to be 1-based!\n" ); printf ( " This will be converted to 0-based.\n" ); for ( element = 0; element < element_num; element++ ) { for ( order = 0; order < element_order; order++ ) { element_node[order+element*element_order] = element_node[order+element*element_order] - 1; } } } else if ( node_min == 0 && node_max == node_num - 1 ) { printf ( "\n" ); printf ( "MESH_BASE_ZERO:\n" ); printf ( " The element indexing appears to be 0-based!\n" ); printf ( " No conversion is necessary.\n" ); } else { printf ( "\n" ); printf ( "MESH_BASE_ZERO - Warning!\n" ); printf ( " The element indexing is not of a recognized type.\n" ); printf ( " NODE_MIN = %d\n", node_min ); printf ( " NODE_MAX = %d\n", node_max ); printf ( " NODE_NUM = %d\n", node_num ); } return; } /******************************************************************************/ int r4_nint ( float x ) /******************************************************************************/ /* Purpose: R4_NINT returns the nearest integer to an R4. Example: X R4_NINT 1.3 1 1.4 1 1.5 1 or 2 1.6 2 0.0 0 -0.7 -1 -1.1 -1 -1.6 -2 Licensing: This code is distributed under the MIT license. Modified: 14 November 2006 Author: John Burkardt Parameters: Input, float X, the value. Output, int R4_NINT, the nearest integer to X. */ { int value; if ( x < 0.0 ) { value = - ( int ) ( fabs ( x ) + 0.5 ); } else { value = ( int ) ( fabs ( x ) + 0.5 ); } return value; } /******************************************************************************/ void r8_swap ( double *x, double *y ) /******************************************************************************/ /* Purpose: R8_SWAP switches two R8's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2004 Author: John Burkardt Parameters: Input/output, double *X, *Y. On output, the values of X and Y have been interchanged. */ { double z; z = *x; *x = *y; *y = z; return; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: R8_UNIFORM_01 returns a pseudorandom R8 scaled to [0,1]. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) r8_uniform_01 = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. If the initial seed is 12345, then the first three computations are Input Output R8_UNIFORM_01 SEED SEED 12345 207482415 0.096616 207482415 1790989824 0.833995 1790989824 2035175616 0.947702 Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation edited by Jerry Banks, Wiley Interscience, page 95, 1998. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. P A Lewis, A S Goodman, J M Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input/output, int *SEED, the "seed" value. Normally, this value should not be 0. On output, SEED has been updated. Output, double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { const int i4_huge = 2147483647; int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( ( double ) ( *seed ) ) * 4.656612875E-10; return r; } /******************************************************************************/ double r8mat_det_4d ( double a[4*4] ) /******************************************************************************/ /* Purpose: R8MAT_DET_4D computes the determinant of a 4 by 4 R8MAT. Discussion: The two dimensional array is stored as a one dimensional vector, by COLUMNS. Licensing: This code is distributed under the MIT license. Modified: 10 September 2003 Author: John Burkardt Parameters: Input, double A[4*4], the matrix whose determinant is desired. Output, double R8MAT_DET_4D, the determinant of the matrix. */ { double det; det = a[0+0*4] * ( a[1+1*4] * ( a[2+2*4] * a[3+3*4] - a[2+3*4] * a[3+2*4] ) - a[1+2*4] * ( a[2+1*4] * a[3+3*4] - a[2+3*4] * a[3+1*4] ) + a[1+3*4] * ( a[2+1*4] * a[3+2*4] - a[2+2*4] * a[3+1*4] ) ) - a[0+1*4] * ( a[1+0*4] * ( a[2+2*4] * a[3+3*4] - a[2+3*4] * a[3+2*4] ) - a[1+2*4] * ( a[2+0*4] * a[3+3*4] - a[2+3*4] * a[3+0*4] ) + a[1+3*4] * ( a[2+0*4] * a[3+2*4] - a[2+2*4] * a[3+0*4] ) ) + a[0+2*4] * ( a[1+0*4] * ( a[2+1*4] * a[3+3*4] - a[2+3*4] * a[3+1*4] ) - a[1+1*4] * ( a[2+0*4] * a[3+3*4] - a[2+3*4] * a[3+0*4] ) + a[1+3*4] * ( a[2+0*4] * a[3+1*4] - a[2+1*4] * a[3+0*4] ) ) - a[0+3*4] * ( a[1+0*4] * ( a[2+1*4] * a[3+2*4] - a[2+2*4] * a[3+1*4] ) - a[1+1*4] * ( a[2+0*4] * a[3+2*4] - a[2+2*4] * a[3+0*4] ) + a[1+2*4] * ( a[2+0*4] * a[3+1*4] - a[2+1*4] * a[3+0*4] ) ); return det; } /******************************************************************************/ double *r8mat_mv_new ( int m, int n, double a[], double x[] ) /******************************************************************************/ /* Purpose: R8MAT_MV_NEW multiplies a matrix times a vector. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. For this routine, the result is returned as the function value. Licensing: This code is distributed under the MIT license. Modified: 11 April 2007 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns of the matrix. Input, double A[M,N], the M by N matrix. Input, double X[N], the vector to be multiplied by A. Output, double R8MAT_MV_NEW[M], the product A*X. */ { int i; int j; double *y; y = ( double * ) malloc ( m * sizeof ( double ) ); for ( i = 0; i < m; i++ ) { y[i] = 0.0; for ( j = 0; j < n; j++ ) { y[i] = y[i] + a[i+j*m] * x[j]; } } return y; } /******************************************************************************/ void r8mat_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT prints an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*M] Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, the number of rows in A. Input, int N, the number of columns in A. Input, double A[M*N], the M by N matrix. Input, char *TITLE, a title. */ { r8mat_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_PRINT_SOME prints some of an R8MAT. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 26 June 2013 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, double 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 5 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 5. */ for ( j2lo = jlo; j2lo <= jhi; j2lo = j2lo + INCX ) { j2hi = j2lo + INCX - 1; if ( n < j2hi ) { j2hi = n; } if ( jhi < j2hi ) { j2hi = jhi; } fprintf ( stdout, "\n" ); /* For each column J in the current range... Write the header. */ fprintf ( stdout, " Col: "); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %7d ", j - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Row\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ if ( 1 < ilo ) { i2lo = ilo; } else { i2lo = 1; } if ( m < ihi ) { i2hi = m; } else { i2hi = ihi; } for ( i = i2lo; i <= i2hi; i++ ) { /* Print out (up to) 5 entries in row I, that lie in the current strip. */ fprintf ( stdout, "%5d:", i - 1 ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, " %14g", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ int r8mat_solve ( int n, int rhs_num, double a[] ) /******************************************************************************/ /* Purpose: R8MAT_SOLVE uses Gauss-Jordan elimination to solve an N by N linear system. Discussion: A R8MAT is a doubly dimensioned array of double precision values, which may be stored as a vector in column-major order. Entry A(I,J) is stored as A[I+J*N] Licensing: This code is distributed under the MIT license. Modified: 29 August 2003 Author: John Burkardt Parameters: Input, int N, the order of the matrix. Input, int RHS_NUM, the number of right hand sides. RHS_NUM must be at least 0. Input/output, double A[N*(N+RHS_NUM)], contains in rows and columns 1 to N the coefficient matrix, and in columns N+1 through N+RHS_NUM, the right hand sides. On output, the coefficient matrix area has been destroyed, while the right hand sides have been overwritten with the corresponding solutions. Output, int R8MAT_SOLVE, singularity flag. 0, the matrix was not singular, the solutions were computed; J, factorization failed on step J, and the solutions could not be computed. */ { double apivot; double factor; int i; int ipivot; int j; int k; double temp; for ( j = 0; j < n; j++ ) { /* Choose a pivot row. */ ipivot = j; apivot = a[j+j*n]; for ( i = j; i < n; i++ ) { if ( fabs ( apivot ) < fabs ( a[i+j*n] ) ) { apivot = a[i+j*n]; ipivot = i; } } if ( apivot == 0.0 ) { return j; } /* Interchange. */ for ( i = 0; i < n + rhs_num; i++ ) { temp = a[ipivot+i*n]; a[ipivot+i*n] = a[j+i*n]; a[j+i*n] = temp; } /* A(J,J) becomes 1. */ a[j+j*n] = 1.0; for ( k = j; k < n + rhs_num; k++ ) { a[j+k*n] = a[j+k*n] / apivot; } /* A(I,J) becomes 0. */ for ( i = 0; i < n; i++ ) { if ( i != j ) { factor = a[i+j*n]; a[i+j*n] = 0.0; for ( k = j; k < n + rhs_num; k++ ) { a[i+k*n] = a[i+k*n] - factor * a[j+k*n]; } } } } return 0; } /******************************************************************************/ void r8mat_transpose_print ( int m, int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, char *TITLE, a title. */ { r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ); return; } /******************************************************************************/ void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ) /******************************************************************************/ /* Purpose: R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. Licensing: This code is distributed under the MIT license. Modified: 10 September 2013 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A[M*N], an M by N matrix to be printed. Input, int ILO, JLO, the first row and column to print. Input, int IHI, JHI, the last row and column to print. Input, char *TITLE, a title. */ { # define INCX 5 int i; int i2; int i2hi; int i2lo; int i2lo_hi; int i2lo_lo; int inc; 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; } if ( ilo < 1 ) { i2lo_lo = 1; } else { i2lo_lo = ilo; } if ( ihi < m ) { i2lo_hi = m; } else { i2lo_hi = ihi; } for ( i2lo = i2lo_lo; i2lo <= i2lo_hi; i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; if ( m < i2hi ) { i2hi = m; } if ( ihi < i2hi ) { i2hi = ihi; } inc = i2hi + 1 - i2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, " Row:" ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, " %7d ", i - 1 ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); if ( jlo < 1 ) { j2lo = 1; } else { j2lo = jlo; } if ( n < jhi ) { j2hi = n; } else { j2hi = jhi; } for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, "%5d:", j - 1 ); for ( i2 = 1; i2 <= inc; i2++ ) { i = i2lo - 1 + i2; fprintf ( stdout, " %14g", a[(i-1)+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ double *r8mat_uniform_01_new ( int m, int n, int *seed ) /******************************************************************************/ /* Purpose: R8MAT_UNIFORM_01_NEW fills an R8MAT with pseudorandom values scaled to [0,1]. Discussion: An R8MAT is a doubly dimensioned array of R8 values, stored as a vector in column-major order. This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 30 June 2009 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Philip Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, pages 136-143, 1969. Parameters: Input, int M, N, the number of rows and columns. Input/output, int *SEED, the "seed" value. Normally, this value should not be 0, otherwise the output value of SEED will still be 0, and R8_UNIFORM will be 0. On output, SEED has been updated. Output, double R8MAT_UNIFORM_01_NEW[M*N], a matrix of pseudorandom values. */ { int i; const int i4_huge = 2147483647; int j; int k; double *r; r = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r[i+j*m] = ( double ) ( *seed ) * 4.656612875E-10; } } return r; } /******************************************************************************/ double *r8vec_cross_3d ( double v1[3], double v2[3] ) /******************************************************************************/ /* Purpose: R8VEC_CROSS_3D computes the cross product of two vectors in 3D. Licensing: This code is distributed under the MIT license. Modified: 07 August 2005 Author: John Burkardt Parameters: Input, double V1[3], V2[3], the coordinates of the vectors. Output, double R8VEC_CROSS_3D[3], the cross product vector. */ { double *v3; v3 = ( double * ) malloc ( 3 * sizeof ( double ) ); v3[0] = v1[1] * v2[2] - v1[2] * v2[1]; v3[1] = v1[2] * v2[0] - v1[0] * v2[2]; v3[2] = v1[0] * v2[1] - v1[1] * v2[0]; return v3; } /******************************************************************************/ int r8vec_is_nonnegative ( int n, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_IS_NONNEGATIVE is true if all entries in an R8VEC are nonnegative. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 04 August 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double X[N], the vector to be checked. Output, int R8VEC_IS_NONNEGATIVE is true if all elements of X are nonnegative. */ { int i; for ( i = 0; i < n; i++ ) { if ( x[i] < 0.0 ) { return 0; } } return 1; } /******************************************************************************/ int r8vec_is_zero ( int n, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_IS_ZERO is true if the entries in an R8VEC are all zero. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 30 July 2009 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double X[N], the vector to be checked. Output, int R8VEC_IS_ZERO is true if all N elements of X are zero. */ { int i; for ( i = 0; i < n; i++ ) { if ( x[i] != 0.0 ) { return 0; } } return 1; } /******************************************************************************/ double r8vec_length ( int dim_num, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_LENGTH returns the Euclidean length of a R8VEC Licensing: This code is distributed under the MIT license. Modified: 08 August 2005 Author: John Burkardt Parameters: Input, int DIM_NUM, the spatial dimension. Input, double X[DIM_NUM], the vector. Output, double R8VEC_LENGTH, the Euclidean length of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < dim_num; i++ ) { value = value + pow ( x[i], 2 ); } value = sqrt ( value ); return value; } /******************************************************************************/ double r8vec_max ( int n, double dvec[] ) /******************************************************************************/ /* Purpose: R8VEC_MAX returns the maximum element in an R8VEC. Licensing: This code is distributed under the MIT license. Modified: 02 July 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double DVEC[N], a pointer to the first entry of the array. Output, double R8VEC_MAX, the value of the maximum element. This is set to 0.0 if N <= 0. */ { int i; double value; value = - HUGE_VAL; if ( n <= 0 ) { return value; } for ( i = 0; i < n; i++ ) { if ( value < dvec[i] ) { value = dvec[i]; } } return value; } /******************************************************************************/ double r8vec_mean ( int n, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_MEAN returns the mean of a R8VEC. Licensing: This code is distributed under the MIT license. Modified: 02 December 2004 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double X[N], the vector whose mean is desired. Output, double R8VEC_MEAN, the mean, or average, of the vector entries. */ { int i; double mean; mean = 0.0; for ( i = 0; i < n; i++ ) { mean = mean + x[i]; } mean = mean / ( double ) n; return mean; } /******************************************************************************/ double r8vec_min ( int n, double dvec[] ) /******************************************************************************/ /* Purpose: R8VEC_MIN returns the minimum element in an R8VEC. Licensing: This code is distributed under the MIT license. Modified: 02 July 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the array. Input, double DVEC[N], the array to be checked. Output, double R8VEC_MIN, the value of the minimum element. */ { int i; double value; value = HUGE_VAL; if ( n <= 0 ) { return value; } for ( i = 0; i < n; i++ ) { if ( dvec[i] < value ) { value = dvec[i]; } } return value; } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: R8VEC_PRINT prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Parameters: Input, int N, the number of components of the vector. Input, double A[N], the vector to be printed. Input, char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14g\n", i, a[i] ); } return; } /******************************************************************************/ double r8vec_sum ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_SUM returns the sum of an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 15 October 2004 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double A[N], the vector. Output, double R8VEC_SUM, the sum of the vector. */ { int i; double value; value = 0.0; for ( i = 0; i < n; i++ ) { value = value + a[i]; } return value; } /******************************************************************************/ double *r8vec_uniform_01_new ( int n, int *seed ) /******************************************************************************/ /* Purpose: R8VEC_UNIFORM_01_NEW returns a unit pseudorandom R8VEC. Discussion: An R8VEC is a vector of R8's. This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 19 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Parameters: Input, int N, the number of entries in the vector. Input/output, int *SEED, a seed for the random number generator. Output, double R8VEC_UNIFORM_01_NEW[N], the vector of pseudorandom values. */ { int i; int k; double *r; r = ( double * ) malloc ( n * sizeof ( double ) ); for ( i = 0; i < n; i++ ) { k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } r[i] = ( double ) ( *seed ) * 4.656612875E-10; } return r; } /******************************************************************************/ double r8vec_variance ( int n, double x[] ) /******************************************************************************/ /* Purpose: R8VEC_VARIANCE returns the variance of a double vector. Licensing: This code is distributed under the MIT license. Modified: 01 May 1999 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, double X[N], the vector whose variance is desired. Output, double R8VEC_VARIANCE, the variance of the vector entries. */ { int i; double mean; double variance; mean = r8vec_mean ( n, x ); variance = 0.0; for ( i = 0; i < n; i++ ) { variance = variance + ( x[i] - mean ) * ( x[i] - mean ); } if ( 1 < n ) { variance = variance / ( double ) ( n - 1 ); } else { variance = 0.0; } return variance; } /******************************************************************************/ void r8vec_zero ( int n, double a[] ) /******************************************************************************/ /* Purpose: R8VEC_ZERO zeroes a real vector. Licensing: This code is distributed under the MIT license. Modified: 03 July 2005 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Output, double A[N], a vector of zeroes. */ { int i; for ( i = 0; i < n; i++ ) { a[i] = 0.0; } return; } /******************************************************************************/ int s_len_trim ( char *s ) /******************************************************************************/ /* Purpose: S_LEN_TRIM returns the length of a string to the last nonblank. Discussion: It turns out that I also want to ignore the '\n' character! Licensing: This code is distributed under the MIT license. Modified: 05 October 2014 Author: John Burkardt Parameters: Input, char *S, a pointer to a string. Output, int S_LEN_TRIM, the length of the string to the last nonblank. If S_LEN_TRIM is 0, then the string is entirely blank. */ { int n; char *t; n = strlen ( s ); t = s + strlen ( s ) - 1; while ( 0 < n ) { if ( *t != ' ' && *t != '\n' ) { return n; } t--; n--; } return n; } /******************************************************************************/ void sort_heap_external ( int n, int *indx, int *i, int *j, int isgn ) /******************************************************************************/ /* Purpose: SORT_HEAP_EXTERNAL externally sorts a list of items into ascending order. Discussion: The actual list is not passed to the routine. Hence it may consist of integers, reals, numbers, names, etc. The user, after each return from the routine, will be asked to compare or interchange two items. The current version of this code mimics the FORTRAN version, so the values of I and J, in particular, are FORTRAN indices. Licensing: This code is distributed under the MIT license. Modified: 05 February 2004 Author: Original FORTRAN77 version by Albert Nijenhuis, Herbert Wilf. C++ version by John Burkardt. Reference: Albert Nijenhuis, Herbert Wilf, Combinatorial Algorithms, Academic Press, 1978, second edition, ISBN 0-12-519260-6. Parameters: Input, int N, the length of the input list. Input/output, int *INDX. The user must set INDX to 0 before the first call. On return, if INDX is greater than 0, the user must interchange items I and J and recall the routine. If INDX is less than 0, the user is to compare items I and J and return in ISGN a negative value if I is to precede J, and a positive value otherwise. If INDX is 0, the sorting is done. Output, int *I, *J. On return with INDX positive, elements I and J of the user's list should be interchanged. On return with INDX negative, elements I and J are to be compared by the user. Input, int ISGN. On return with INDX negative, the user should compare elements I and J of the list. If item I is to precede item J, set ISGN negative, otherwise set ISGN positive. */ { static int i_save = 0; static int j_save = 0; static int k = 0; static int k1 = 0; static int n1 = 0; /* INDX = 0: This is the first call. */ if ( *indx == 0 ) { i_save = 0; j_save = 0; k = n / 2; k1 = k; n1 = n; } /* INDX < 0: The user is returning the results of a comparison. */ else if ( *indx < 0 ) { if ( *indx == -2 ) { if ( isgn < 0 ) { i_save = i_save + 1; } j_save = k1; k1 = i_save; *indx = -1; *i = i_save; *j = j_save; return; } if ( 0 < isgn ) { *indx = 2; *i = i_save; *j = j_save; return; } if ( k <= 1 ) { if ( n1 == 1 ) { i_save = 0; j_save = 0; *indx = 0; } else { i_save = n1; j_save = 1; n1 = n1 - 1; *indx = 1; } *i = i_save; *j = j_save; return; } k = k - 1; k1 = k; } /* 0 < INDX: the user was asked to make an interchange. */ else if ( *indx == 1 ) { k1 = k; } for ( ; ; ) { i_save = 2 * k1; if ( i_save == n1 ) { j_save = k1; k1 = i_save; *indx = -1; *i = i_save; *j = j_save; return; } else if ( i_save <= n1 ) { j_save = i_save + 1; *indx = -2; *i = i_save; *j = j_save; return; } if ( k <= 1 ) { break; } k = k - 1; k1 = k; } if ( n1 == 1 ) { i_save = 0; j_save = 0; *indx = 0; *i = i_save; *j = j_save; } else { i_save = n1; j_save = 1; n1 = n1 - 1; *indx = 1; *i = i_save; *j = j_save; } return; } /******************************************************************************/ int *tet_mesh_neighbor_tets ( int tetra_order, int tetra_num, int tetra_node[] ) /******************************************************************************/ /* Purpose: TET_MESH_NEIGHBOR_TETS determines tetrahedron neighbors. Discussion: A tet mesh of a set of nodes can be completely described by the coordinates of the nodes, and the list of nodes that make up each tetrahedron. In the most common case, four nodes are used. There is also a 10 node case, where nodes are also placed on the midsides of the tetrahedral edges. This routine can handle 4 or 10-node tetrahedral meshes. The 10-node case is handled simply by ignoring the six midside nodes, which are presumed to be listed after the vertices. The tetrahedron adjacency information records which tetrahedron is adjacent to a given tetrahedron on a particular face. This routine creates a data structure recording this information. The primary amount of work occurs in sorting a list of 4 * TETRA_NUM data items. The neighbor tetrahedrons are indexed by the face they share with the tetrahedron. Each face of the tetrahedron is indexed by the node which is NOT part of the face. That is: * Neighbor 1 shares face 1 defined by nodes 2, 3, 4. * Neighbor 2 shares face 2 defined by nodes 1, 3, 4; * Neighbor 3 shares face 3 defined by nodes 1, 2, 4; * Neighbor 4 shares face 4 defined by nodes 1, 2, 3. For instance, if the (transposed) TETRA_NODE array was: Row 1 2 3 4 Col 1 4 3 5 1 2 4 2 5 1 3 4 7 3 5 4 4 7 8 5 5 4 6 2 5 6 4 6 8 5 then the (transposed) TETRA_NEIGHBOR array should be: Row 1 2 3 4 Col 1 -1 2 -1 3 2 -1 1 -1 5 3 -1 1 4 -1 4 -1 6 3 -1 5 -1 2 6 -1 6 -1 4 5 -1 Licensing: This code is distributed under the MIT license. Modified: 12 December 2006 Author: John Burkardt Parameters: Input, int TETRA_ORDER, the order of the tetrahedrons. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[TETRA_ORDER*TETRA_NUM], the indices of the nodes. Output, int TET_MESH_NEIGHBORS[4*TETRA_NUM], the four tetrahedrons that are direct neighbors of a given tetrahedron. If there is no neighbor sharing a given face, the index is set to -1. */ { int a; int b; int c; int face; int face1; int face2; int *faces; int i; int j; int k; int l; int tetra; int *tetra_neighbor; int tetra1; int tetra2; faces = ( int * ) malloc ( 5 * 4 * tetra_num * sizeof ( int ) ); tetra_neighbor = ( int * ) malloc ( 4 * tetra_num * sizeof ( int ) ); /* Step 1. From the list of nodes for tetrahedron T, of the form: (I,J,K,L) construct the four face relations: (J,K,L,1,T) (I,K,L,2,T) (I,J,L,3,T) (I,J,K,4,T) In order to make matching easier, we reorder each triple of nodes into ascending order. */ for ( tetra = 0; tetra < tetra_num; tetra++ ) { i = tetra_node[0+tetra*tetra_order]; j = tetra_node[1+tetra*tetra_order]; k = tetra_node[2+tetra*tetra_order]; l = tetra_node[3+tetra*tetra_order]; i4i4i4_sort_a ( j, k, l, &a, &b, &c ); faces[0+0*5+tetra*5*4] = a; faces[1+0*5+tetra*5*4] = b; faces[2+0*5+tetra*5*4] = c; faces[3+0*5+tetra*5*4] = 0; faces[4+0*5+tetra*5*4] = tetra; i4i4i4_sort_a ( i, k, l, &a, &b, &c ); faces[0+1*5+tetra*5*4] = a; faces[1+1*5+tetra*5*4] = b; faces[2+1*5+tetra*5*4] = c; faces[3+1*5+tetra*5*4] = 1; faces[4+1*5+tetra*5*4] = tetra; i4i4i4_sort_a ( i, j, l, &a, &b, &c ); faces[0+2*5+tetra*5*4] = a; faces[1+2*5+tetra*5*4] = b; faces[2+2*5+tetra*5*4] = c; faces[3+2*5+tetra*5*4] = 2; faces[4+2*5+tetra*5*4] = tetra; i4i4i4_sort_a ( i, j, k, &a, &b, &c ); faces[0+3*5+tetra*5*4] = a; faces[1+3*5+tetra*5*4] = b; faces[2+3*5+tetra*5*4] = c; faces[3+3*5+tetra*5*4] = 3; faces[4+3*5+tetra*5*4] = tetra; } /* Step 2. Perform an ascending dictionary sort on the neighbor relations. We only intend to sort on rows 1:3; the routine we call here sorts on rows 1 through 5 but that won't hurt us. What we need is to find cases where two tetrahedrons share a face. By sorting the columns of the FACES array, we will put shared faces next to each other. */ i4col_sort_a ( 5, 4*tetra_num, faces ); /* Step 3. Neighboring tetrahedrons show up as consecutive columns with identical first three entries. Whenever you spot this happening, make the appropriate entries in TETRA_NEIGHBOR. */ for ( j = 0; j < tetra_num; j++ ) { for ( i = 0; i < 4; i++ ) { tetra_neighbor[i+j*4] = -1; } } face = 0; for ( ; ; ) { if ( 4 * tetra_num - 1 <= face ) { break; } if ( faces[0+face*5] == faces[0+(face+1)*5] && faces[1+face*5] == faces[1+(face+1)*5] && faces[2+face*5] == faces[2+(face+1)*5] ) { face1 = faces[3+face*5]; tetra1 = faces[4+face*5]; face2 = faces[3+(face+1)*5]; tetra2 = faces[4+(face+1)*5]; tetra_neighbor[face1+tetra1*4] = tetra2; tetra_neighbor[face2+tetra2*4] = tetra1; face = face + 2; } else { face = face + 1; } } free ( faces ); return tetra_neighbor; } /******************************************************************************/ int *tet_mesh_node_order ( int tetra_order, int tetra_num, int tetra_node[], int node_num ) /******************************************************************************/ /* Purpose: TET_MESH_NODE_ORDER: determines the order of nodes. Discussion: The order of a node is the number of tetrahedrons that use that node as a vertex. Licensing: This code is distributed under the MIT license. Modified: 27 October 2005 Author: John Burkardt Parameters: Input, int TETRA_ORDER, the order of the tetrahedrons. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[TETRA_ORDER*TETRA_NUM], the indices of the nodes. Input, int NODE_NUM, the number of nodes. Output, int TET_MESH_NODE_ORDER[NODE_NUM], the order of each node. */ { int i; int node; int *node_order; int tetra; node_order = ( int * ) malloc ( node_num * sizeof ( int ) ); i4vec_zero ( node_num, node_order ); for ( tetra = 0; tetra < tetra_num; tetra++ ) { for ( i = 0; i < tetra_order; i++ ) { node = tetra_node[i+tetra*tetra_order]; if ( node < 0 || node_num <= node ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TET_MESH_NODE_ORDER - Fatal error!\n" ); fprintf ( stderr, " Illegal entry in TETRA_NODE.\n" ); exit ( 1 ); } else { node_order[node] = node_order[node] + 1; } } } return node_order; } /******************************************************************************/ void tet_mesh_order4_adj_count ( int node_num, int tetra_num, int tetra_node[], int *adj_num, int adj_row[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_ADJ_COUNT counts the number of nodal adjacencies. Discussion: Assuming that the tet mesh is to be used in a finite element computation, we declare that two distinct nodes are "adjacent" if and only if they are both included in some tetrahedron. It is the purpose of this routine to determine the number of such adjacency relationships. The initial count gets only the (I,J) relationships, for which node I is strictly less than node J. This value is doubled to account for symmetry. Licensing: This code is distributed under the MIT license. Modified: 13 January 2007 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[4*TETRA_NUM], the indices of the nodes. Output, int *ADJ_NUM, the total number of adjacency relationships, Output, int ADJ_ROW[NODE_NUM+1], the ADJ pointer array. */ { int i; int j; int k; int node; int *pair; int pair_num; int pair_unique_num; int tetra; /* Each order 4 tetrahedron defines 6 adjacency pairs. */ pair = ( int * ) malloc ( 2 * 6 * tetra_num * sizeof ( int ) ); for ( tetra = 0; tetra < tetra_num; tetra++ ) { pair[0+ tetra *2] = tetra_node[0+tetra*4]; pair[1+ tetra *2] = tetra_node[1+tetra*4]; pair[0+( tetra_num+tetra)*2] = tetra_node[0+tetra*4]; pair[1+( tetra_num+tetra)*2] = tetra_node[2+tetra*4]; pair[0+(2*tetra_num+tetra)*2] = tetra_node[0+tetra*4]; pair[1+(2*tetra_num+tetra)*2] = tetra_node[3+tetra*4]; pair[0+(3*tetra_num+tetra)*2] = tetra_node[1+tetra*4]; pair[1+(3*tetra_num+tetra)*2] = tetra_node[2+tetra*4]; pair[0+(4*tetra_num+tetra)*2] = tetra_node[1+tetra*4]; pair[1+(4*tetra_num+tetra)*2] = tetra_node[3+tetra*4]; pair[0+(5*tetra_num+tetra)*2] = tetra_node[2+tetra*4]; pair[1+(5*tetra_num+tetra)*2] = tetra_node[3+tetra*4]; } pair_num = 6 * tetra_num; /* Force the nodes of each pair to be listed in ascending order. */ i4mat_transpose_print_some ( 2, pair_num, pair, 1, 1, 2, pair_num, "DEBUG: PAIR before first sort" ); i4col_sort2_a ( 2, pair_num, pair ); i4mat_transpose_print_some ( 2, pair_num, pair, 1, 1, 2, pair_num, "DEBUG: PAIR after first sort" ); /* Rearrange the columns in ascending order. */ i4col_sort_a ( 2, pair_num, pair ); /* Get the number of unique columns. */ pair_unique_num = i4col_sorted_unique_count ( 2, pair_num, pair ); /* The number of adjacencies is TWICE this value, plus the number of nodes. */ *adj_num = 2 * pair_unique_num; /* Now set up the ADJ_ROW counts. */ for ( node = 0; node < node_num; node++ ) { adj_row[node] = 0; } for ( k = 0; k < pair_num; k++ ) { if ( 0 < k ) { if ( pair[0+(k-1)*2] == pair[0+k*2] && pair[1+(k-1)*2] == pair[1+k*2] ) { continue; } } i = pair[0+k*2]; j = pair[1+k*2]; adj_row[i-1] = adj_row[i-1] + 1; adj_row[j-1] = adj_row[j-1] + 1; } /* We used ADJ_ROW to count the number of entries in each row. Convert it to pointers into the ADJ array. */ for ( node = node_num-1; 0 <= node; node-- ) { adj_row[node] = adj_row[node+1]; } adj_row[0] = 1; for ( node = 1; node <= node_num; node++ ) { adj_row[node] = adj_row[node-1] + adj_row[i]; } free ( pair ); return; } /******************************************************************************/ int *tet_mesh_order4_adj_set ( int node_num, int element_num, int element_node[], int adj_num, int adj_row[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_ADJ_SET sets the nodal adjacency matrix. Discussion: A compressed format is used for the nodal adjacency matrix. It is assumed that we know ADJ_NUM, the number of adjacency entries and the ADJ_ROW array, which keeps track of the list of slots in ADJ where we can store adjacency information for each row. We essentially repeat the work of TET_MESH_ORDER4_ADJ_COUNT, but now we have a place to store the adjacency information. A copy of the ADJ_ROW array is useful, as we can use it to keep track of the next available entry in ADJ for adjacencies associated with a given row. Licensing: This code is distributed under the MIT license. Modified: 13 January 2007 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[4*TETRA_NUM], the indices of the nodes. Input, int ADJ_NUM, the total number of adjacency relationships, Input, int ADJ_ROW[NODE_NUM+1], the ADJ pointer array. Output, int TET_MESH_ORDER4_ADJ_SET[ADJ_NUM], the adjacency information. */ { int *adj; int *adj_row_copy; int i; int j; int k; int node; int *pair; int pair_num; int tetra; /* Each order 4 tetrahedron defines 6 adjacency pairs. */ pair = ( int * ) malloc ( 2 * 6 * element_num * sizeof ( int ) ); for ( tetra = 0; tetra < element_num; tetra++ ) { pair[0+ tetra *2] = element_node[0+tetra*4]; pair[1+ tetra *2] = element_node[1+tetra*4]; pair[0+( element_num+tetra)*2] = element_node[0+tetra*4]; pair[1+( element_num+tetra)*2] = element_node[2+tetra*4]; pair[0+(2*element_num+tetra)*2] = element_node[0+tetra*4]; pair[1+(2*element_num+tetra)*2] = element_node[3+tetra*4]; pair[0+(3*element_num+tetra)*2] = element_node[1+tetra*4]; pair[1+(3*element_num+tetra)*2] = element_node[2+tetra*4]; pair[0+(4*element_num+tetra)*2] = element_node[1+tetra*4]; pair[1+(4*element_num+tetra)*2] = element_node[3+tetra*4]; pair[0+(5*element_num+tetra)*2] = element_node[2+tetra*4]; pair[1+(5*element_num+tetra)*2] = element_node[3+tetra*4]; } pair_num = 6 * element_num; /* Force the nodes of each pair to be listed in ascending order. */ i4col_sort2_a ( 2, pair_num, pair ); /* Rearrange the columns in ascending order. */ i4col_sort_a ( 2, pair_num, pair ); /* Mark all entries of ADJ so we will know later if we missed one. */ adj = ( int * ) malloc ( adj_num * sizeof ( int ) ); for ( i = 0; i < adj_num; i++ ) { adj[i] = -1; } /* Copy the ADJ_ROW array and use it to keep track of the next free entry for each row. */ adj_row_copy = ( int * ) malloc ( node_num * sizeof ( int ) ); for ( node = 0; node < node_num; node++ ) { adj_row_copy[node] = adj_row[node]; } /* Now set up the ADJ_ROW counts. */ for ( k = 0; k < pair_num; k++ ) { if ( 0 < k ) { if ( pair[0+(k-1)*2] == pair[0+k*2] && pair[1+(k-1)*2] == pair[1+k*2] ) { continue; } } i = pair[0+k*2]; j = pair[1+k*2]; adj[adj_row_copy[i]] = j; adj_row_copy[i] = adj_row_copy[i] + 1; adj[adj_row_copy[j]] = i; adj_row_copy[j] = adj_row_copy[j] + 1; } free ( adj_row_copy ); free ( pair ); return adj; } /******************************************************************************/ int tet_mesh_order4_boundary_face_count ( int tetra_num, int tetra_node[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_BOUNDARY_FACE_COUNT counts the number of boundary faces. Discussion: This routine is given a tet mesh, an abstract list of quadruples of nodes. It is assumed that the nodes forming each face of each tetrahedron are listed in a counterclockwise order, although the routine should work if the nodes are consistently listed in a clockwise order as well. It is assumed that each face of the tet mesh is either * an INTERIOR face, which is listed twice, once with positive orientation and once with negative orientation, or; * a BOUNDARY face, which will occur only once. This routine should work even if the region has holes. Licensing: This code is distributed under the MIT license. Modified: 11 January 2007 Author: John Burkardt Parameters: Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[4*TETRA_NUM], the indices of the nodes. Output, int TET_MESH_ORDER4_BOUNDARY_FACE_COUNT, the number of boundary faces. */ { int boundary_face_num; int *face; int face_num; int interior_face_num; int m; int tet; int unique_face_num; face = ( int * ) malloc ( 3 * 4 * tetra_num * sizeof ( int ) ); m = 3; face_num = 4 * tetra_num; /* Set up the face array: (Omit node 1) (Omit node 2) (Omit node 3) (Omit node 4) */ for ( tet = 0; tet < tetra_num; tet++ ) { face[0+( tet)*3] = tetra_node[1+tet*4]; face[1+( tet)*3] = tetra_node[2+tet*4]; face[2+( tet)*3] = tetra_node[3+tet*4]; face[0+( tetra_num+tet)*3] = tetra_node[0+tet*4]; face[1+( tetra_num+tet)*3] = tetra_node[2+tet*4]; face[2+( tetra_num+tet)*3] = tetra_node[3+tet*4]; face[0+(2*tetra_num+tet)*3] = tetra_node[0+tet*4]; face[1+(2*tetra_num+tet)*3] = tetra_node[1+tet*4]; face[2+(2*tetra_num+tet)*3] = tetra_node[3+tet*4]; face[0+(3*tetra_num+tet)*3] = tetra_node[0+tet*4]; face[1+(3*tetra_num+tet)*3] = tetra_node[1+tet*4]; face[2+(3*tetra_num+tet)*3] = tetra_node[2+tet*4]; } /* Force the nodes of each face to be listed in ascending order. */ i4col_sort2_a ( m, face_num, face ); /* Ascending sort the columns. */ i4col_sort_a ( m, face_num, face ); /* Get the number of unique columns. */ unique_face_num = i4col_sorted_unique_count ( m, face_num, face ); /* Determine the number of interior and boundary faces. */ interior_face_num = 4 * tetra_num - unique_face_num; boundary_face_num = 4 * tetra_num - 2 * interior_face_num; free ( face ); return boundary_face_num; } /******************************************************************************/ int tet_mesh_order4_edge_count ( int tetra_num, int tetra_node[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_EDGE_COUNT counts the number of edges. Discussion: This routine is given a tet mesh, an abstract list of quadruples of nodes. Each tetrahedron defines 6 edges; however, assuming that tetrahedrons are touching each other, most edges will be used more than once. This routine determines the actual number of "geometric" edges associated with the tet mesh. Licensing: This code is distributed under the MIT license. Modified: 11 January 2007 Author: John Burkardt Parameters: Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[4*TETRA_NUM], the indices of the nodes. Output, int TET_MESH_ORDER4_EDGE_COUNT, the number of edges. */ { int *edge; int edge_num; int edge_num_raw; int m; int tet; edge = ( int * ) malloc ( 2 * 6 * tetra_num * sizeof ( int ) ); m = 3; edge_num_raw = 6 * tetra_num; /* Set up the raw edge array: */ for ( tet = 0; tet < tetra_num; tet++ ) { edge[0+ tet *2] = tetra_node[0+tet*4]; edge[1+ tet *2] = tetra_node[1+tet*4]; edge[0+( tetra_num+tet)*2] = tetra_node[0+tet*4]; edge[1+( tetra_num+tet)*2] = tetra_node[2+tet*4]; edge[0+(2*tetra_num+tet)*2] = tetra_node[0+tet*4]; edge[1+(2*tetra_num+tet)*2] = tetra_node[3+tet*4]; edge[0+(3*tetra_num+tet)*2] = tetra_node[1+tet*4]; edge[1+(3*tetra_num+tet)*2] = tetra_node[2+tet*4]; edge[0+(4*tetra_num+tet)*2] = tetra_node[1+tet*4]; edge[1+(4*tetra_num+tet)*2] = tetra_node[3+tet*4]; edge[0+(5*tetra_num+tet)*2] = tetra_node[2+tet*4]; edge[1+(5*tetra_num+tet)*2] = tetra_node[3+tet*4]; } /* Force the nodes of each face to be listed in ascending order. */ i4col_sort2_a ( m, edge_num_raw, edge ); /* Ascending sort the columns. */ i4col_sort_a ( m, edge_num_raw, edge ); /* Get the number of unique columns. */ edge_num = i4col_sorted_unique_count ( m, edge_num_raw, edge ); free ( edge ); return edge_num; } /******************************************************************************/ void tet_mesh_order4_example_set ( int node_num, int tetra_num, double node_xyz[], int tetra_node[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_EXAMPLE_SET sets an example linear tet mesh. Licensing: This code is distributed under the MIT license. Modified: 19 August 2009 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int TETRA_NUM, the number of tetrahedrons. Output, double NODE_XYZ[3*NODE_NUM], the node coordinates. Output, int TETRA_NODE[4*TETRA_NUM], the nodes forming each tet. */ { int i; int j; double node_xyz_save[3*63] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.5, 0.0, 0.0, 1.0, 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.0, 0.5, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 0.5, 0.0, 1.0, 1.0, 0.5, 0.0, 0.0, 0.5, 0.0, 0.5, 0.5, 0.0, 1.0, 0.5, 0.5, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 0.5, 1.0, 0.0, 0.5, 1.0, 0.5, 0.5, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.5, 1.0, 0.0, 1.0, 1.0, 0.5, 0.0, 1.0, 0.5, 0.5, 1.0, 0.5, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.5, 1.0, 1.0, 1.0, 1.5, 0.0, 0.0, 1.5, 0.0, 0.5, 1.5, 0.0, 1.0, 1.5, 0.5, 0.0, 1.5, 0.5, 0.5, 1.5, 0.5, 1.0, 1.5, 1.0, 0.0, 1.5, 1.0, 0.5, 1.5, 1.0, 1.0, 2.0, 0.0, 0.0, 2.0, 0.0, 0.5, 2.0, 0.0, 1.0, 2.0, 0.5, 0.0, 2.0, 0.5, 0.5, 2.0, 0.5, 1.0, 2.0, 1.0, 0.0, 2.0, 1.0, 0.5, 2.0, 1.0, 1.0, 2.5, 0.0, 0.0, 2.5, 0.0, 0.5, 2.5, 0.0, 1.0, 2.5, 0.5, 0.0, 2.5, 0.5, 0.5, 2.5, 0.5, 1.0, 2.5, 1.0, 0.0, 2.5, 1.0, 0.5, 2.5, 1.0, 1.0, 3.0, 0.0, 0.0, 3.0, 0.0, 0.5, 3.0, 0.0, 1.0, 3.0, 0.5, 0.0, 3.0, 0.5, 0.5, 3.0, 0.5, 1.0, 3.0, 1.0, 0.0, 3.0, 1.0, 0.5, 3.0, 1.0, 1.0 }; int tetra_node_save[4*144] = { 1, 2, 4, 10, 2, 4, 5, 10, 2, 5, 10, 11, 2, 3, 5, 11, 4, 5, 10, 13, 3, 5, 6, 11, 5, 10, 11, 13, 4, 5, 7, 13, 5, 6, 8, 14, 5, 7, 8, 13, 6, 8, 9, 14, 11, 13, 14, 19, 12, 14, 15, 20, 3, 6, 11, 12, 5, 6, 11, 14, 6, 9, 14, 15, 6, 11, 12, 14, 6, 12, 14, 15, 7, 8, 13, 16, 5, 8, 13, 14, 10, 11, 13, 19, 8, 9, 14, 17, 11, 12, 14, 20, 5, 11, 13, 14, 8, 13, 14, 16, 9, 14, 15, 17, 13, 14, 16, 22, 8, 14, 16, 17, 14, 15, 17, 23, 14, 16, 17, 22, 9, 15, 17, 18, 15, 17, 18, 23, 14, 17, 22, 23, 13, 14, 19, 22, 11, 14, 19, 20, 14, 15, 20, 23, 15, 20, 21, 23, 21, 23, 24, 29, 20, 22, 23, 28, 14, 19, 20, 22, 15, 18, 23, 24, 12, 15, 20, 21, 15, 21, 23, 24, 16, 17, 22, 25, 19, 20, 22, 28, 17, 18, 23, 26, 20, 21, 23, 29, 14, 20, 22, 23, 17, 22, 23, 25, 18, 23, 24, 26, 22, 23, 25, 31, 17, 23, 25, 26, 23, 24, 26, 32, 23, 25, 26, 31, 18, 24, 26, 27, 24, 26, 27, 32, 23, 26, 31, 32, 22, 23, 28, 31, 20, 23, 28, 29, 23, 24, 29, 32, 24, 29, 30, 32, 30, 32, 33, 38, 29, 31, 32, 37, 23, 28, 29, 31, 24, 27, 32, 33, 21, 24, 29, 30, 24, 30, 32, 33, 25, 26, 31, 34, 28, 29, 31, 37, 26, 27, 32, 35, 29, 30, 32, 38, 23, 29, 31, 32, 26, 31, 32, 34, 27, 32, 33, 35, 31, 32, 34, 40, 26, 32, 34, 35, 32, 33, 35, 41, 32, 34, 35, 40, 27, 33, 35, 36, 33, 35, 36, 41, 32, 35, 40, 41, 31, 32, 37, 40, 29, 32, 37, 38, 32, 33, 38, 41, 33, 38, 39, 41, 39, 41, 42, 47, 38, 40, 41, 46, 32, 37, 38, 40, 33, 36, 41, 42, 30, 33, 38, 39, 33, 39, 41, 42, 34, 35, 40, 43, 37, 38, 40, 46, 35, 36, 41, 44, 38, 39, 41, 47, 32, 38, 40, 41, 35, 40, 41, 43, 36, 41, 42, 44, 40, 41, 43, 49, 35, 41, 43, 44, 41, 42, 44, 50, 41, 43, 44, 49, 36, 42, 44, 45, 42, 44, 45, 50, 41, 44, 49, 50, 40, 41, 46, 49, 38, 41, 46, 47, 41, 42, 47, 50, 42, 47, 48, 50, 48, 50, 51, 56, 47, 49, 50, 55, 41, 46, 47, 49, 42, 45, 50, 51, 39, 42, 47, 48, 42, 48, 50, 51, 43, 44, 49, 52, 46, 47, 49, 55, 44, 45, 50, 53, 47, 48, 50, 56, 41, 47, 49, 50, 44, 49, 50, 52, 45, 50, 51, 53, 49, 50, 52, 58, 44, 50, 52, 53, 50, 51, 53, 59, 50, 52, 53, 58, 45, 51, 53, 54, 51, 53, 54, 59, 50, 53, 58, 59, 49, 50, 55, 58, 47, 50, 55, 56, 50, 51, 56, 59, 51, 56, 57, 59, 50, 55, 56, 58, 51, 54, 59, 60, 48, 51, 56, 57, 51, 57, 59, 60, 52, 53, 58, 61, 53, 54, 59, 62, 50, 56, 58, 59, 53, 58, 59, 61, 54, 59, 60, 62, 53, 59, 61, 62, 54, 60, 62, 63 }; for ( j = 0; j < node_num; j++ ) { for ( i = 0; i < 3; i++ ) { node_xyz[i+j*3] = node_xyz_save[i+j*3]; } } for ( j = 0; j < tetra_num; j++ ) { for ( i = 0; i < 4; i++ ) { tetra_node[i+j*4] = tetra_node_save[i+j*4] - 1; } } return; } /******************************************************************************/ void tet_mesh_order4_example_size ( int *node_num, int *tetra_num ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_EXAMPLE_SIZE sizes an example linear tet mesh. Licensing: This code is distributed under the MIT license. Modified: 03 August 2009 Author: John Burkardt Parameters: Output, int *NODE_NUM, the number of nodes. Output, int *TETRA_NUM, the number of tetrahedrons. */ { *node_num = 63; *tetra_num = 144; return; } /******************************************************************************/ void tet_mesh_order4_refine_compute ( int node_num1, int tetra_num1, double node_xyz1[], int tetra_node1[], int node_num2, int tetra_num2, int edge_data[], double node_xyz2[], int tetra_node2[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_REFINE_COMPUTE computes a refined order 4 tet mesh Discussion: A refined 4-node tet mesh can be derived from a given 4-node tet mesh by interpolating nodes at the midpoint of every edge of the mesh. The mesh is described indirectly, as the sum of individual tetrahedrons. A single physical edge may be a logical edge of any number of tetrahedrons. It is important, however, that a new node be created exactly once for each edge, assigned an index, and associated with every tetrahedron that shares this edge. This routine handles that problem. The primary amount of work occurs in sorting a list of 6 * TETRA_NUM data items, one item for every edge of every tetrahedron. Each data item records, for a given tetrahedron edge, the global indices of the two endpoints, the local indices of the two endpoints, and the index of the tetrahedron. Through careful sorting, it is possible to arrange this data in a way that allows the proper generation of the interpolated nodes. Let us add the new nodes and temporarily assign them local indices 5 through X, based on the following ordering: 1, 2, 3, 4, (1+2), (1+3), (1+4), (2+3), (2+4), (3+4). Then let us assign these nodes to eight subtetrahedrons as follows: 1, 5, 6, 7 2, 5, 8, 9 3, 6, 8, 9 4, 7, 9, X 5, 6, 7, 9 5, 6, 8, 9 6, 7, 9, X 6, 8, 9, X Licensing: This code is distributed under the MIT license. Modified: 25 January 2007 Author: John Burkardt Reference: Anwei Liu, Barry Joe, Quality Local Refinement of Tetrahedral Meshes Based on 8-Subtetrahedron Subdivision, Mathematics of Computation, Volume 65, Number 215, July 1996, pages 1183-1200. Parameters: Input, int NODE_NUM1, the number of nodes in the input mesh. Input, int TETRA_NUM1, the number of tetrahedrons in the input mesh. Input, double NODE_XYZ1[3*NODE_NUM1], the coordinates of the nodes that make up the input mesh. Input, int TETRA_NODE1[4*TETRA_NUM], the indices of the nodes in the input mesh. Input, int NODE_NUM2, the number of nodes in the refined mesh. Input, int TETRA_NUM2, the number of tetrahedrons in the refined mesh. Input, int EDGE_DATA[5*(6*TETRA_NUM1)], edge data. Output, double NODE_XYZ2[3*NODE_NUM2], the coordinates of the nodes that make up the refined mesh. Output, int TETRA_NODE2[4*TETRA_NUM2], the indices of the nodes in the refined mesh. */ { int dim_num = 3; int edge; int i; int j; int n1; int n1_old; int n2; int n2_old; int node; int tetra_order = 4; int tetra1; int v1; int v2; /* Generate the index and coordinates of the new midside nodes, and update the tetradehron-node data. */ for ( j = 0; j < node_num1; j++ ) { for ( i = 0; i < dim_num; i++ ) { node_xyz2[i+j*dim_num] = node_xyz1[i+j*dim_num]; } } for ( j = 0; j < tetra_num2; j++ ) { for ( i = 0; i < tetra_order; i++ ) { tetra_node2[i+j*tetra_order] = -1; } } /* The vertices of the input tetrahedron can be assigned now. */ for ( tetra1 = 0; tetra1 < tetra_num1; tetra1++ ) { tetra_node2[0+(tetra1*8+0)*tetra_order] = tetra_node1[0+tetra1*tetra_order]; tetra_node2[0+(tetra1*8+1)*tetra_order] = tetra_node1[1+tetra1*tetra_order]; tetra_node2[0+(tetra1*8+2)*tetra_order] = tetra_node1[2+tetra1*tetra_order]; tetra_node2[0+(tetra1*8+3)*tetra_order] = tetra_node1[3+tetra1*tetra_order]; } node = node_num1; n1_old = -1; n2_old = -1; for ( edge = 0; edge < 6 * tetra_num1; edge++ ) { /* Read the data defining the edge. */ n1 = edge_data[0+edge*5]; n2 = edge_data[1+edge*5]; /* If this edge is new, create the coordinates and index. */ if ( n1 != n1_old || n2 != n2_old ) { if ( node_num2 <= node ) { printf ( "\n" ); printf ( "TET_MESH_ORDER4_REFINE_COMPUTE - Fatal error!\n" ); printf ( " Node index exceeds NODE_NUM2.\n" ); exit ( 1 ); } for ( i = 0; i < dim_num; i++ ) { node_xyz2[i+node*dim_num] = ( node_xyz2[i+(n1-1)*dim_num] + node_xyz2[i+(n2-1)*dim_num] ) / 2.0; } node = node + 1; n1_old = n1; n2_old = n2; } /* Assign the node to the tetrahedron. */ v1 = edge_data[2+edge*5]; v2 = edge_data[3+edge*5]; tetra1 = edge_data[4+edge*5]; /* We know the two vertices that bracket this new node. This tells us whether it is new node number 5, 6, 7, 8, 9 or 10. This tells us which of the new subtetrahedrons it belongs to, and what position it occupies. */ if ( v1 == 1 && v2 == 2 ) { tetra_node2[1+(tetra1*8+0)*tetra_order] = node; tetra_node2[1+(tetra1*8+1)*tetra_order] = node; tetra_node2[0+(tetra1*8+4)*tetra_order] = node; tetra_node2[0+(tetra1*8+5)*tetra_order] = node; } else if ( v1 == 1 && v2 == 3 ) { tetra_node2[2+(tetra1*8+0)*tetra_order] = node; tetra_node2[1+(tetra1*8+2)*tetra_order] = node; tetra_node2[1+(tetra1*8+4)*tetra_order] = node; tetra_node2[1+(tetra1*8+5)*tetra_order] = node; tetra_node2[0+(tetra1*8+6)*tetra_order] = node; tetra_node2[0+(tetra1*8+7)*tetra_order] = node; } else if ( v1 == 1 && v2 == 4 ) { tetra_node2[3+(tetra1*8+0)*tetra_order] = node; tetra_node2[1+(tetra1*8+3)*tetra_order] = node; tetra_node2[2+(tetra1*8+4)*tetra_order] = node; tetra_node2[1+(tetra1*8+6)*tetra_order] = node; } else if ( v1 == 2 && v2 == 3 ) { tetra_node2[2+(tetra1*8+1)*tetra_order] = node; tetra_node2[2+(tetra1*8+2)*tetra_order] = node; tetra_node2[2+(tetra1*8+5)*tetra_order] = node; tetra_node2[1+(tetra1*8+7)*tetra_order] = node; } else if ( v1 == 2 && v2 == 4 ) { tetra_node2[3+(tetra1*8+1)*tetra_order] = node; tetra_node2[3+(tetra1*8+2)*tetra_order] = node; tetra_node2[2+(tetra1*8+3)*tetra_order] = node; tetra_node2[3+(tetra1*8+4)*tetra_order] = node; tetra_node2[3+(tetra1*8+5)*tetra_order] = node; tetra_node2[2+(tetra1*8+6)*tetra_order] = node; tetra_node2[2+(tetra1*8+7)*tetra_order] = node; } else if ( v1 == 3 && v2 == 4 ) { tetra_node2[3+(tetra1*8+3)*tetra_order] = node; tetra_node2[3+(tetra1*8+6)*tetra_order] = node; tetra_node2[3+(tetra1*8+7)*tetra_order] = node; } } return; } /******************************************************************************/ void tet_mesh_order4_refine_size ( int node_num1, int tetra_num1, int tetra_node1[], int *node_num2, int *tetra_num2, int edge_data[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_REFINE_SIZE sizes a refined order 4 tet mesh. Discussion: A refined tet mesh can be derived from an existing one by interpolating nodes at the midpoint of every edge of the mesh. The mesh is described indirectly, as the sum of individual tetrahedrons. A single physical edge may be a logical edge of any number of tetrahedrons. It is important, however, that a new node be created exactly once for each edge, assigned an index, and associated with every tetrahedron that shares this edge. This routine handles that problem. The primary amount of work occurs in sorting a list of 6 * TETRA_NUM data items, one item for every edge of every tetrahedron. Each data item records, for a given tetrahedron edge, the global indices of the two endpoints, the local indices of the two endpoints, and the index of the tetrahedron. Through careful sorting, it is possible to arrange this data in a way that allows the proper generation of the interpolated nodes. Licensing: This code is distributed under the MIT license. Modified: 25 January 2007 Author: John Burkardt Parameters: Input, int NODE_NUM1, the number of nodes in the original mesh. Input, int TETRA_NUM1, the number of tetrahedrons in the original mesh. Input, int TETRA_NODE1[4*TETRA_NUM1], the indices of the nodes in the original mesh. Output, int *NODE_NUM2, the number of nodes in the refined mesh. Output, int *TETRA_NUM2, the number of tetrahedrons in the refined mesh. Output, int EDGE_DATA[5*(6*TETRA_NUM1)], edge data. */ { int a; int b; int edge; int i; int j; int k; int l; int n1; int n1_old; int n2; int n2_old; int tetra; int tetra_order = 4; /* Step 1. From the list of nodes for tetrahedron T, of the form: (I,J,K,L) construct the six edge relations: (I,J,1,2,T) (I,K,1,3,T) (I,L,1,4,T) (J,K,2,3,T) (J,L,2,4,T) (K,L,3,4,T) In order to make matching easier, we reorder each pair of nodes into ascending order. */ for ( tetra = 0; tetra < tetra_num1; tetra++ ) { i = tetra_node1[0+tetra*tetra_order]; j = tetra_node1[1+tetra*tetra_order]; k = tetra_node1[2+tetra*tetra_order]; l = tetra_node1[3+tetra*tetra_order]; i4i4_sort_a ( i, j, &a, &b ); edge_data[0+(6*tetra)*5] = a; edge_data[1+(6*tetra)*5] = b; edge_data[2+(6*tetra)*5] = 1; edge_data[3+(6*tetra)*5] = 2; edge_data[4+(6*tetra)*5] = tetra; i4i4_sort_a ( i, k, &a, &b ); edge_data[0+(6*tetra+1)*5] = a; edge_data[1+(6*tetra+1)*5] = b; edge_data[2+(6*tetra+1)*5] = 1; edge_data[3+(6*tetra+1)*5] = 3; edge_data[4+(6*tetra+1)*5] = tetra; i4i4_sort_a ( i, l, &a, &b ); edge_data[0+(6*tetra+2)*5] = a; edge_data[1+(6*tetra+2)*5] = b; edge_data[2+(6*tetra+2)*5] = 1; edge_data[3+(6*tetra+2)*5] = 4; edge_data[4+(6*tetra+2)*5] = tetra; i4i4_sort_a ( j, k, &a, &b ); edge_data[0+(6*tetra+3)*5] = a; edge_data[1+(6*tetra+3)*5] = b; edge_data[2+(6*tetra+3)*5] = 2; edge_data[3+(6*tetra+3)*5] = 3; edge_data[4+(6*tetra+3)*5] = tetra; i4i4_sort_a ( j, l, &a, &b ); edge_data[0+(6*tetra+4)*5] = a; edge_data[1+(6*tetra+4)*5] = b; edge_data[2+(6*tetra+4)*5] = 2; edge_data[3+(6*tetra+4)*5] = 4; edge_data[4+(6*tetra+4)*5] = tetra; i4i4_sort_a ( k, l, &a, &b ); edge_data[0+(6*tetra+5)*5] = a; edge_data[1+(6*tetra+5)*5] = b; edge_data[2+(6*tetra+5)*5] = 3; edge_data[3+(6*tetra+5)*5] = 4; edge_data[4+(6*tetra+5)*5] = tetra; } /* Step 2. Perform an ascending dictionary sort on the neighbor relations. We only intend to sort on rows 1:2; the routine we call here sorts on the full column but that won't hurt us. What we need is to find all cases where tetrahedrons share an edge. By sorting the columns of the EDGE_DATA array, we will put shared edges next to each other. */ i4col_sort_a ( 5, 6*tetra_num1, edge_data ); /* Step 3. All the tetrahedrons which share an edge show up as consecutive columns with identical first two entries. Figure out how many new nodes there are, and allocate space for their coordinates. */ *node_num2 = node_num1; n1_old = -1; n2_old = -1; for ( edge = 0; edge < 6 * tetra_num1; edge++ ) { n1 = edge_data[0+edge*5]; n2 = edge_data[1+edge*5]; if ( n1 != n1_old || n2 != n2_old ) { *node_num2 = *node_num2 + 1; n1_old = n1; n2_old = n2; } } *tetra_num2 = 8 * tetra_num1; return; } /******************************************************************************/ void tet_mesh_order4_to_order10_compute ( int tetra_num, int tetra_node1[], int node_num1, double node_xyz1[], int edge_data[], int tetra_node2[], int node_num2, double node_xyz2[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_TO_ORDER10_COMPUTE computes a quadratic tet mesh from a linear one. Discussion: A quadratic (10 node) tet mesh can be derived from a linear (4 node) tet mesh by interpolating nodes at the midpoint of every edge of the mesh. The mesh is described indirectly, as the sum of individual tetrahedrons. A single physical edge may be a logical edge of any number of tetrahedrons. It is important, however, that a new node be created exactly once for each edge, assigned an index, and associated with every tetrahedron that shares this edge. This routine handles that problem. The primary amount of work occurs in sorting a list of 6 * TETRA_NUM data items, one item for every edge of every tetrahedron. Each data item records, for a given tetrahedron edge, the global indices of the two endpoints, the local indices of the two endpoints, and the index of the tetrahedron. Through careful sorting, it is possible to arrange this data in a way that allows the proper generation of the interpolated nodes. The node ordering for the quadratic tetrahedron is somewhat arbitrary. In the current scheme, the vertices are listed first, followed by the 6 midside nodes. Each midside node may be identified by the two vertices that bracket it. Thus, the node ordering may be suggested by: 1 2 3 4 (1+2) (1+3) (1+4) (2+3) (2+4) (3+4) Licensing: This code is distributed under the MIT license. Modified: 23 January 2007 Author: John Burkardt Parameters: Input, int TETRA_NUM, the number of tetrahedrons in the linear mesh. Input, int TETRA_NODE1[4*TETRA_NUM], the indices of the nodes in the linear mesh. Input, int NODE_NUM1, the number of nodes for the linear mesh. Input, double NODE_XYZ1[3*NODE_NUM1], the coordinates of the nodes that make up the linear mesh. Input, int EDGE_DATA[5*(6*TETRA_NUM)], edge data. Output, int TETRA_NODE2[10*TETRA_NUM], the indices of the nodes in the quadratic mesh. Input, int NODE_NUM2, the number of nodes for the quadratic mesh. Output, double NODE_XYZ2[3*NODE_NUM2], the coordinates of the nodes that make up the quadratic mesh. */ { int dim_num = 3; int edge; int i; int j; int n1; int n1_old; int n2; int n2_old; int node; int tetra; int tetra_order1 = 4; int tetra_order2 = 10; int v; int v1; int v2; /* Generate the index and coordinates of the new midside nodes, and update the tetradehron-node data. */ for ( j = 0; j < node_num1; j++ ) { for ( i = 0; i < dim_num; i++ ) { node_xyz2[i+j*dim_num] = node_xyz1[i+j*dim_num]; } } for ( j = 0; j < tetra_num; j++ ) { for ( i = 0; i < tetra_order1; i++ ) { tetra_node2[i+j*tetra_order2] = tetra_node1[i+j*tetra_order1]; } } node = node_num1; n1_old = -1; n2_old = -1; for ( edge = 0; edge < 6 * tetra_num; edge++ ) { /* Read the data defining the edge. */ n1 = edge_data[0+edge*5]; n2 = edge_data[1+edge*5]; /* If this edge is new, create the coordinates and index. */ if ( n1 != n1_old || n2 != n2_old ) { if ( node_num2 <= node ) { printf ( "\n" ); printf ( "TET_MESH_ORDER4_TO_ORDER10_COMPUTE - Fatal error!\n" ); printf ( " Node index exceeds NODE_NUM2.\n" ); exit ( 1 ); } for ( i = 0; i < dim_num; i++ ) { node_xyz2[i+node*dim_num] = ( node_xyz2[i+(n1-1)*dim_num] + node_xyz2[i+(n2-1)*dim_num] ) / 2.0; } node = node + 1; n1_old = n1; n2_old = n2; } /* Assign the node to the tetrahedron. */ v1 = edge_data[2+edge*5]; v2 = edge_data[3+edge*5]; /* Here is where the local ordering of the nodes is effected: */ if ( v1 == 1 && v2 == 2 ) { v = 5; } else if ( v1 == 1 && v2 == 3 ) { v = 6; } else if ( v1 == 1 && v2 == 4 ) { v = 7; } else if ( v1 == 2 && v2 == 3 ) { v = 8; } else if ( v1 == 2 && v2 == 4 ) { v = 9; } else if ( v1 == 3 && v2 == 4 ) { v = 10; } tetra = edge_data[4+edge*5]; tetra_node2[v-1+tetra*tetra_order2] = node; } return; } /******************************************************************************/ void tet_mesh_order4_to_order10_size ( int tetra_num, int tetra_node1[], int node_num1, int edge_data[], int *node_num2 ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER4_TO_ORDER10_SIZE sizes a quadratic tet mesh from a linear one. Discussion: A quadratic (10 node) tet mesh can be derived from a linear (4 node) tet mesh by interpolating nodes at the midpoint of every edge of the mesh. The mesh is described indirectly, as the sum of individual tetrahedrons. A single physical edge may be a logical edge of any number of tetrahedrons. It is important, however, that a new node be created exactly once for each edge, assigned an index, and associated with every tetrahedron that shares this edge. This routine handles that problem. The primary amount of work occurs in sorting a list of 6 * TETRA_NUM data items, one item for every edge of every tetrahedron. Each data item records, for a given tetrahedron edge, the global indices of the two endpoints, the local indices of the two endpoints, and the index of the tetrahedron. Through careful sorting, it is possible to arrange this data in a way that allows the proper generation of the interpolated nodes. Licensing: This code is distributed under the MIT license. Modified: 09 December 2006 Author: John Burkardt Parameters: Input, int TETRA_NUM, the number of tetrahedrons in the linear mesh. Input, int TETRA_NODE1[4*TETRA_NUM], the indices of the nodes in the linear mesh. Input, int NODE_NUM1, the number of nodes for the linear mesh. Output, int EDGE_DATA[5*(6*TETRA_NUM)], edge data. Output, int *NODE_NUM2, the number of nodes for the quadratic mesh. */ { int a; int b; int edge; int i; int j; int k; int l; int n1; int n1_old; int n2; int n2_old; int tetra; int tetra_order1 = 4; /* Step 1. From the list of nodes for tetrahedron T, of the form: (I,J,K,L) construct the six edge relations: (I,J,1,2,T) (I,K,1,3,T) (I,L,1,4,T) (J,K,2,3,T) (J,L,2,4,T) (K,L,3,4,T) In order to make matching easier, we reorder each pair of nodes into ascending order. */ for ( tetra = 0; tetra < tetra_num; tetra++ ) { i = tetra_node1[0+tetra*tetra_order1]; j = tetra_node1[1+tetra*tetra_order1]; k = tetra_node1[2+tetra*tetra_order1]; l = tetra_node1[3+tetra*tetra_order1]; i4i4_sort_a ( i, j, &a, &b ); edge_data[0+(6*tetra)*5] = a; edge_data[1+(6*tetra)*5] = b; edge_data[2+(6*tetra)*5] = 1; edge_data[3+(6*tetra)*5] = 2; edge_data[4+(6*tetra)*5] = tetra; i4i4_sort_a ( i, k, &a, &b ); edge_data[0+(6*tetra+1)*5] = a; edge_data[1+(6*tetra+1)*5] = b; edge_data[2+(6*tetra+1)*5] = 1; edge_data[3+(6*tetra+1)*5] = 3; edge_data[4+(6*tetra+1)*5] = tetra; i4i4_sort_a ( i, l, &a, &b ); edge_data[0+(6*tetra+2)*5] = a; edge_data[1+(6*tetra+2)*5] = b; edge_data[2+(6*tetra+2)*5] = 1; edge_data[3+(6*tetra+2)*5] = 4; edge_data[4+(6*tetra+2)*5] = tetra; i4i4_sort_a ( j, k, &a, &b ); edge_data[0+(6*tetra+3)*5] = a; edge_data[1+(6*tetra+3)*5] = b; edge_data[2+(6*tetra+3)*5] = 2; edge_data[3+(6*tetra+3)*5] = 3; edge_data[4+(6*tetra+3)*5] = tetra; i4i4_sort_a ( j, l, &a, &b ); edge_data[0+(6*tetra+4)*5] = a; edge_data[1+(6*tetra+4)*5] = b; edge_data[2+(6*tetra+4)*5] = 2; edge_data[3+(6*tetra+4)*5] = 4; edge_data[4+(6*tetra+4)*5] = tetra; i4i4_sort_a ( k, l, &a, &b ); edge_data[0+(6*tetra+5)*5] = a; edge_data[1+(6*tetra+5)*5] = b; edge_data[2+(6*tetra+5)*5] = 3; edge_data[3+(6*tetra+5)*5] = 4; edge_data[4+(6*tetra+5)*5] = tetra; } /* Step 2. Perform an ascending dictionary sort on the neighbor relations. We only intend to sort on rows 1:2; the routine we call here sorts on the full column but that won't hurt us. What we need is to find all cases where tetrahedrons share an edge. By sorting the columns of the EDGE_DATA array, we will put shared edges next to each other. */ i4col_sort_a ( 5, 6*tetra_num, edge_data ); /* Step 3. All the tetrahedrons which share an edge show up as consecutive columns with identical first two entries. Figure out how many new nodes there are, and allocate space for their coordinates. */ *node_num2 = node_num1; n1_old = -1; n2_old = -1; for ( edge = 0; edge < 6 * tetra_num; edge++ ) { n1 = edge_data[0+edge*5]; n2 = edge_data[1+edge*5]; if ( n1 != n1_old || n2 != n2_old ) { *node_num2 = *node_num2 + 1; n1_old = n1; n2_old = n2; } } return; } /******************************************************************************/ void tet_mesh_order10_adj_count ( int node_num, int tet_num, int tet_node[], int *adj_num, int adj_row[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER10_ADJ_COUNT counts the number of nodal adjacencies. Discussion: Assuming that the tet mesh is to be used in a finite element computation, we declare that two distinct nodes are "adjacent" if and only if they are both included in some tetrahedron. It is the purpose of this routine to determine the number of such adjacency relationships. The initial count gets only the (I,J) relationships, for which node I is strictly less than node J. This value is doubled to account for symmetry. Licensing: This code is distributed under the MIT license. Modified: 08 March 2013 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int TET_NUM, the number of tetrahedrons. Input, int TET_NODE[10*TET_NUM], the indices of the nodes. Output, int *ADJ_NUM, the total number of adjacency relationships, Output, int ADJ_ROW[NODE_NUM+1], the ADJ pointer array. */ { int i; int j; int k; int l; int node; int *pair; int pair_num; int pair_unique_num; /* Each order 10 tetrahedron defines 45 adjacency pairs. */ pair = ( int * ) malloc ( 2 * 45 * tet_num * sizeof ( int ) ); k = 0; for ( i = 0; i < 9; i++ ) { for ( j = i + 1; j < 10; j++ ) { for ( l = 0; l < tet_num; l++ ) { pair[0+(k*tet_num+l)*2] = tet_node[i+l*10]; pair[1+(k*tet_num+l)*2] = tet_node[j+l*10]; } k = k + 1; } } /* Force the nodes of each pair to be listed in ascending order. */ pair_num = 45 * tet_num; i4col_sort2_a ( 2, pair_num, pair ); /* Rearrange the columns in ascending order. */ i4col_sort_a ( 2, pair_num, pair ); /* Get the number of unique columns. */ pair_unique_num = i4col_sorted_unique_count ( 2, pair_num, pair ); /* The number of adjacencies is TWICE this value, plus the number of nodes. */ *adj_num = 2 * pair_unique_num; /* Now set up the ADJ_ROW counts. */ for ( node = 0; node < node_num; node++ ) { adj_row[node] = 0; } for ( k = 0; k < pair_num; k++ ) { if ( 0 < k ) { if ( pair[0+(k-1)*2] == pair[0+k*2] && pair[1+(k-1)*2] == pair[1+k*2] ) { continue; } } i = pair[0+k*2]; j = pair[1+k*2]; adj_row[i-1] = adj_row[i-1] + 1; adj_row[j-1] = adj_row[j-1] + 1; } /* We used ADJ_ROW to count the number of entries in each row. Convert it to pointers into the ADJ array. */ for ( node = node_num-1; 0 <= node; node-- ) { adj_row[node] = adj_row[node+1]; } adj_row[0] = 1; for ( node = 1; node <= node_num; node++ ) { adj_row[node] = adj_row[node-1] + adj_row[i]; } free ( pair ); return; } /******************************************************************************/ int *tet_mesh_order10_adj_set ( int node_num, int tet_num, int tet_node[], int adj_num, int adj_row[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER10_ADJ_SET sets the nodal adjacency matrix. Discussion: A compressed format is used for the nodal adjacency matrix. It is assumed that we know ADJ_NUM, the number of adjacency entries and the ADJ_ROW array, which keeps track of the list of slots in ADJ where we can store adjacency information for each row. We essentially repeat the work of TET_MESH_ORDER4_ADJ_COUNT, but now we have a place to store the adjacency information. A copy of the ADJ_ROW array is useful, as we can use it to keep track of the next available entry in ADJ for adjacencies associated with a given row. Licensing: This code is distributed under the MIT license. Modified: 08 March 2013 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int TET_NUM, the number of tetrahedrons. Input, int TET_NODE[10*TET_NUM], the indices of the nodes. Input, int ADJ_NUM, the total number of adjacency relationships, Input, int ADJ_ROW[NODE_NUM+1], the ADJ pointer array. Output, int TET_MESH_ORDER4_ADJ_SET[ADJ_NUM], the adjacency information. */ { int *adj; int *adj_row_copy; int i; int j; int k; int l; int node; int *pair; int pair_num; /* Each order 10 tetrahedron defines 45 adjacency pairs. */ pair = ( int * ) malloc ( 2 * 45 * tet_num * sizeof ( int ) ); k = 0; for ( i = 0; i < 9; i++ ) { for ( j = i + 1; j < 10; j++ ) { for ( l = 0; l < tet_num; l++ ) { pair[0+(k*tet_num+l)*2] = tet_node[i+l*10]; pair[1+(k*tet_num+l)*2] = tet_node[j+l*10]; } k = k + 1; } } /* Force the nodes of each pair to be listed in ascending order. */ pair_num = 45 * tet_num; i4col_sort2_a ( 2, pair_num, pair ); /* Rearrange the columns in ascending order. */ i4col_sort_a ( 2, pair_num, pair ); /* Mark all entries of ADJ so we will know later if we missed one. */ adj = ( int * ) malloc ( adj_num * sizeof ( int ) ); for ( i = 0; i < adj_num; i++ ) { adj[i] = -1; } /* Copy the ADJ_ROW array and use it to keep track of the next free entry for each row. */ adj_row_copy = ( int * ) malloc ( node_num * sizeof ( int ) ); for ( node = 0; node < node_num; node++ ) { adj_row_copy[node] = adj_row[node]; } /* Now set up the ADJ_ROW counts. */ for ( k = 0; k < pair_num; k++ ) { if ( 0 < k ) { if ( pair[0+(k-1)*2] == pair[0+k*2] && pair[1+(k-1)*2] == pair[1+k*2] ) { continue; } } i = pair[0+k*2]; j = pair[1+k*2]; adj[adj_row_copy[i]] = j; adj_row_copy[i] = adj_row_copy[i] + 1; adj[adj_row_copy[j]] = i; adj_row_copy[j] = adj_row_copy[j] + 1; } free ( adj_row_copy ); free ( pair ); return adj; } /******************************************************************************/ void tet_mesh_order10_example_set ( int node_num, int tetra_num, double node_xyz[], int tetra_node[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER10_EXAMPLE_SET sets an example quadratic tet mesh. Licensing: This code is distributed under the MIT license. Modified: 03 August 2009 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, int TETRA_NUM, the number of tetrahedrons. Output, double NODE_XYZ[3*NODE_NUM], the node coordinates. Output, int TETRA_NODE[10*TETRA_NUM], the nodes forming each tet. */ { int i; int j; double node_xyz_save[3*27] = { 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 0.0, 0.5, 0.0, 0.5, 0.0, 0.0, 0.5, 0.5, 0.5, 0.0, 0.0, 0.0, 0.5, 1.0, 0.5, 0.0, 0.5, 0.5, 0.0, 1.0, 0.0, 1.0, 0.5, 0.5, 0.5, 0.0, 0.5, 1.0, 0.0, 0.5, 0.5, 0.5, 0.5, 0.5, 1.0, 0.5, 1.0, 0.5, 0.5, 1.0, 1.0, 1.0, 0.0, 0.5, 1.0, 0.5, 0.0, 1.0, 0.5, 0.5, 1.0, 0.5, 1.0, 1.0, 1.0, 0.5 }; int tetra_node_save[10*6] = { 4, 3, 5, 1, 16, 19, 17, 11, 10, 12, 4, 2, 5, 1, 13, 19, 14, 11, 9, 12, 4, 7, 3, 5, 21, 16, 18, 19, 24, 17, 4, 7, 8, 5, 21, 22, 27, 19, 24, 25, 4, 6, 2, 5, 20, 13, 15, 19, 23, 14, 4, 6, 8, 5, 20, 22, 26, 19, 23, 25 }; for ( j = 0; j < node_num; j++ ) { for ( i = 0; i < 3; i++ ) { node_xyz[i+j*3] = node_xyz_save[i+j*3]; } } for ( j = 0; j < tetra_num; j++ ) { for ( i = 0; i < 10; i++ ) { tetra_node[i+j*10] = tetra_node_save[i+j*10] - 1; } } return; } /******************************************************************************/ void tet_mesh_order10_example_size ( int *node_num, int *tetra_num ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER10_EXAMPLE_SIZE sizes an example quadratic tet mesh. Licensing: This code is distributed under the MIT license. Modified: 03 August 2009 Author: John Burkardt Parameters: Output, int *NODE_NUM, the number of nodes. Output, int *TETRA_NUM, the number of tetrahedrons. */ { *node_num = 27; *tetra_num = 6; return; } /******************************************************************************/ void tet_mesh_order10_to_order4_compute ( int tetra_num1, int tetra_node1[], int tetra_num2, int tetra_node2[] ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER10_TO_ORDER4_COMPUTE linearizes a quadratic tet mesh. Discussion: A quadratic tet mesh is assumed to consist of 10-node tetrahedrons. This routine rearranges the information so as to define a 4-node tet mesh. The same nodes are used, but there are 8 times as many tetrahedrons. The node ordering for the quadratic tetrahedron is somewhat arbitrary. In the current scheme, the vertices are listed first, followed by the 6 midside nodes. Each midside node may be identified by the two vertices that bracket it. Thus, the node ordering may be suggested by: 1 2 3 4 (1+2) (1+3) (1+4) (2+3) (2+4) (3+4) Licensing: This code is distributed under the MIT license. Modified: 22 January 2007 Author: John Burkardt Reference: Anwei Liu, Barry Joe, Quality Local Refinement of Tetrahedral Meshes Based on 8-Subtetrahedron Subdivision, Mathematics of Computation, Volume 65, Number 215, July 1996, pages 1183-1200. Parameters: Input, int TETRA_NUM1, the number of tetrahedrons in the quadratic tet mesh. Input, int TETRA_NODE1[10*TETRA_NUM1], the indices of the nodes that made up the quadratic mesh. Input, int TETRA_NUM2, the number of tetrahedrons in the linear tet mesh. TETRA_NUM2 = 8 * TETRA_NUM1. Output, int TETRA_NODE2[4*TETRA_NUM2], the indices of the nodes that make up the linear mesh. */ { int n1; int n2; int n3; int n4; int n5; int n6; int n7; int n8; int n9; int nx; int tetra1; int tetra2; tetra2 = 0; for ( tetra1 = 0; tetra1 < tetra_num1; tetra1++ ) { n1 = tetra_node1[0+tetra1*10]; n2 = tetra_node1[1+tetra1*10]; n3 = tetra_node1[2+tetra1*10]; n4 = tetra_node1[3+tetra1*10]; n5 = tetra_node1[4+tetra1*10]; n6 = tetra_node1[5+tetra1*10]; n7 = tetra_node1[6+tetra1*10]; n8 = tetra_node1[7+tetra1*10]; n9 = tetra_node1[8+tetra1*10]; nx = tetra_node1[9+tetra1*10]; tetra_node2[0+tetra2*4] = n1; tetra_node2[1+tetra2*4] = n5; tetra_node2[2+tetra2*4] = n6; tetra_node2[3+tetra2*4] = n7; tetra2 = tetra2 + 1; tetra_node2[0+tetra2*4] = n2; tetra_node2[1+tetra2*4] = n5; tetra_node2[2+tetra2*4] = n8; tetra_node2[3+tetra2*4] = n9; tetra2 = tetra2 + 1; tetra_node2[0+tetra2*4] = n3; tetra_node2[1+tetra2*4] = n6; tetra_node2[2+tetra2*4] = n8; tetra_node2[3+tetra2*4] = n9; tetra2 = tetra2 + 1; tetra_node2[0+tetra2*4] = n4; tetra_node2[1+tetra2*4] = n7; tetra_node2[2+tetra2*4] = n9; tetra_node2[3+tetra2*4] = nx; tetra2 = tetra2 + 1; tetra_node2[0+tetra2*4] = n5; tetra_node2[1+tetra2*4] = n6; tetra_node2[2+tetra2*4] = n7; tetra_node2[3+tetra2*4] = n9; tetra2 = tetra2 + 1; tetra_node2[0+tetra2*4] = n5; tetra_node2[1+tetra2*4] = n6; tetra_node2[2+tetra2*4] = n8; tetra_node2[3+tetra2*4] = n9; tetra2 = tetra2 + 1; tetra_node2[0+tetra2*4] = n6; tetra_node2[1+tetra2*4] = n7; tetra_node2[2+tetra2*4] = n9; tetra_node2[3+tetra2*4] = nx; tetra2 = tetra2 + 1; tetra_node2[0+tetra2*4] = n6; tetra_node2[1+tetra2*4] = n8; tetra_node2[2+tetra2*4] = n9; tetra_node2[3+tetra2*4] = nx; tetra2 = tetra2 + 1; } return; } /******************************************************************************/ void tet_mesh_order10_to_order4_size ( int node_num1, int tetra_num1, int *node_num2, int *tetra_num2 ) /******************************************************************************/ /* Purpose: TET_MESH_ORDER10_TO_ORDER4_SIZE sizes a linear tet mesh from a quadratic one. Discussion: A linear (4 node) tet mesh can be derived from a quadratic (10 node) tet mesh using the same set of nodes, but reassigning the nodes of each quadratic tet among 8 linear subtets. This routine returns the number of nodes and tetrahedra in the linear mesh. Licensing: This code is distributed under the MIT license. Modified: 06 December 2006 Author: John Burkardt Reference: Anwei Liu, Barry Joe, Quality Local Refinement of Tetrahedral Meshes Based on 8-Subtetrahedron Subdivision, Mathematics of Computation, Volume 65, Number 215, July 1996, pages 1183-1200. Parameters: Input, int NODE_NUM1, the number of nodes in the quadratic mesh. Input, int TETRA_NUM1, the number of tetrahedrons in the quadratic mesh. Output, int *NODE_NUM2, the number of nodes for the linear mesh. Output, int *TETRA_NUM2, the number of tetrahedrons in the linear mesh. */ { *node_num2 = node_num1; *tetra_num2 = 8 * tetra_num1; return; } /******************************************************************************/ void tet_mesh_quad ( int node_num, double node_xyz[], int tetra_order, int tetra_num, int tetra_node[], void quad_fun ( int n, double xyz_vec[], double fvec[] ), int quad_num, double quad_xyz[], double quad_w[], double *quad_value, double *region_volume ) /******************************************************************************/ /* Purpose: TET_MESH_QUAD approximates an integral over a tet mesh. Discussion: The routine will accept tetrahedral meshes of order higher than 4. However, only the first four nodes (the vertices) of each tetrahedron will be used. This will still produce correct results for higher order tet meshes, as long as the sides of each tetrahedron are flat (linear). We assume that the vertices of each tetrahedron are listed first in the description of higher order tetrahedrons. The approximation of the integral is made using a quadrature rule defined on the unit tetrahedron, and supplied by the user. The user also supplies the name of a subroutine, here called "QUAD_FUN", which evaluates the integrand at a set of points. The form is: void quad_fun ( int n, double xyz_vec[3*n], double f_vec[n] ) and it returns in each entry F_VEC(1:N), the value of the integrand at XYZ_VEC(1:3,1:N). Licensing: This code is distributed under the MIT license. Modified: 06 December 2006 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes in the tet mesh. Input, double NODE_XYZ[3*NODE_NUM], the coordinates of the nodes. Input, int TETRA_ORDER, the order of tetrahedrons in the tet mesh. Input, int TETRA_NUM, the number of tetrahedrons in the tet mesh. Input, int TETRA_NODE[TETRA_ORDER*TETRA_NUM], indices of the nodes. Input, void QUAD_FUN ( int N, double XYZ_VEC[3*N], F_VEC[N] ), the name of the routine that evaluates the integrand. Input, int QUAD_NUM, the order of the quadrature rule. Input, double QUAD_XYZ[3*QUAD_NUM], the abscissas of the quadrature rule, in the unit tetrahedron. Input, double QUAD_W[QUAD_NUM], the weights of the quadrature rule. Output, double *QUAD_VALUE, the estimate of the integral of F(X,Y) over the region covered by the tet mesh. Output, double *REGION_VOLUME, the volume of the region. */ { int i; int j; int quad; double quad_f[quad_num]; double quad2_xyz[3*quad_num]; double temp; int tet; double tetra_volume; double tetra_xyz[3*4]; *quad_value = 0.0; *region_volume = 0.0; for ( tet = 0; tet < tetra_num; tet++ ) { for ( j = 0; j < 4; j++ ) { for ( i = 0; i < 3; i++ ) { tetra_xyz[i+j*3] = node_xyz[i+(tetra_node[j+tet*4]-1)*3]; } } tetra_volume = tetrahedron_volume ( tetra_xyz ); tetrahedron_order4_reference_to_physical ( tetra_xyz, quad_num, quad_xyz, quad2_xyz ); quad_fun ( quad_num, quad2_xyz, quad_f ); temp = 0.0; for ( quad = 0; quad < quad_num; quad++ ) { temp = temp + quad_w[quad] * quad_f[quad]; } *quad_value = *quad_value + tetra_volume * temp; *region_volume = *region_volume + tetra_volume; } return; } /******************************************************************************/ void tet_mesh_quality1 ( int node_num, double node_xyz[], int tetra_order, int tetra_num, int tetra_node[], double *value_min, double *value_mean, double *value_max, double *value_var ) /******************************************************************************/ /* Purpose: TET_MESH_QUALITY1 returns a tet mesh quality factor. Discussion: The tet mesh quality measure is the minimum of the corresponding tetrahedron quality measure, over all tetrahedrons in the tet mesh. This routine is designed for a 4-node tet mesh. It can handle a 10-node tet mesh, but it simply ignores the extra nodes. Licensing: This code is distributed under the MIT license. Modified: 27 October 2005 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, double NODE_XYZ[3*NODE_NUM], the coordinates of the nodes. Input, int TETRA_ORDER, the order of the tetrahedrons. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[TETRA_ORDER*TETRA_NUM], the indices of the nodes. Output, double *VALUE_MIN, *VALUE_MEAN, *VALUE_MAX, *VALUE_VAR, the minimum, mean, maximum and variance of the quality measure. */ { # define DIM_NUM 3 int i; int j; int node; int tetra; double tetrahedron[DIM_NUM*4]; double *tetrahedron_quality; tetrahedron_quality = ( double * ) malloc ( tetra_num * sizeof ( double ) ); for ( tetra = 0; tetra < tetra_num; tetra++ ) { for ( j = 0; j < 4; j++ ) { node = tetra_node[j+tetra*tetra_order]; for ( i = 0; i < DIM_NUM; i++ ) { tetrahedron[i+j*DIM_NUM] = node_xyz[i+(node-1)*DIM_NUM]; } } tetrahedron_quality[tetra] = tetrahedron_quality1_3d ( tetrahedron ); } *value_max = r8vec_max ( tetra_num, tetrahedron_quality ); *value_min = r8vec_min ( tetra_num, tetrahedron_quality ); *value_mean = r8vec_mean ( tetra_num, tetrahedron_quality ); *value_var = r8vec_variance ( tetra_num, tetrahedron_quality ); free ( tetrahedron_quality ); return; # undef DIM_NUM } /******************************************************************************/ void tet_mesh_quality2 ( int node_num, double node_xyz[], int tetra_order, int tetra_num, int tetra_node[], double *value_min, double *value_mean, double *value_max, double *value_var ) /******************************************************************************/ /* Purpose: TET_MESH_QUALITY2 returns a tet mesh quality factor. Discussion: The tet mesh quality measure is the minimum of the corresponding tetrahedron quality measure, over all tetrahedrons in the tet mesh. This routine is designed for a 4-node tet mesh. It can handle a 10-node tet mesh, but it simply ignores the extra nodes. Licensing: This code is distributed under the MIT license. Modified: 27 October 2005 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, double NODE_XYZ[3*NODE_NUM], the coordinates of the nodes. Input, int TETRA_ORDER, the order of the tetrahedrons. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[TETRA_ORDER*TETRA_NUM], the indices of the nodes. Output, double *VALUE_MIN, *VALUE_MEAN, *VALUE_MAX, *VALUE_VAR, the minimum, mean, maximum and variance of the quality measure. */ { # define DIM_NUM 3 int i; int j; int node; int tetra; double tetrahedron[DIM_NUM*4]; double *tetrahedron_quality; tetrahedron_quality = ( double * ) malloc ( tetra_num * sizeof ( double ) ); for ( tetra = 0; tetra < tetra_num; tetra++ ) { for ( j = 0; j < 4; j++ ) { node = tetra_node[j+tetra*tetra_order]; for ( i = 0; i < DIM_NUM; i++ ) { tetrahedron[i+j*DIM_NUM] = node_xyz[i+(node-1)*DIM_NUM]; } } tetrahedron_quality[tetra] = tetrahedron_quality2_3d ( tetrahedron ); } *value_max = r8vec_max ( tetra_num, tetrahedron_quality ); *value_min = r8vec_min ( tetra_num, tetrahedron_quality ); *value_mean = r8vec_mean ( tetra_num, tetrahedron_quality ); *value_var = r8vec_variance ( tetra_num, tetrahedron_quality ); free ( tetrahedron_quality ); return; # undef DIM_NUM } /******************************************************************************/ void tet_mesh_quality3 ( int node_num, double node_xyz[], int tetra_order, int tetra_num, int tetra_node[], double *value_min, double *value_mean, double *value_max, double *value_var ) /******************************************************************************/ /* Purpose: TET_MESH_QUALITY3 returns a tet mesh quality factor. Discussion: The tet mesh quality measure is the minimum of the corresponding tetrahedron quality measure, over all tetrahedrons in the tet mesh. This routine is designed for a 4-node tet mesh. It can handle a 10-node tet mesh, but it simply ignores the extra nodes. Licensing: This code is distributed under the MIT license. Modified: 27 October 2005 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, double NODE_XYZ[3*NODE_NUM], the coordinates of the nodes. Input, int TETRA_ORDER, the order of the tetrahedrons. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[TETRA_ORDER*TETRA_NUM], the indices of the nodes. Output, double *VALUE_MIN, *VALUE_MEAN, *VALUE_MAX, *VALUE_VAR, the minimum, mean, maximum and variance of the quality measure. */ { # define DIM_NUM 3 int i; int j; int node; int tetra; double tetrahedron[DIM_NUM*4]; double *tetrahedron_quality; tetrahedron_quality = ( double * ) malloc ( tetra_num * sizeof ( double ) ); for ( tetra = 0; tetra < tetra_num; tetra++ ) { for ( j = 0; j < 4; j++ ) { node = tetra_node[j+tetra*tetra_order]; for ( i = 0; i < DIM_NUM; i++ ) { tetrahedron[i+j*DIM_NUM] = node_xyz[i+(node-1)*DIM_NUM]; } } tetrahedron_quality[tetra] = tetrahedron_quality3_3d ( tetrahedron ); } *value_max = r8vec_max ( tetra_num, tetrahedron_quality ); *value_min = r8vec_min ( tetra_num, tetrahedron_quality ); *value_mean = r8vec_mean ( tetra_num, tetrahedron_quality ); *value_var = r8vec_variance ( tetra_num, tetrahedron_quality ); free ( tetrahedron_quality ); return; # undef DIM_NUM } /******************************************************************************/ void tet_mesh_quality4 ( int node_num, double node_xyz[], int tetra_order, int tetra_num, int tetra_node[], double *value_min, double *value_mean, double *value_max, double *value_var ) /******************************************************************************/ /* Purpose: TET_MESH_QUALITY4 returns a tet mesh quality factor. Discussion: The tet mesh quality measure is the minimum of the corresponding tetrahedron quality measure, over all tetrahedrons in the tet mesh. This routine is designed for a 4-node tet mesh. It can handle a 10-node tet mesh, but it simply ignores the extra nodes. Licensing: This code is distributed under the MIT license. Modified: 27 October 2005 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, double NODE_XYZ[3*NODE_NUM], the coordinates of the nodes. Input, int TETRA_ORDER, the order of the tetrahedrons. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[TETRA_ORDER*TETRA_NUM], the indices of the nodes. Output, double *VALUE_MIN, *VALUE_MEAN, *VALUE_MAX, *VALUE_VAR, the minimum, mean, maximum and variance of the quality measure. */ { # define DIM_NUM 3 int i; int j; int node; int tetra; double tetrahedron[DIM_NUM*4]; double *tetrahedron_quality; tetrahedron_quality = ( double * ) malloc ( tetra_num * sizeof ( double ) ); for ( tetra = 0; tetra < tetra_num; tetra++ ) { for ( j = 0; j < 4; j++ ) { node = tetra_node[j+tetra*tetra_order]; for ( i = 0; i < DIM_NUM; i++ ) { tetrahedron[i+j*DIM_NUM] = node_xyz[i+(node-1)*DIM_NUM]; } } tetrahedron_quality[tetra] = tetrahedron_quality4_3d ( tetrahedron ); } *value_max = r8vec_max ( tetra_num, tetrahedron_quality ); *value_min = r8vec_min ( tetra_num, tetrahedron_quality ); *value_mean = r8vec_mean ( tetra_num, tetrahedron_quality ); *value_var = r8vec_variance ( tetra_num, tetrahedron_quality ); free ( tetrahedron_quality ); return; # undef DIM_NUM } /******************************************************************************/ void tet_mesh_quality5 ( int node_num, double node_xyz[], int tetra_order, int tetra_num, int tetra_node[], double *value_min, double *value_mean, double *value_max, double *value_var ) /******************************************************************************/ /* Purpose: TET_MESH_QUALITY5 returns a tet mesh quality factor. Discussion: The tet mesh quality measure is the ratio of the minimum tetrahedron volume to the maximum tetrahedron volume. This routine is designed for a 4-node tet mesh. It can handle a 10-node tet mesh, but it simply ignores the extra nodes. Licensing: This code is distributed under the MIT license. Modified: 27 October 2005 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, double NODE_XYZ[3*NODE_NUM], the coordinates of the nodes. Input, int TETRA_ORDER, the order of the tetrahedrons. Input, int TETRA_NUM, the number of tetrahedrons. Input, int TETRA_NODE[TETRA_ORDER*TETRA_NUM], the indices of the nodes. Output, double *VALUE_MIN, *VALUE_MEAN, *VALUE_MAX, *VALUE_VAR, the minimum, mean, maximum and variance of the quality measure. */ { # define DIM_NUM 3 int i; int j; int node; int tetra; double tetrahedron[DIM_NUM*4]; double *tetrahedron_quality; double volume_max; tetrahedron_quality = ( double * ) malloc ( tetra_num * sizeof ( double ) ); for ( tetra = 0; tetra < tetra_num; tetra++ ) { for ( j = 0; j < 4; j++ ) { node = tetra_node[j+tetra*tetra_order]; for ( i = 0; i < DIM_NUM; i++ ) { tetrahedron[i+j*DIM_NUM] = node_xyz[i+(node-1)*DIM_NUM]; } } tetrahedron_quality[tetra] = tetrahedron_volume ( tetrahedron ); } volume_max = r8vec_max ( tetra_num, tetrahedron_quality ); for ( tetra = 0; tetra < tetra_num; tetra++ ) { tetrahedron_quality[tetra] = tetrahedron_quality[tetra] / volume_max; } *value_max = r8vec_max ( tetra_num, tetrahedron_quality ); *value_min = r8vec_min ( tetra_num, tetrahedron_quality ); *value_mean = r8vec_mean ( tetra_num, tetrahedron_quality ); *value_var = r8vec_variance ( tetra_num, tetrahedron_quality ); free ( tetrahedron_quality ); return; # undef DIM_NUM } /******************************************************************************/ int tet_mesh_search_delaunay ( int node_num, double node_xyz[], int tet_order, int tet_num, int tet_node[], int tet_neighbor[], double p[], int *face, int *step_num ) /******************************************************************************/ /* Purpose: TET_MESH_SEARCH_DELAUNAY searches a Delaunay tet mesh for a point. Discussion: The algorithm "walks" from one tetrahedron to its neighboring tetrahedron, and so on, until a tetrahedron 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 tetrahedron. If all 4 quantities are positive, the point is contained in the tetrahedron. If the I-th coordinate is negative, then P 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 tet mesh, the search is guaranteed to terminate. For other meshes, a continue may occur. Note the surprising fact that, even for a Delaunay tet mesh of a set of nodes, the nearest node to P need not be one of the vertices of the tetrahedron containing P. The code can be called for tet meshes of any order, but only the first 4 nodes in each tetrahedron are considered. Thus, if higher order tetrahedrons are used, and the extra nodes are intended to give the tetrahedron 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: 19 August 2009 Author: 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_XYZ[3*NODE_NUM], the coordinates of the nodes. Input, int TET_ORDER, the order of the tetrahedrons. Input, int TET_NUM, the number of tetrahedrons. Input, int TET_NODE[TET_ORDER*TET_NUM], the nodes that make up each tetrahedron. Input, int TET_NEIGHBOR[4*TET_NUM], the tetrahedron neighbor list. Input, double P[3], the coordinates of a point. Output, int *FACE, indicates the position of the point P in face TET_INDEX: 0, the interior or boundary of the tetrahedron; -1, outside the convex hull of the tet mesh, past face 1; -2, outside the convex hull of the tet mesh, past face 2; -3, outside the convex hull of the tet mesh, past face 3. -4, outside the convex hull of the tet mesh, past face 4. Output, int *STEP_NUM, the number of steps taken. Output, int TET_MESH_SEARCH_DELAUNAY, the index of the tetrahedron where the search ended. If a cycle occurred, then -1 is returned. */ { double *alpha; int i; int j; int k; int tet_index; double tet_xyz[3*4]; static int tet_index_save = -1; /* If possible, start with the previous successful value of TET_INDEX. */ if ( tet_index_save < 1 || tet_num < tet_index_save ) { tet_index = ( tet_num + 1 ) / 2; } else { tet_index = tet_index_save; } *step_num = -1; *face = 0; for ( ; ; ) { *step_num = *step_num + 1; if ( tet_num < *step_num ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TET_MESH_SEARCH_DELAUNAY - Fatal error!\n" ); fprintf ( stderr, " The algorithm seems to be cycling.\n" ); tet_index = -1; *face = -1; exit ( 1 ); } for ( j = 0; j < 4; j++ ) { k = tet_node[j+tet_index*4]; for ( i = 0; i < 3; i++ ) { tet_xyz[i+j*3] = node_xyz[i+k*3]; } } alpha = tetrahedron_barycentric ( tet_xyz, p ); /* If the barycentric coordinates are all positive, then the point is inside the tetrahedron and we're done. */ if ( 0.0 <= alpha[0] && 0.0 <= alpha[1] && 0.0 <= alpha[2] && 0.0 <= alpha[3] ) { break; } /* At least one barycentric coordinate is negative. If there is a negative barycentric coordinate for which there exists an opposing tetrahedron neighbor closer to the point, move to that tetrahedron. */ if ( alpha[0] < 0.0 && 0 < tet_neighbor[0+tet_index*4] ) { tet_index = tet_neighbor[0+tet_index*4]; continue; } else if ( alpha[1] < 0.0 && 0 < tet_neighbor[1+tet_index*4] ) { tet_index = tet_neighbor[1+tet_index*4]; continue; } else if ( alpha[2] < 0.0 && 0 < tet_neighbor[2+tet_index*4] ) { tet_index = tet_neighbor[2+tet_index*4]; continue; } else if ( alpha[3] < 0.0 && 0 < tet_neighbor[3+tet_index*4] ) { tet_index = tet_neighbor[3+tet_index*4]; continue; } /* All negative barycentric coordinates correspond to vertices opposite faces on the convex hull. Note the face and exit. */ if ( alpha[0] < 0.0 ) { *face = -1; break; } else if ( alpha[1] < 0.0 ) { *face = -2; break; } else if ( alpha[2] < 0.0 ) { *face = -3; break; } else if ( alpha[3] < 0.0 ) { *face = -4; break; } } tet_index_save = tet_index; return tet_index; } /******************************************************************************/ int tet_mesh_search_naive ( int node_num, double node_xyz[], int tet_order, int tet_num, int tet_node[], double p[], int *step_num ) /******************************************************************************/ /* Purpose: TET_MESH_SEARCH_NAIVE naively searches a tet mesh. Discussion: The algorithm simply checks each tetrahedron to see if point P is contained in it. Licensing: This code is distributed under the MIT license. Modified: 19 August 2009 Author: John Burkardt Parameters: Input, int NODE_NUM, the number of nodes. Input, double NODE_XYZ[3*NODE_NUM], the coordinates of the nodes. Input, int TET_ORDER, the order of the tetrahedrons. Input, int TET_NUM, the number of tetrahedrons in the mesh. Input, int TET_NODE[TET_ORDER*TET_NUM], the nodes that make up each tetrahedron. Input, double P[3], the coordinates of a point. Output, int TET_MESH_ORDER4_SEARCH_NAIE, the index of the tetrahedron where the search ended, or -1 if no tetrahedron was found containing the point. Output, int *STEP_NUM, the number of tetrahedrons examined.*/ { double *alpha; int i; int j; int tet; int tet_index; double tet_xyz[3*4]; tet_index = -1; *step_num = 0; for ( tet = 0; tet < tet_num; tet++ ) { for ( j = 0; j < 4; j++ ) { for ( i = 0; i < 3; i++ ) { tet_xyz[i+j*3] = node_xyz[i+tet_node[j+tet*4]*3]; } } alpha = tetrahedron_barycentric ( tet_xyz, p ); if ( r8vec_is_nonnegative ( 4, alpha ) ) { tet_index = tet; *step_num = tet; return tet_index; } free ( alpha ); } return tet_index; } /******************************************************************************/ double *tetrahedron_barycentric ( double tetra[3*4], double p[3] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_BARYCENTRIC returns the barycentric coordinates of a point. Discussion: The barycentric coordinates of a point P with respect to a tetrahedron are a set of four values C(1:4), each associated with a vertex of the tetrahedron. The values must sum to 1. If all the values are between 0 and 1, the point is contained within the tetrahedron. The barycentric coordinate of point X related to vertex A can be interpreted as the ratio of the volume of the tetrahedron with vertex A replaced by vertex X to the volume of the original tetrahedron. Licensing: This code is distributed under the MIT license. Modified: 12 August 2005 Author: John Burkardt Parameters: Input, double TETRA[3*4], the vertices of the tetrahedron. Input, double P[3], the point to be checked. Output, double C[4], the barycentric coordinates of the point with respect to the tetrahedron. */ { # define N 3 # define RHS_NUM 1 double a[N*(N+RHS_NUM)]; double *c; int info; /* Set up the linear system ( X2-X1 X3-X1 X4-X1 ) C1 X - X1 ( Y2-Y1 Y3-Y1 Y4-Y1 ) C2 = Y - Y1 ( Z2-Z1 Z3-Z1 Z4-Z1 ) C3 Z - Z1 which is satisfied by the barycentric coordinates. */ a[0+0*N] = tetra[0+1*3] - tetra[0+0*3]; a[1+0*N] = tetra[1+1*3] - tetra[1+0*3]; a[2+0*N] = tetra[2+1*3] - tetra[2+0*3]; a[0+1*N] = tetra[0+2*3] - tetra[0+0*3]; a[1+1*N] = tetra[1+2*3] - tetra[1+0*3]; a[2+1*N] = tetra[2+2*3] - tetra[2+0*3]; a[0+2*N] = tetra[0+3*3] - tetra[0+0*3]; a[1+2*N] = tetra[1+3*3] - tetra[1+0*3]; a[2+2*N] = tetra[2+3*3] - tetra[2+0*3]; a[0+3*N] = p[0] - tetra[0+0*3]; a[1+3*N] = p[1] - tetra[1+0*3]; a[2+3*N] = p[2] - tetra[2+0*3]; /* Solve the linear system. */ info = r8mat_solve ( N, RHS_NUM, a ); if ( info != 0 ) { printf ( "\n" ); printf ( "TETRAHEDRON_BARYCENTRIC - Fatal error!\n" ); printf ( " The linear system is singular.\n" ); printf ( " The input data does not form a proper tetrahedron.\n" ); exit ( 1 ); } c = ( double * ) malloc ( 4 * sizeof ( double ) ); c[1] = a[0+3*N]; c[2] = a[1+3*N]; c[3] = a[2+3*N]; c[0] = 1.0 - c[1] - c[2] - c[3]; return c; # undef N # undef RHS_NUM } /******************************************************************************/ void tetrahedron_circumsphere_3d ( double tetra[3*4], double *r, double pc[3] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_CIRCUMSPHERE_3D computes the circumsphere of a tetrahedron in 3D. Discussion: The circumsphere, or circumscribed sphere, of a tetrahedron is the sphere that passes through the four vertices. The circumsphere is not necessarily the smallest sphere that contains the tetrahedron. Surprisingly, the diameter of the sphere can be found by solving a 3 by 3 linear system. This is because the vectors P2 - P1, P3 - P1 and P4 - P1 are secants of the sphere, and each forms a right triangle with the diameter through P1. Hence, the dot product of P2 - P1 with that diameter is equal to the square of the length of P2 - P1, and similarly for P3 - P1 and P4 - P1. This determines the diameter vector originating at P1, and hence the radius and center. Licensing: This code is distributed under the MIT license. Modified: 10 August 2005 Author: John Burkardt Reference: Adrian Bowyer, John Woodwark, A Programmer's Geometry, Butterworths, 1983. Parameters: Input, double TETRA[3*4], the coordinates of the vertices. Output, double *R, PC[3], the coordinates of the center of the circumscribed sphere, and its radius. If the linear system is singular, then R = -1, PC[] = 0. */ { # define DIM_NUM 3 # define RHS_NUM 1 double a[DIM_NUM*(DIM_NUM+RHS_NUM)]; int info; /* Set up the linear system. */ a[0+0*3] = tetra[0+1*3] - tetra[0+0*3]; a[0+1*3] = tetra[1+1*3] - tetra[1+0*3]; a[0+2*3] = tetra[2+1*3] - tetra[2+0*3]; a[0+3*3] = pow ( tetra[0+1*3] - tetra[0+0*3], 2 ) + pow ( tetra[1+1*3] - tetra[1+0*3], 2 ) + pow ( tetra[2+1*3] - tetra[2+0*3], 2 ); a[1+0*3] = tetra[0+2*3] - tetra[0+0*3]; a[1+1*3] = tetra[1+2*3] - tetra[1+0*3]; a[1+2*3] = tetra[2+2*3] - tetra[2+0*3]; a[1+3*3] = pow ( tetra[0+2*3] - tetra[0+0*3], 2 ) + pow ( tetra[1+2*3] - tetra[1+0*3], 2 ) + pow ( tetra[2+2*3] - tetra[2+0*3], 2 ); a[2+0*3] = tetra[0+3*3] - tetra[0+0*3]; a[2+1*3] = tetra[1+3*3] - tetra[1+0*3]; a[2+2*3] = tetra[2+3*3] - tetra[2+0*3]; a[2+3*3] = pow ( tetra[0+3*3] - tetra[0+0*3], 2 ) + pow ( tetra[1+3*3] - tetra[1+0*3], 2 ) + pow ( tetra[2+3*3] - tetra[2+0*3], 2 ); /* Solve the linear system. */ info = r8mat_solve ( DIM_NUM, RHS_NUM, a ); /* If the system was singular, return a consolation prize. */ if ( info != 0 ) { *r = -1.0; r8vec_zero ( DIM_NUM, pc ); return; } /* Compute the radius and center. */ *r = 0.5 * sqrt ( a[0+3*3] * a[0+3*3] + a[1+3*3] * a[1+3*3] + a[2+3*3] * a[2+3*3] ); pc[0] = tetra[0+0*3] + 0.5 * a[0+3*3]; pc[1] = tetra[1+0*3] + 0.5 * a[1+3*3]; pc[2] = tetra[2+0*3] + 0.5 * a[2+3*3]; return; # undef DIM_NUM # undef RHS_NUM } /******************************************************************************/ double *tetrahedron_edge_length_3d ( double tetra[3*4] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_EDGE_LENGTH_3D returns edge lengths of a tetrahedron in 3D. Licensing: This code is distributed under the MIT license. Modified: 10 August 2005 Author: John Burkardt Parameters: Input, double TETRA[3*4], the coordinates of the vertices. Output, double EDGE_LENGTH[6], the length of the edges. */ { # define DIM_NUM 3 double *edge_length; int i; int j1; int j2; int k; double v[DIM_NUM]; edge_length = ( double * ) malloc ( 6 * sizeof ( double ) ); k = 0; for ( j1 = 0; j1 < 3; j1++ ) { for ( j2 = j1 + 1; j2 < 4; j2++ ) { for ( i = 0; i < DIM_NUM; i++ ) { v[i] = tetra[i+j2*DIM_NUM] - tetra[i+j1*DIM_NUM]; } edge_length[k] = r8vec_length ( DIM_NUM, v ); k = k + 1; } } return edge_length; # undef DIM_NUM } /******************************************************************************/ void tetrahedron_insphere_3d ( double tetra[3*4], double *r, double pc[3] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_INSPHERE_3D finds the insphere of a tetrahedron in 3D. Discussion: The insphere of a tetrahedron is the inscribed sphere, which touches each face of the tetrahedron at a single point. The points of contact are the centroids of the triangular faces of the tetrahedron. Therefore, the point of contact for a face can be computed as the average of the vertices of that face. The sphere can then be determined as the unique sphere through the four given centroids. Licensing: This code is distributed under the MIT license. Modified: 08 August 2005 Author: John Burkardt Reference: Philip Schneider, David Eberly, Geometric Tools for Computer Graphics, Elsevier, 2002, ISBN: 1558605940, LC: T385.G6974. Parameters: Input, double TETRA[3*4], the coordinates of the vertices. Output, double *R, PC[3], the radius and the center of the sphere. */ { # define DIM_NUM 3 double b[4*4]; double gamma; int i; int j; double l123; double l124; double l134; double l234; double *n123; double *n124; double *n134; double *n234; double v21[DIM_NUM]; double v31[DIM_NUM]; double v41[DIM_NUM]; double v32[DIM_NUM]; double v42[DIM_NUM]; for ( i = 0; i < DIM_NUM; i++ ) { v21[i] = tetra[i+1*DIM_NUM] - tetra[i+0*DIM_NUM]; } for ( i = 0; i < DIM_NUM; i++ ) { v31[i] = tetra[i+2*DIM_NUM] - tetra[i+0*DIM_NUM]; } for ( i = 0; i < DIM_NUM; i++ ) { v41[i] = tetra[i+3*DIM_NUM] - tetra[i+0*DIM_NUM]; } for ( i = 0; i < DIM_NUM; i++ ) { v32[i] = tetra[i+2*DIM_NUM] - tetra[i+1*DIM_NUM]; } for ( i = 0; i < DIM_NUM; i++ ) { v42[i] = tetra[i+3*DIM_NUM] - tetra[i+1*DIM_NUM]; } n123 = r8vec_cross_3d ( v21, v31 ); n124 = r8vec_cross_3d ( v41, v21 ); n134 = r8vec_cross_3d ( v31, v41 ); n234 = r8vec_cross_3d ( v42, v32 ); l123 = r8vec_length ( DIM_NUM, n123 ); l124 = r8vec_length ( DIM_NUM, n124 ); l134 = r8vec_length ( DIM_NUM, n134 ); l234 = r8vec_length ( DIM_NUM, n234 ); free ( n123 ); free ( n124 ); free ( n134 ); free ( n234 ); for ( i = 0; i < DIM_NUM; i++ ) { pc[i] = ( l234 * tetra[i+0*DIM_NUM] + l134 * tetra[i+1*DIM_NUM] + l124 * tetra[i+2*DIM_NUM] + l123 * tetra[i+3*DIM_NUM] ) / ( l234 + l134 + l124 + l123 ); } for ( j = 0; j < 4; j++ ) { for ( i = 0; i < DIM_NUM; i++ ) { b[i+j*4] = tetra[i+j*DIM_NUM]; } b[3+j*4] = 1.0; } gamma = fabs ( r8mat_det_4d ( b ) ); *r = gamma / ( l234 + l134 + l124 + l123 ); return; # undef DIM_NUM } /******************************************************************************/ void tetrahedron_order4_physical_to_reference ( double tetra[], int n, double phy[], double ref[] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_ORDER4_PHYSICAL_TO_REFERENCE maps physical points to reference points. Discussion: Given the vertices of an order 4 physical tetrahedron and a point (X,Y,Z) in the physical tetrahedron, the routine computes the value of the corresponding image point (R,S,T) in reference space. This routine may be appropriate for an order 10 tetrahedron, if the mapping between reference and physical space is linear. This implies, in particular, that the edges of the image tetrahedron are straight, the faces are flat, and the "midside" nodes in the physical tetrahedron are halfway along the sides of the physical tetrahedron. Licensing: This code is distributed under the MIT license. Modified: 06 December 2006 Author: John Burkardt Parameters: Input, double TETRA[3*4], the coordinates of the vertices. The vertices are assumed to be the images of (0,0,0), (1,0,0), (0,1,0) and (0,0,1) respectively. Input, int N, the number of points to transform. Input, double PHY[3*N], the coordinates of physical points to be transformed. Output, double REF[3*N], the coordinates of the corresponding points in the reference space. */ { double a[3*3]; double det; int i; int j; /* Set up the matrix. */ for ( i = 0; i < 3; i++ ) { a[i+0*3] = tetra[i+1*3] - tetra[i+0*3]; a[i+1*3] = tetra[i+2*3] - tetra[i+0*3]; a[i+2*3] = tetra[i+3*3] - tetra[i+0*3]; } /* Compute the determinant. */ det = a[0+0*3] * ( a[1+1*3] * a[2+2*3] - a[1+2*3] * a[2+1*3] ) + a[0+1*3] * ( a[1+2*3] * a[2+0*3] - a[1+0*3] * a[2+2*3] ) + a[0+2*3] * ( a[1+0*3] * a[2+1*3] - a[1+1*3] * a[2+0*3] ); /* If the determinant is zero, bail out. */ if ( det == 0.0 ) { for ( j = 0; j < n; j++ ) { for ( i = 0; i < 3; i++ ) { ref[i+j*3] = 0.0; } } return; } /* Compute the solution. */ for ( j = 0; j < n; j++ ) { ref[0+j*3] = ( ( a[1+1*3] * a[2+2*3] - a[1+2*3] * a[2+1*3] ) * ( phy[0+j*3] - tetra[0+0*3] ) - ( a[0+1*3] * a[2+2*3] - a[0+2*3] * a[2+1*3] ) * ( phy[1+j*3] - tetra[1+0*3] ) + ( a[0+1*3] * a[1+2*3] - a[0+2*3] * a[1+1*3] ) * ( phy[2+j*3] - tetra[2+0*3] ) ) / det; ref[1+j*3] = ( - ( a[1+0*3] * a[2+2*3] - a[1+2*3] * a[2+0*3] ) * ( phy[0+j*3] - tetra[0+0*3] ) + ( a[0+0*3] * a[2+2*3] - a[0+2*3] * a[2+0*3] ) * ( phy[1+j*3] - tetra[1+0*3] ) - ( a[0+0*3] * a[1+2*3] - a[0+2*3] * a[1+0*3] ) * ( phy[2+j*3] - tetra[2+0*3] ) ) / det; ref[2+j*3] = ( ( a[1+0*3] * a[2+1*3] - a[1+1*3] * a[2+0*3] ) * ( phy[0+j*3] - tetra[0+0*3] ) - ( a[0+0*3] * a[2+1*3] - a[0+1*3] * a[2+0*3] ) * ( phy[1+j*3] - tetra[1+0*3] ) + ( a[0+0*3] * a[1+1*3] - a[0+1*3] * a[1+0*3] ) * ( phy[2+j*3] - tetra[2+0*3] ) ) / det; } return; } /******************************************************************************/ void tetrahedron_order4_reference_to_physical ( double tetra[], int n, double ref[], double phy[] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_ORDER4_REFERENCE_TO_PHYSICAL maps reference points to physical points. Discussion: Given the vertices of an order 4 physical tetrahedron and a point (R,S,T) in the reference triangle, the routine computes the value of the corresponding image point (X,Y,Z) in physical space. This routine will also be correct for an order 10 tetrahedron, if the mapping between reference and physical space is linear. This implies, in particular, that the sides of the image tetrahedron are straight, the faces are flat, and the "midside" nodes in the physical tetrahedron are halfway along the edges of the physical tetrahedron. Licensing: This code is distributed under the MIT license. Modified: 06 December 2006 Author: John Burkardt Parameters: Input, double TETRA[3*4], the coordinates of the vertices. The vertices are assumed to be the images of (0,0,0), (1,0,0), (0,1,0) and (0,0,1) respectively. Input, int N, the number of points to transform. Input, double REF[3*N], points in the reference tetrahedron Output, double PHY[3*N], corresponding points in the physical tetrahedron. */ { int i; int j; for ( i = 0; i < 3; i++ ) { for ( j = 0; j < n; j++ ) { phy[i+j*3] = tetra[i+0*3] * ( 1.0 - ref[0+j*3] - ref[1+j*3] - ref[2+j*3] ) + tetra[i+1*3] * + ref[0+j*3] + tetra[i+2*3] * + ref[1+j*3] + tetra[i+3*3] * + ref[2+j*3]; } } return; } /******************************************************************************/ double tetrahedron_quality1_3d ( double tetra[3*4] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_QUALITY1_3D: "quality" of a tetrahedron in 3D. Discussion: The quality of a tetrahedron is 3.0 times the ratio of the radius of the inscribed sphere divided by that of the circumscribed sphere. An equilateral tetrahredron achieves the maximum possible quality of 1. Licensing: This code is distributed under the MIT license. Modified: 20 September 2005 Author: John Burkardt Parameters: Input, double TETRA[3*4], the coordinates of the vertices. Output, double TETRAHEDRON_QUALITY1_3D, the quality of the tetrahedron. */ { # define DIM_NUM 3 double pc[DIM_NUM]; double quality; double r_in; double r_out; tetrahedron_circumsphere_3d ( tetra, &r_out, pc ); tetrahedron_insphere_3d ( tetra, &r_in, pc ); quality = 3.0 * r_in / r_out; return quality; # undef DIM_NUM } /******************************************************************************/ double tetrahedron_quality2_3d ( double tetra[3*4] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_QUALITY2_3D: "quality" of a tetrahedron in 3D. Discussion: The quality measure #2 of a tetrahedron is: QUALITY2 = 2 * sqrt ( 6 ) * RIN / LMAX where RIN = radius of the inscribed sphere; LMAX = length of longest side of the tetrahedron. An equilateral tetrahredron achieves the maximum possible quality of 1. Licensing: This code is distributed under the MIT license. Modified: 16 August 2005 Author: John Burkardt Reference: Qiang Du, Desheng Wang, The Optimal Centroidal Voronoi Tesselations and the Gersho's Conjecture in the Three-Dimensional Space, Computers and Mathematics with Applications, Volume 49, 2005, pages 1355-1373. Parameters: Input, double TETRA[3*4], the coordinates of the vertices. Output, double TETRAHEDRON_QUALITY2_3D, the quality of the tetrahedron. */ { # define DIM_NUM 3 double *edge_length; double l_max; double pc[DIM_NUM]; double quality2; double r_in; edge_length = tetrahedron_edge_length_3d ( tetra ); l_max = r8vec_max ( 6, edge_length ); tetrahedron_insphere_3d ( tetra, &r_in, pc ); quality2 = 2.0 * sqrt ( 6.0 ) * r_in / l_max; free ( edge_length ); return quality2; # undef DIM_NUM } /******************************************************************************/ double tetrahedron_quality3_3d ( double tetra[3*4] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_QUALITY3_3D computes the mean ratio of a tetrahedron. Discussion: This routine computes QUALITY3, the eigenvalue or mean ratio of a tetrahedron. QUALITY3 = 12 * ( 3 * volume )**(2/3) / (sum of square of edge lengths). This value may be used as a shape quality measure for the tetrahedron. For an equilateral tetrahedron, the value of this quality measure will be 1. For any other tetrahedron, the value will be between 0 and 1. Licensing: This code is distributed under the MIT license. Modified: 17 August 2005 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 TETRA(3,4), the coordinates of the vertices. Output, double TETRAHEDRON_QUALITY3_3D, the mean ratio of the tetrahedron. */ { # define DIM_NUM 3 double ab[DIM_NUM]; double ac[DIM_NUM]; double ad[DIM_NUM]; double bc[DIM_NUM]; double bd[DIM_NUM]; double cd[DIM_NUM]; double denom; int i; double lab; double lac; double lad; double lbc; double lbd; double lcd; double quality3; double volume; /* Compute the vectors representing the sides of the tetrahedron. */ for ( i = 0; i < DIM_NUM; i++ ) { ab[i] = tetra[i+1*DIM_NUM] - tetra[i+0*DIM_NUM]; ac[i] = tetra[i+2*DIM_NUM] - tetra[i+0*DIM_NUM]; ad[i] = tetra[i+3*DIM_NUM] - tetra[i+0*DIM_NUM]; bc[i] = tetra[i+2*DIM_NUM] - tetra[i+1*DIM_NUM]; bd[i] = tetra[i+3*DIM_NUM] - tetra[i+1*DIM_NUM]; cd[i] = tetra[i+3*DIM_NUM] - tetra[i+2*DIM_NUM]; } /* Compute the squares of the lengths of the sides. */ lab = pow ( ab[0], 2 ) + pow ( ab[1], 2 ) + pow ( ab[2], 2 ); lac = pow ( ac[0], 2 ) + pow ( ac[1], 2 ) + pow ( ac[2], 2 ); lad = pow ( ad[0], 2 ) + pow ( ad[1], 2 ) + pow ( ad[2], 2 ); lbc = pow ( bc[0], 2 ) + pow ( bc[1], 2 ) + pow ( bc[2], 2 ); lbd = pow ( bd[0], 2 ) + pow ( bd[1], 2 ) + pow ( bd[2], 2 ); lcd = pow ( cd[0], 2 ) + pow ( cd[1], 2 ) + pow ( cd[2], 2 ); /* Compute the volume. */ volume = fabs ( ab[0] * ( ac[1] * ad[2] - ac[2] * ad[1] ) + ab[1] * ( ac[2] * ad[0] - ac[0] * ad[2] ) + ab[2] * ( ac[0] * ad[1] - ac[1] * ad[0] ) ) / 6.0; denom = lab + lac + lad + lbc + lbd + lcd; if ( denom == 0.0 ) { quality3 = 0.0; } else { quality3 = 12.0 * pow ( 3.0 * volume, 2.0 / 3.0 ) / denom; } return quality3; # undef DIM_NUM } /******************************************************************************/ double tetrahedron_quality4_3d ( double tetra[3*4] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_QUALITY4_3D computes the minimum solid angle of a tetrahedron. Discussion: This routine computes a quality measure for a tetrahedron, based on the sine of half the minimum of the four solid angles. Modified: 17 August 2005 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 TETRA[3*4], the coordinates of the vertices. Output, double QUALITY4, the value of the quality measure. */ { # define DIM_NUM 3 double ab[DIM_NUM]; double ac[DIM_NUM]; double ad[DIM_NUM]; double bc[DIM_NUM]; double bd[DIM_NUM]; double cd[DIM_NUM]; double denom; int i; double l1; double l2; double l3; double lab; double lac; double lad; double lbc; double lbd; double lcd; double quality4; double volume; /* Compute the vectors that represent the sides. */ for ( i = 0; i < DIM_NUM; i++ ) { ab[i] = tetra[i+1*DIM_NUM] - tetra[i+0*DIM_NUM]; ac[i] = tetra[i+2*DIM_NUM] - tetra[i+0*DIM_NUM]; ad[i] = tetra[i+3*DIM_NUM] - tetra[i+0*DIM_NUM]; bc[i] = tetra[i+2*DIM_NUM] - tetra[i+1*DIM_NUM]; bd[i] = tetra[i+3*DIM_NUM] - tetra[i+1*DIM_NUM]; cd[i] = tetra[i+3*DIM_NUM] - tetra[i+2*DIM_NUM]; } /* Compute the lengths of the sides. */ lab = r8vec_length ( DIM_NUM, ab ); lac = r8vec_length ( DIM_NUM, ac ); lad = r8vec_length ( DIM_NUM, ad ); lbc = r8vec_length ( DIM_NUM, bc ); lbd = r8vec_length ( DIM_NUM, bd ); lcd = r8vec_length ( DIM_NUM, cd ); /* Compute the volume. */ volume = fabs ( ab[0] * ( ac[1] * ad[2] - ac[2] * ad[1] ) + ab[1] * ( ac[2] * ad[0] - ac[0] * ad[2] ) + ab[2] * ( ac[0] * ad[1] - ac[1] * ad[0] ) ) / 6.0; quality4 = 1.0; l1 = lab + lac; l2 = lab + lad; l3 = lac + lad; denom = ( l1 + lbc ) * ( l1 - lbc ) * ( l2 + lbd ) * ( l2 - lbd ) * ( l3 + lcd ) * ( l3 - lcd ); if ( denom <= 0.0 ) { quality4 = 0.0; } else { quality4 = fmin ( quality4, 12.0 * volume / sqrt ( denom ) ); } l1 = lab + lbc; l2 = lab + lbd; l3 = lbc + lbd; denom = ( l1 + lac ) * ( l1 - lac ) * ( l2 + lad ) * ( l2 - lad ) * ( l3 + lcd ) * ( l3 - lcd ); if ( denom <= 0.0 ) { quality4 = 0.0; } else { quality4 = fmin ( quality4, 12.0 * volume / sqrt ( denom ) ); } l1 = lac + lbc; l2 = lac + lcd; l3 = lbc + lcd; denom = ( l1 + lab ) * ( l1 - lab ) * ( l2 + lad ) * ( l2 - lad ) * ( l3 + lbd ) * ( l3 - lbd ); if ( denom <= 0.0 ) { quality4 = 0.0; } else { quality4 = fmin ( quality4, 12.0 * volume / sqrt ( denom ) ); } l1 = lad + lbd; l2 = lad + lcd; l3 = lbd + lcd; denom = ( l1 + lab ) * ( l1 - lab ) * ( l2 + lac ) * ( l2 - lac ) * ( l3 + lbc ) * ( l3 - lbc ); if ( denom <= 0.0 ) { quality4 = 0.0; } else { quality4 = fmin ( quality4, 12.0 * volume / sqrt ( denom ) ); } quality4 = quality4 * 1.5 * sqrt ( 6.0 ); return quality4; # undef DIM_NUM } /******************************************************************************/ void tetrahedron_reference_sample ( int n, int *seed, double p[] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_REFERENCE_SAMPLE samples points in the reference tetrahedron. Licensing: This code is distributed under the MIT license. Modified: 08 December 2006 Author: John Burkardt Parameters: Input, int N, the number of points to sample. Input/output, int *SEED, a seed for the random number generator. Output, double P[3*N], random points in the tetrahedron. */ { double alpha; double beta; double gamma; int j; double r; for ( j = 0; j < n; j++ ) { r = r8_uniform_01 ( seed ); /* Interpret R as a percentage of the tetrahedron's volume. Imagine a plane, parallel to face 1, so that the volume between vertex 1 and the plane is R percent of the full tetrahedron volume. The plane will intersect sides 12, 13, and 14 at a fraction ALPHA = R^1/3 of the distance from vertex 1 to vertices 2, 3, and 4. */ alpha = pow ( r, 1.0 / 3.0 ); /* Determine the coordinates of the points on sides 12, 13 and 14 intersected by the plane, which form a triangle TR. Now choose, uniformly at random, a point in this triangle. */ r = r8_uniform_01 ( seed ); /* Interpret R as a percentage of the triangle's area. Imagine a line L, parallel to side 1, so that the area between vertex 1 and line L is R percent of the full triangle's area. The line L will intersect sides 2 and 3 at a fraction ALPHA = SQRT ( R ) of the distance from vertex 1 to vertices 2 and 3. */ beta = sqrt ( r ); /* Determine the coordinates of the points on sides 2 and 3 intersected by line L. Now choose, uniformly at random, a point on the line L. */ gamma = r8_uniform_01 ( seed ); p[0+j*3] = alpha * ( 1.0 - beta ) * gamma; p[1+j*3] = alpha * beta * ( 1.0 - gamma ); p[2+j*3] = alpha * beta * gamma; } return; } /******************************************************************************/ void tetrahedron_sample ( double tetra[3*4], int n, int *seed, double p[] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_SAMPLE returns random points in a tetrahedron. Licensing: This code is distributed under the MIT license. Modified: 06 December 2006 Author: John Burkardt Parameters: Input, double TETRA[3*4], the coordinates of the vertices. Input/output, int *SEED, a seed for the random number generator. Output, double P[3*N], random points in the tetrahedron. */ { # define DIM_NUM 3 double alpha; double beta; double gamma; int i; int j; int k; double *p12; double *p13; double r; double *t; p12 = ( double * ) malloc ( DIM_NUM * sizeof ( double ) ); p13 = ( double * ) malloc ( DIM_NUM * sizeof ( double ) ); t = ( double * ) malloc ( DIM_NUM * 3 * sizeof ( double ) ); for ( k = 0; k < n; k++ ) { r = r8_uniform_01 ( seed ); /* Interpret R as a percentage of the tetrahedron's volume. Imagine a plane, parallel to face 1, so that the volume between vertex 1 and the plane is R percent of the full tetrahedron volume. The plane will intersect sides 12, 13, and 14 at a fraction ALPHA = R^1/3 of the distance from vertex 1 to vertices 2, 3, and 4. */ alpha = pow ( r, 1.0 / 3.0 ); /* Determine the coordinates of the points on sides 12, 13 and 14 intersected by the plane, which form a triangle TR. */ for ( i = 0; i < DIM_NUM; i++ ) { for ( j = 0; j < 3; j++ ) { t[i+j*3] = ( 1.0 - alpha ) * tetra[i+0*3] + alpha * tetra[i+(j+1)*3]; } } /* Now choose, uniformly at random, a point in this triangle. */ r = r8_uniform_01 ( seed ); /* Interpret R as a percentage of the triangle's area. Imagine a line L, parallel to side 1, so that the area between vertex 1 and line L is R percent of the full triangle's area. The line L will intersect sides 2 and 3 at a fraction ALPHA = SQRT ( R ) of the distance from vertex 1 to vertices 2 and 3. */ beta = sqrt ( r ); /* Determine the coordinates of the points on sides 2 and 3 intersected by line L. */ for ( i = 0; i < DIM_NUM; i++ ) { p12[i] = ( 1.0 - beta ) * t[i+0*3] + beta * t[i+1*3]; p13[i] = ( 1.0 - beta ) * t[i+0*3] + beta * t[i+2*3]; } /* Now choose, uniformly at random, a point on the line L. */ gamma = r8_uniform_01 ( seed ); for ( i = 0; i < DIM_NUM; i++ ) { p[i+k*3] = gamma * p12[i] + ( 1.0 - gamma ) * p13[i]; } } free ( p12 ); free ( p13 ); free ( t ); return; # undef DIM_NUM } /******************************************************************************/ double tetrahedron_volume ( double tetra[3*4] ) /******************************************************************************/ /* Purpose: TETRAHEDRON_VOLUME computes the volume of a tetrahedron in 3D. Licensing: This code is distributed under the MIT license. Modified: 06 August 2005 Author: John Burkardt Parameters: Input, double TETRA[3*4], the coordinates of the vertices. Output, double TETRAHEDRON_VOLUME, the volume of the tetrahedron. */ { double a[4*4]; int i; int j; double volume; for ( i = 0; i < 3; i++ ) { for ( j = 0; j < 4; j++ ) { a[i+j*4] = tetra[i+j*3]; } } i = 3; for ( j = 0; j < 4; j++ ) { a[i+j*4] = 1.0; } volume = fabs ( r8mat_det_4d ( a ) ) / 6.0; return volume; } /******************************************************************************/ void timestamp ( ) /******************************************************************************/ /* Purpose: TIMESTAMP prints the current YMDHMS date as a time stamp. Example: 31 May 2001 09:45:54 AM Licensing: This code is distributed under the MIT license. Modified: 24 September 2003 Author: John Burkardt Parameters: None */ { # define TIME_SIZE 40 static char time_buffer[TIME_SIZE]; const struct tm *tm; time_t now; now = time ( NULL ); tm = localtime ( &now ); strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm ); fprintf ( stdout, "%s\n", time_buffer ); return; # undef TIME_SIZE }