# include # include # include # include # include int main ( int argc, char *argv[] ); char ch_cap ( char ch ); int ch_eqi ( char ch1, char ch2 ); int ch_to_digit ( char ch ); int file_column_count ( char *filename ); int file_row_count ( char *filename ); int i4_max ( int i1, int i2 ); int i4_min ( int i1, int i2 ); int i4col_compare ( int m, int n, int a[], int i, int j ); void i4col_sort_a ( int m, int n, int a[] ); void i4col_swap ( int m, int n, int a[], int irow1, int irow2 ); int *i4mat_data_read ( char *input_filename, int m, int n ); void i4mat_header_read ( char *input_filename, int *m, int *n ); void i4mat_transpose_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi, int jhi, char *title ); void i4mat_write ( char *output_filename, int m, int n, int table[] ); void mesh_base_zero ( int node_num, int element_order, int element_num, int element_node[] ); double *r8mat_data_read ( char *input_filename, int m, int n ); void r8mat_header_read ( char *input_filename, int *m, int *n ); void r8mat_transpose_print_some ( int m, int n, double a[], int ilo, int jlo, int ihi, int jhi, char *title ); int s_len_trim ( char *s ); int s_to_i4 ( char *s, int *last, int *error ); int s_to_i4vec ( char *s, int n, int ivec[] ); double s_to_r8 ( char *s, int *lchar, int *error ); int s_to_r8vec ( char *s, int n, double rvec[] ); int s_word_count ( char *s ); void sort_heap_external ( int n, int *indx, int *i, int *j, int isgn ); void timestamp ( ); int *triangulation_neighbor_triangles ( int triangle_order, int triangle_num, int triangle_node[] ); /******************************************************************************/ int main ( int argc, char *argv[] ) /******************************************************************************/ /* Purpose: MAIN is the main program for TRIANGULATION_TRIANGLE_NEIGHBORS. Discussion: TRIANGULATION_TRIANGLE_NEIGHBORS determines the neighbor triangles of each triangle in a triangulation. The user supplies a node file and a triangle file, containing the coordinates of the nodes, and the indices of the nodes that make up each triangle. Either 3-node or 6-node triangles may be used. The program reads the node and triangle data, computes the triangle neighbor information, and writes it to a file. Usage: triangulation_triangle_neighbors prefix where 'prefix' is the common filename prefix: * prefix_nodes.txt contains the node coordinates, * prefix_elements.txt contains the element definitions. * prefix_element_neighbors.txt will contain the triangle neighbors. Licensing: This code is distributed under the MIT license. Modified: 08 October 2010 Author: John Burkardt */ { int dim_num; char neighbor_filename[255]; char node_filename[255]; int node_num; double *node_xy; char prefix[255]; char element_filename[255]; int *triangle_neighbor; int *triangle_node; int triangle_num; int triangle_order; printf ( "\n" ); timestamp ( ); printf ( "\n" ); printf ( "TRIANGULATION_TRIANGLE_NEIGHBORS.\n" ); printf ( " C version:\n" ); printf ( " Read a node dataset of NODE_NUM points in 2 dimensions.\n" ); printf ( " Read an associated triangulation dataset of \n" ); printf ( " TRIANGLE_NUM triangles using 3 or 6 nodes.\n" ); printf ( "\n" ); printf ( " For each triangle, determine the indices of the\n" ); printf ( " triangles opposite vertices 1, 2 and 3.\n" ); printf ( "\n" ); printf ( " Write this triangle neighbor data to files.\n" ); printf ( "\n" ); printf ( " Compiled on %s at %s\n", __DATE__, __TIME__ ); /* Get the filename prefix. */ if ( argc <= 1 ) { printf ( "\n" ); printf ( "TRIANGULATION_TRIANGLE_NEIGHBORS:\n" ); printf ( " Please enter the filename prefix.\n" ); scanf ( "%s", prefix ); } else { strcpy ( prefix, argv[1] ); } /* Create the filenames. */ strcpy ( node_filename, prefix ); strcat ( node_filename, "_nodes.txt" ); strcpy ( element_filename, prefix ); strcat ( element_filename, "_elements.txt" ); strcpy ( neighbor_filename, prefix ); strcat ( neighbor_filename, "_neighbors.txt" ); /* Read the node data. */ r8mat_header_read ( node_filename, &dim_num, &node_num ); printf ( "\n" ); printf ( " Read the header of \"%s\".\n", node_filename ); printf ( "\n" ); printf ( " Spatial dimension DIM_NUM = %d\n", dim_num ); printf ( " Number of nodes NODE_NUM = %d\n", node_num ); node_xy = r8mat_data_read ( node_filename, dim_num, node_num ); printf ( "\n" ); printf ( " Read the data in \"%s\".\n", node_filename ); r8mat_transpose_print_some ( dim_num, node_num, node_xy, 1, 1, dim_num, 5, " Portion of coordinate data from file:" ); /* Read the element data. */ i4mat_header_read ( element_filename, &triangle_order, &triangle_num ); printf ( "\n" ); printf ( " Read the header of \"%s\".\n", element_filename ); printf ( "\n" ); printf ( " Triangle order TRIANGLE_ORDER = %d\n", triangle_order ); printf ( " Number of triangles TRIANGLE_NUM = %d\n", triangle_num ); triangle_node = i4mat_data_read ( element_filename, triangle_order, triangle_num ); printf ( "\n" ); printf ( " Read the data in \"%s\".\n", element_filename ); i4mat_transpose_print_some ( triangle_order, triangle_num, triangle_node, 1, 1, triangle_order, 5, " Portion of data read from file:" ); /* Detect and correct 1-based node indexing. */ mesh_base_zero ( node_num, triangle_order, triangle_num, triangle_node ); /* Create the triangle neighbors. */ triangle_neighbor = triangulation_neighbor_triangles ( triangle_order, triangle_num, triangle_node ); /* Write the output file. */ i4mat_write ( neighbor_filename, 3, triangle_num, triangle_neighbor ); printf ( "\n" ); printf ( " Created the triangle neighbor file \"%s\"\n", neighbor_filename ); free ( node_xy ); free ( triangle_neighbor ); free ( triangle_node ); /* Terminate. */ printf ( "\n" ); printf ( "TRIANGULATION_TRIANGLE_NEIGHBORS.:\n" ); printf ( " Normal end of execution.\n" ); printf ( "\n" ); timestamp ( ); return 0; } /******************************************************************************/ char ch_cap ( char ch ) /******************************************************************************/ /* Purpose: CH_CAP capitalizes a single character. Discussion: This routine should be equivalent to the library "toupper" function. Licensing: This code is distributed under the MIT license. Modified: 19 July 1998 Author: John Burkardt Parameters: Input, char CH, the character to capitalize. Output, char CH_CAP, the capitalized character. */ { if ( 97 <= ch && ch <= 122 ) { ch = ch - 32; } return ch; } /******************************************************************************/ int ch_eqi ( char ch1, char ch2 ) /******************************************************************************/ /* Purpose: CH_EQI is TRUE (1) if two characters are equal, disregarding case. Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char CH1, CH2, the characters to compare. Output, int CH_EQI, is TRUE (1) if the two characters are equal, disregarding case and FALSE (0) otherwise. */ { int value; if ( 97 <= ch1 && ch1 <= 122 ) { ch1 = ch1 - 32; } if ( 97 <= ch2 && ch2 <= 122 ) { ch2 = ch2 - 32; } if ( ch1 == ch2 ) { value = 1; } else { value = 0; } return value; } /******************************************************************************/ int ch_to_digit ( char ch ) /******************************************************************************/ /* Purpose: CH_TO_DIGIT returns the integer value of a base 10 digit. Example: CH DIGIT --- ----- '0' 0 '1' 1 ... ... '9' 9 ' ' 0 'X' -1 Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char CH, the decimal digit, '0' through '9' or blank are legal. Output, int CH_TO_DIGIT, the corresponding integer value. If the character was 'illegal', then DIGIT is -1. */ { int digit; if ( '0' <= ch && ch <= '9' ) { digit = ch - '0'; } else if ( ch == ' ' ) { digit = 0; } else { digit = -1; } return digit; } /******************************************************************************/ int file_column_count ( char *input_filename ) /******************************************************************************/ /* Purpose: FILE_COLUMN_COUNT counts the number of columns in the first line of a file. Discussion: The file is assumed to be a simple text file. Most lines of the file is presumed to consist of COLUMN_NUM words, separated by spaces. There may also be some blank lines, and some comment lines, which have a "#" in column 1. The routine tries to find the first non-comment non-blank line and counts the number of words in that line. If all lines are blanks or comments, it goes back and tries to analyze a comment line. Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the file. Output, int FILE_COLUMN_COUNT, the number of columns assumed to be in the file. */ { # define LINE_MAX 255 int column_num; char *error; FILE *input; int got_one; char line[LINE_MAX]; /* Open the file. */ input = fopen ( input_filename, "r" ); if ( !input ) { column_num = -1; printf ( "\n" ); printf ( "FILE_COLUMN_COUNT - Fatal error!\n" ); printf ( " Could not open the input file: \"%s\"\n", input_filename ); return column_num; } /* Read one line, but skip blank lines and comment lines. */ got_one = 0; for ( ; ; ) { error = fgets ( line, LINE_MAX, input ); if ( !error ) { break; } if ( s_len_trim ( line ) == 0 ) { continue; } if ( line[0] == '#' ) { continue; } got_one = 1; break; } if ( got_one == 0 ) { fclose ( input ); input = fopen ( input_filename, "r" ); for ( ; ; ) { error = fgets ( line, LINE_MAX, input ); if ( !error ) { break; } if ( s_len_trim ( line ) == 0 ) { continue; } got_one = 1; break; } } fclose ( input ); if ( got_one == 0 ) { printf ( "\n" ); printf ( "FILE_COLUMN_COUNT - Warning!\n" ); printf ( " The file does not seem to contain any data.\n" ); return -1; } column_num = s_word_count ( line ); return column_num; # undef LINE_MAX } /******************************************************************************/ int file_row_count ( char *input_filename ) /******************************************************************************/ /* Purpose: FILE_ROW_COUNT counts the number of row records in a file. Discussion: It does not count lines that are blank, or that begin with a comment symbol '#'. Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Output, int FILE_ROW_COUNT, the number of rows found. */ { # define LINE_MAX 255 int comment_num; char *error; FILE *input; char line[LINE_MAX]; int record_num; int row_num; row_num = 0; comment_num = 0; record_num = 0; input = fopen ( input_filename, "r" ); if ( !input ) { printf ( "\n" ); printf ( "FILE_ROW_COUNT - Fatal error!\n" ); printf ( " Could not open the input file: \"%s\"\n", input_filename ); return (-1); } for ( ; ; ) { error = fgets ( line, LINE_MAX, input ); if ( !error ) { break; } record_num = record_num + 1; if ( line[0] == '#' ) { comment_num = comment_num + 1; continue; } if ( s_len_trim ( line ) == 0 ) { comment_num = comment_num + 1; continue; } row_num = row_num + 1; } fclose ( input ); return row_num; # undef LINE_MAX } /******************************************************************************/ 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: 29 August 2006 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 smaller of two I4's. Licensing: This code is distributed under the MIT license. Modified: 29 August 2006 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; } /******************************************************************************/ 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_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 } /******************************************************************************/ int *i4mat_data_read ( char *input_filename, int m, int n ) /******************************************************************************/ /* Purpose: I4MAT_DATA_READ reads the data from an I4MAT file. Discussion: An I4MAT is an array of I4's. The file is assumed to contain one record per line. Records beginning with the '#' character are comments, and are ignored. Blank lines are also ignored. Each line that is not ignored is assumed to contain exactly (or at least) M real numbers, representing the coordinates of a point. There are assumed to be exactly (or at least) N such records. Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Input, int M, the number of spatial dimensions. Input, int N, the number of points. The program will stop reading data once N values have been read. Output, int I4MAT_DATA_READ[M*N], the table data. */ { # define LINE_MAX 255 int error; char *got_string; FILE *input; int i; int j; char line[255]; int *table; int *x; input = fopen ( input_filename, "r" ); if ( !input ) { printf ( "\n" ); printf ( "I4MAT_DATA_READ - Fatal error!\n" ); printf ( " Could not open the input file: \"%s\"\n", input_filename ); return NULL; } table = ( int * ) malloc ( m * n * sizeof ( int ) ); x = ( int * ) malloc ( m * sizeof ( int ) ); j = 0; while ( j < n ) { got_string = fgets ( line, LINE_MAX, input ); if ( !got_string ) { break; } if ( line[0] == '#' || s_len_trim ( line ) == 0 ) { continue; } error = s_to_i4vec ( line, m, x ); if ( error == 1 ) { continue; } for ( i = 0; i < m; i++ ) { table[i+j*m] = x[i]; } j = j + 1; } fclose ( input ); free ( x ); return table; # undef LINE_MAX } /******************************************************************************/ void i4mat_header_read ( char *input_filename, int *m, int *n ) /******************************************************************************/ /* Purpose: I4MAT_HEADER_READ reads the header from an I4MAT file. Discussion: An I4MAT is an array of I4's. Licensing: This code is distributed under the MIT license. Modified: 28 May 2008 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Output, int *M, the number of spatial dimensions. Output, int *N, the number of points. */ { *m = file_column_count ( input_filename ); if ( *m <= 0 ) { printf ( "\n" ); printf ( "I4MAT_HEADER_READ - Fatal error!\n" ); printf ( " FILE_COLUMN_COUNT failed.\n" ); *n = -1; return; } *n = file_row_count ( input_filename ); if ( *n <= 0 ) { printf ( "\n" ); printf ( "I4MAT_HEADER_READ - Fatal error!\n" ); printf ( " FILE_ROW_COUNT failed.\n" ); return; } 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: 14 June 2005 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 ); /* Print the columns of the matrix, in strips of INCX. */ for ( i2lo = ilo; i2lo <= ihi; i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( i2hi, ihi ); fprintf ( stdout, "\n" ); /* For each row I in the current range... Write the header. */ fprintf ( stdout, " Row: " ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, "%6d ", i ); } fprintf ( stdout, "\n" ); fprintf ( stdout, " Col\n" ); fprintf ( stdout, "\n" ); /* Determine the range of the rows in this strip. */ j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { /* Print out (up to INCX) entries in column J, that lie in the current strip. */ fprintf ( stdout, "%5d: ", j ); for ( i = i2lo; i <= i2hi; i++ ) { fprintf ( stdout, "%6d ", a[i-1+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ void i4mat_write ( char *output_filename, int m, int n, int table[] ) /******************************************************************************/ /* Purpose: I4MAT_WRITE writes an I4MAT file. Discussion: An I4MAT is an array of I4's. Licensing: This code is distributed under the MIT license. Modified: 01 June 2009 Author: John Burkardt Parameters: Input, char *OUTPUT_FILENAME, the output filename. Input, int M, the spatial dimension. Input, int N, the number of points. Input, int TABLE[M*N], the table data. */ { int i; int j; FILE *output; /* Open the file. */ output = fopen ( output_filename, "wt" ); if ( !output ) { printf ( "\n" ); printf ( "I4MAT_WRITE - Fatal error!\n" ); printf ( " Could not open the output file.\n" ); return; } /* Write the data. */ for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { fprintf ( output, " %d", table[i+j*m] ); } fprintf ( output, "\n" ); } /* Close the file. */ fclose ( output ); 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 or 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: 08 October 2010 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; int node; int node_max; int node_min; int order; node_min = node_num + 1; node_max = -1; for ( element = 0; element < element_num; element++ ) { for ( order = 0; order < element_order; order++ ) { node = element_node[order+element*element_order]; node_min = i4_min ( node_min, node ); node_max = i4_max ( 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; } /******************************************************************************/ double *r8mat_data_read ( char *input_filename, int m, int n ) /******************************************************************************/ /* Purpose: R8MAT_DATA_READ reads the data from an R8MAT file. Discussion: An R8MAT is an array of R8's. The file is assumed to contain one record per line. Records beginning with the '#' character are comments, and are ignored. Blank lines are also ignored. Each line that is not ignored is assumed to contain exactly (or at least) M real numbers, representing the coordinates of a point. There are assumed to be exactly (or at least) N such records. Licensing: This code is distributed under the MIT license. Modified: 27 January 2005 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Input, int M, the number of spatial dimensions. Input, int N, the number of points. The program will stop reading data once N values have been read. Output, double R8MAT_DATA_READ[M*N], the table data. */ { # define LINE_MAX 255 int error; char *got_string; FILE *input; int i; int j; char line[255]; double *table; double *x; input = fopen ( input_filename, "r" ); if ( !input ) { printf ( "\n" ); printf ( "R8MAT_DATA_READ - Fatal error!\n" ); printf ( " Could not open the input file: \"%s\"\n", input_filename ); return NULL; } table = ( double * ) malloc ( m * n * sizeof ( double ) ); x = ( double * ) malloc ( m * sizeof ( double ) ); j = 0; while ( j < n ) { got_string = fgets ( line, LINE_MAX, input ); if ( !got_string ) { break; } if ( line[0] == '#' || s_len_trim ( line ) == 0 ) { continue; } error = s_to_r8vec ( line, m, x ); if ( error == 1 ) { continue; } for ( i = 0; i < m; i++ ) { table[i+j*m] = x[i]; } j = j + 1; } fclose ( input ); free ( x ); return table; # undef LINE_MAX } /******************************************************************************/ void r8mat_header_read ( char *input_filename, int *m, int *n ) /******************************************************************************/ /* Purpose: R8MAT_HEADER_READ reads the header from an R8MAT file. Discussion: An R8MAT is an array of R8's. Licensing: This code is distributed under the MIT license. Modified: 04 June 2004 Author: John Burkardt Parameters: Input, char *INPUT_FILENAME, the name of the input file. Output, int *M, the number of spatial dimensions. Output, int *N, the number of points. */ { *m = file_column_count ( input_filename ); if ( *m <= 0 ) { printf ( "\n" ); printf ( "R8MAT_HEADER_READ - Fatal error!\n" ); printf ( " FILE_COLUMN_COUNT failed.\n" ); *n = -1; return; } *n = file_row_count ( input_filename ); if ( *n <= 0 ) { printf ( "\n" ); printf ( "R8MAT_HEADER_READ - Fatal error!\n" ); printf ( " FILE_ROW_COUNT failed.\n" ); return; } 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. Licensing: This code is distributed under the MIT license. Modified: 20 August 2010 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 inc; int j; int j2hi; int j2lo; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); for ( i2lo = i4_max ( ilo, 1 ); i2lo <= i4_min ( ihi, m ); i2lo = i2lo + INCX ) { i2hi = i2lo + INCX - 1; i2hi = i4_min ( i2hi, m ); i2hi = i4_min ( 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" ); j2lo = i4_max ( jlo, 1 ); j2hi = i4_min ( jhi, n ); for ( j = j2lo; j <= j2hi; j++ ) { fprintf ( stdout, "%5d:", j - 1 ); for ( i2 = 1; i2 <= inc; i2++ ) { i = i2lo - 1 + i2; fprintf ( stdout, " %14f", a[(i-1)+(j-1)*m] ); } fprintf ( stdout, "\n" ); } } return; # undef INCX } /******************************************************************************/ int s_len_trim ( char *s ) /******************************************************************************/ /* Purpose: S_LEN_TRIM returns the length of a string to the last nonblank. Licensing: This code is distributed under the MIT license. Modified: 26 April 2003 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 != ' ' ) { return n; } t--; n--; } return n; } /******************************************************************************/ int s_to_i4 ( char *s, int *last, int *error ) /******************************************************************************/ /* Purpose: S_TO_I4 reads an I4 from a string. Licensing: This code is distributed under the MIT license. Modified: 13 June 2003 Author: John Burkardt Parameters: Input, char *S, a string to be examined. Output, int *LAST, the last character of S used to make IVAL. Output, int *ERROR is TRUE (1) if an error occurred and FALSE (0) otherwise. Output, int *S_TO_I4, the integer value read from the string. If the string is blank, then IVAL will be returned 0. */ { char c; int i; int isgn; int istate; int ival; *error = 0; istate = 0; isgn = 1; i = 0; ival = 0; while ( *s ) { c = s[i]; i = i + 1; /* Haven't read anything. */ if ( istate == 0 ) { if ( c == ' ' ) { } else if ( c == '-' ) { istate = 1; isgn = -1; } else if ( c == '+' ) { istate = 1; isgn = + 1; } else if ( '0' <= c && c <= '9' ) { istate = 2; ival = c - '0'; } else { *error = 1; return ival; } } /* Have read the sign, expecting digits. */ else if ( istate == 1 ) { if ( c == ' ' ) { } else if ( '0' <= c && c <= '9' ) { istate = 2; ival = c - '0'; } else { *error = 1; return ival; } } /* Have read at least one digit, expecting more. */ else if ( istate == 2 ) { if ( '0' <= c && c <= '9' ) { ival = 10 * (ival) + c - '0'; } else { ival = isgn * ival; *last = i - 1; return ival; } } } /* If we read all the characters in the string, see if we're OK. */ if ( istate == 2 ) { ival = isgn * ival; *last = s_len_trim ( s ); } else { *error = 1; *last = 0; } return ival; } /******************************************************************************/ int s_to_i4vec ( char *s, int n, int ivec[] ) /******************************************************************************/ /* Purpose: S_TO_I4VEC reads an I4VEC from a string. Licensing: This code is distributed under the MIT license. Modified: 19 February 2001 Author: John Burkardt Parameters: Input, char *S, the string to be read. Input, int N, the number of values expected. Output, int IVEC[N], the values read from the string. Output, int S_TO_I4VEC, is TRUE (1) if an error occurred and FALSE (0) otherwise. */ { int error; int i; int lchar; error = 0; for ( i = 0; i < n; i++ ) { ivec[i] = s_to_i4 ( s, &lchar, &error ); if ( error ) { return error; } s = s + lchar; } return error; } /******************************************************************************/ double s_to_r8 ( char *s, int *lchar, int *error ) /******************************************************************************/ /* Purpose: S_TO_R8 reads an R8 value from a string. Discussion: We have had some trouble with input of the form 1.0E-312. For now, let's assume anything less than 1.0E-20 is zero. This routine will read as many characters as possible until it reaches the end of the string, or encounters a character which cannot be part of the real number. Legal input is: 1 blanks, 2 '+' or '-' sign, 2.5 spaces 3 integer part, 4 decimal point, 5 fraction part, 6 'E' or 'e' or 'D' or 'd', exponent marker, 7 exponent sign, 8 exponent integer part, 9 exponent decimal point, 10 exponent fraction part, 11 blanks, 12 final comma or semicolon. with most quantities optional. Example: S R '1' 1.0 ' 1 ' 1.0 '1A' 1.0 '12,34,56' 12.0 ' 34 7' 34.0 '-1E2ABCD' -100.0 '-1X2ABCD' -1.0 ' 2E-1' 0.2 '23.45' 23.45 '-4.2E+2' -420.0 '17d2' 1700.0 '-14e-2' -0.14 'e2' 100.0 '-12.73e-9.23' -12.73 * 10.0**(-9.23) Licensing: This code is distributed under the MIT license. Modified: 24 June 2005 Author: John Burkardt Parameters: Input, char *S, the string containing the data to be read. Reading will begin at position 1 and terminate at the end of the string, or when no more characters can be read to form a legal real. Blanks, commas, or other nonnumeric data will, in particular, cause the conversion to halt. Output, int *LCHAR, the number of characters read from the string to form the number, including any terminating characters such as a trailing comma or blanks. Output, int *ERROR, is TRUE (1) if an error occurred and FALSE (0) otherwise. Output, double S_TO_R8, the value that was read from the string. */ { char c; int ihave; int isgn; int iterm; int jbot; int jsgn; int jtop; int nchar; int ndig; double r; double rbot; double rexp; double rtop; char TAB = 9; nchar = s_len_trim ( s ); *error = 0; r = 0.0; *lchar = -1; isgn = 1; rtop = 0.0; rbot = 1.0; jsgn = 1; jtop = 0; jbot = 1; ihave = 1; iterm = 0; for ( ; ; ) { c = s[*lchar+1]; *lchar = *lchar + 1; /* Blank or TAB character. */ if ( c == ' ' || c == TAB ) { if ( ihave == 2 ) { } else if ( ihave == 6 || ihave == 7 ) { iterm = 1; } else if ( 1 < ihave ) { ihave = 11; } } /* Comma. */ else if ( c == ',' || c == ';' ) { if ( ihave != 1 ) { iterm = 1; ihave = 12; *lchar = *lchar + 1; } } /* Minus sign. */ else if ( c == '-' ) { if ( ihave == 1 ) { ihave = 2; isgn = -1; } else if ( ihave == 6 ) { ihave = 7; jsgn = -1; } else { iterm = 1; } } /* Plus sign. */ else if ( c == '+' ) { if ( ihave == 1 ) { ihave = 2; } else if ( ihave == 6 ) { ihave = 7; } else { iterm = 1; } } /* Decimal point. */ else if ( c == '.' ) { if ( ihave < 4 ) { ihave = 4; } else if ( 6 <= ihave && ihave <= 8 ) { ihave = 9; } else { iterm = 1; } } /* Exponent marker. */ else if ( ch_eqi ( c, 'E' ) || ch_eqi ( c, 'D' ) ) { if ( ihave < 6 ) { ihave = 6; } else { iterm = 1; } } /* Digit. */ else if ( ihave < 11 && '0' <= c && c <= '9' ) { if ( ihave <= 2 ) { ihave = 3; } else if ( ihave == 4 ) { ihave = 5; } else if ( ihave == 6 || ihave == 7 ) { ihave = 8; } else if ( ihave == 9 ) { ihave = 10; } ndig = ch_to_digit ( c ); if ( ihave == 3 ) { rtop = 10.0 * rtop + ( double ) ndig; } else if ( ihave == 5 ) { rtop = 10.0 * rtop + ( double ) ndig; rbot = 10.0 * rbot; } else if ( ihave == 8 ) { jtop = 10 * jtop + ndig; } else if ( ihave == 10 ) { jtop = 10 * jtop + ndig; jbot = 10 * jbot; } } /* Anything else is regarded as a terminator. */ else { iterm = 1; } /* If we haven't seen a terminator, and we haven't examined the entire string, go get the next character. */ if ( iterm == 1 || nchar <= *lchar + 1 ) { break; } } /* If we haven't seen a terminator, and we have examined the entire string, then we're done, and LCHAR is equal to NCHAR. */ if ( iterm != 1 && (*lchar) + 1 == nchar ) { *lchar = nchar; } /* Number seems to have terminated. Have we got a legal number? Not if we terminated in states 1, 2, 6 or 7! */ if ( ihave == 1 || ihave == 2 || ihave == 6 || ihave == 7 ) { *error = 1; return r; } /* Number seems OK. Form it. We have had some trouble with input of the form 1.0E-312. For now, let's assume anything less than 1.0E-20 is zero. */ if ( jtop == 0 ) { rexp = 1.0; } else { if ( jbot == 1 ) { if ( jsgn * jtop < -20 ) { rexp = 0.0; } else { rexp = pow ( ( double ) 10.0, ( double ) ( jsgn * jtop ) ); } } else { if ( jsgn * jtop < -20 * jbot ) { rexp = 0.0; } else { rexp = jsgn * jtop; rexp = rexp / jbot; rexp = pow ( ( double ) 10.0, ( double ) rexp ); } } } r = isgn * rexp * rtop / rbot; return r; } /******************************************************************************/ int s_to_r8vec ( char *s, int n, double rvec[] ) /******************************************************************************/ /* Purpose: S_TO_R8VEC reads an R8VEC from a string. Licensing: This code is distributed under the MIT license. Modified: 19 February 2001 Author: John Burkardt Parameters: Input, char *S, the string to be read. Input, int N, the number of values expected. Output, double RVEC[N], the values read from the string. Output, int S_TO_R8VEC, is TRUE (1) if an error occurred and FALSE (0) otherwise. */ { int error; int i; int lchar; error = 0; for ( i = 0; i < n; i++ ) { rvec[i] = s_to_r8 ( s, &lchar, &error ); if ( error ) { return error; } s = s + lchar; } return error; } /******************************************************************************/ int s_word_count ( char *s ) /******************************************************************************/ /* Purpose: S_WORD_COUNT counts the number of "words" in a string. Licensing: This code is distributed under the MIT license. Modified: 16 September 2015 Author: John Burkardt Parameters: Input, char *S, the string to be examined. Output, int S_WORD_COUNT, the number of "words" in the string. Words are presumed to be separated by one or more blanks. */ { int blank; int word_num; char *t; word_num = 0; blank = 1; t = s; while ( *t ) { if ( *t == ' ' || *t == '\n' ) { blank = 1; } else if ( blank ) { word_num = word_num + 1; blank = 0; } t++; } return word_num; } /******************************************************************************/ 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; } /******************************************************************************/ 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 ); printf ( "%s\n", time_buffer ); return; # undef TIME_SIZE } /******************************************************************************/ int *triangulation_neighbor_triangles ( int triangle_order, int triangle_num, int triangle_node[] ) /******************************************************************************/ /* Purpose: TRIANGULATION_ORDER3_NEIGHBOR_TRIANGLES determines triangle neighbors. Discussion: A triangulation of a set of nodes can be completely described by the coordinates of the nodes, and the list of nodes that make up each triangle. However, in some cases, it is necessary to know triangle adjacency information, that is, which triangle, if any, is adjacent to a given triangle on a particular side. This routine creates a data structure recording this information. The primary amount of work occurs in sorting a list of 3 * TRIANGLE_NUM data items. This routine was modified to work with columns rather than rows. Example: The input information from TRIANGLE_NODE: Triangle Nodes -------- --------------- 1 3 4 1 2 3 1 2 3 3 2 8 4 2 1 5 5 8 2 13 6 8 13 9 7 3 8 9 8 13 2 5 9 9 13 7 10 7 13 5 11 6 7 5 12 9 7 6 13 10 9 6 14 6 5 12 15 11 6 12 16 10 6 11 The output information in TRIANGLE_NEIGHBOR: Triangle Neighboring Triangles -------- --------------------- 1 -1 -1 2 2 1 4 3 3 2 5 7 4 2 -1 8 5 3 8 6 6 5 9 7 7 3 6 -1 8 5 4 10 9 6 10 12 10 9 8 11 11 12 10 14 12 9 11 13 13 -1 12 16 14 11 -1 15 15 16 14 -1 16 13 15 -1 Licensing: This code is distributed under the MIT license. Modified: 08 October 2010 Author: John Burkardt Parameters: Input, int TRIANGLE_ORDER, the order of the triangles. Input, int TRIANGLE_NUM, the number of triangles. Input, int TRIANGLE_NODE[TRIANGLE_ORDER*TRIANGLE_NUM], the nodes that make up each triangle. Output, int TRIANGLE_ORDER3_NEIGHBOR_TRIANGLES[3*TRIANGLE_NUM], the three triangles that are direct neighbors of a given triangle. TRIANGLE_NEIGHBOR(1,I) is the index of the triangle which touches side 1, defined by nodes 2 and 3, and so on. TRIANGLE_NEIGHBOR(1,I) is negative if there is no neighbor on that side. In this case, that side of the triangle lies on the boundary of the triangulation. */ { int *col; int i; int icol; int j; int k; int side1; int side2; int tri; int tri1; int tri2; int *triangle_neighbor; triangle_neighbor = ( int * ) malloc ( 3 * triangle_num * sizeof ( int ) ); col = ( int * ) malloc ( 4 * 3 * triangle_num * sizeof ( int ) ); /* Step 1. From the list of nodes for triangle T, of the form: (I,J,K) construct the three neighbor relations: (I,J,3,T) or (J,I,3,T), (J,K,1,T) or (K,J,1,T), (K,I,2,T) or (I,K,2,T) where we choose (I,J,3,T) if I < J, or else (J,I,3,T) */ for ( tri = 0; tri < triangle_num; tri++ ) { i = triangle_node[0+tri*triangle_order]; j = triangle_node[1+tri*triangle_order]; k = triangle_node[2+tri*triangle_order]; if ( i < j ) { col[0+(3*tri+0)*4] = i; col[1+(3*tri+0)*4] = j; col[2+(3*tri+0)*4] = 3; col[3+(3*tri+0)*4] = tri + 1; } else { col[0+(3*tri+0)*4] = j; col[1+(3*tri+0)*4] = i; col[2+(3*tri+0)*4] = 3; col[3+(3*tri+0)*4] = tri + 1; } if ( j < k ) { col[0+(3*tri+1)*4] = j; col[1+(3*tri+1)*4] = k; col[2+(3*tri+1)*4] = 1; col[3+(3*tri+1)*4] = tri + 1; } else { col[0+(3*tri+1)*4] = k; col[1+(3*tri+1)*4] = j; col[2+(3*tri+1)*4] = 1; col[3+(3*tri+1)*4] = tri + 1; } if ( k < i ) { col[0+(3*tri+2)*4] = k; col[1+(3*tri+2)*4] = i; col[2+(3*tri+2)*4] = 2; col[3+(3*tri+2)*4] = tri + 1; } else { col[0+(3*tri+2)*4] = i; col[1+(3*tri+2)*4] = k; col[2+(3*tri+2)*4] = 2; col[3+(3*tri+2)*4] = tri + 1; } } /* Step 2. Perform an ascending dictionary sort on the neighbor relations. We only intend to sort on rows 1 and 2; the routine we call here sorts on rows 1 through 4 but that won't hurt us. What we need is to find cases where two triangles share an edge. Say they share an edge defined by the nodes I and J. Then there are two columns of COL that start out ( I, J, ?, ? ). By sorting COL, we make sure that these two columns occur consecutively. That will make it easy to notice that the triangles are neighbors. */ i4col_sort_a ( 4, 3*triangle_num, col ); /* Step 3. Neighboring triangles show up as consecutive columns with identical first two entries. Whenever you spot this happening, make the appropriate entries in TRIANGLE_NEIGHBOR. */ for ( j = 0; j < triangle_num; j++ ) { for ( i = 0; i < 3; i++ ) { triangle_neighbor[i+j*3] = -1; } } icol = 1; for ( ; ; ) { if ( 3 * triangle_num <= icol ) { break; } if ( col[0+(icol-1)*4] != col[0+icol*4] || col[1+(icol-1)*4] != col[1+icol*4] ) { icol = icol + 1; continue; } side1 = col[2+(icol-1)*4]; tri1 = col[3+(icol-1)*4]; side2 = col[2+ icol *4]; tri2 = col[3+ icol *4]; triangle_neighbor[side1-1+(tri1-1)*3] = tri2; triangle_neighbor[side2-1+(tri2-1)*3] = tri1; icol = icol + 2; } free ( col ); return triangle_neighbor; }