# 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 *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 ( int m, int n, int a[], char *title ); void i4mat_transpose_print_some ( int m, int n, int a[], int ilo, int jlo, int ihi, int jhi, char *title ); 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 ( int m, int n, double a[], char *title ); 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 timestamp ( ); void xml_mesh1d_write ( char *xml_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ); void xml_mesh2d_write ( char *xml_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ); void xml_mesh3d_write ( char *xml_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ); /******************************************************************************/ int main ( int argc, char *argv[] ) /******************************************************************************/ /* Purpose: MAIN is the main program for FEM_TO_XML. Discussion: FEM_TO_XML converts a 1D, 2D or 3D mesh from FEM to XML format. Usage: fem_to_xml prefix where 'prefix' is the common filename prefix: * 'prefix'_nodes.txt contains the node coordinates, * 'prefix'_elements.txt contains the element connectivity * 'prefix'.xml will contain the DOLFIN XML version of the data. Licensing: This code is distributed under the MIT license. Modified: 05 October 2014 Author: John Burkardt Reference: Anders Logg, Kent-Andre Mardal, Garth Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Lecture Notes in Computational Science and Engineering, Springer, 2011, ISBN13: 978-364223098 */ { char element_filename[255]; int *element_node; int element_num; int element_order; int m; char node_filename[255]; int node_num; double *node_x; char prefix[255]; char xml_filename[255]; timestamp ( ); printf ( "\n" ); printf ( "FEM_TO_XML\n" ); printf ( " C version:\n" ); printf ( " Convert a 1D, 2D or 3D mesh from FEM to DOLFIN XML format.\n" ); printf ( "\n" ); printf ( " Read \"prefix\"_nodes.txt, node coordinates.\n" ); printf ( " Read \"prefix\"_elements.txt, element node connectivity.\n" ); printf ( "\n" ); printf ( " Create \"prefix\".xml, a corresponding DOLFIN XML mesh file.\n" ); /* Get the filename prefix. */ if ( argc <= 1 ) { printf ( "\n" ); printf ( "FEM_TO_XML:\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 ( xml_filename, prefix ); strcat ( xml_filename, ".xml" ); /* Read the node data. */ r8mat_header_read ( node_filename, &m, &node_num ); printf ( "\n" ); printf ( " Read the header of \"%s\".\n", node_filename ); printf ( "\n" ); printf ( " Spatial dimension = %d\n", m ); printf ( " Number of nodes = %d\n", node_num ); node_x = r8mat_data_read ( node_filename, m, node_num ); printf ( "\n" ); printf ( " Read the data in \"%s\".\n", node_filename ); r8mat_transpose_print_some ( m, node_num, node_x, 1, 1, m, 5, " Portion of node coordinate data:" ); /* Read the element data. */ i4mat_header_read ( element_filename, &element_order, &element_num ); if ( m == 1 ) { if ( element_order == 2 ) { } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "FEM_TO_XML - Fatal error!\n" ); fprintf ( stderr, " 1D elements must use two vertices.\n" ); exit ( 1 ); } } else if ( m == 2 ) { if ( element_order == 3 ) { } else if ( element_order == 6 ) { } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "FEM_TO_XML - Fatal error!\n" ); fprintf ( stderr, " 2D elements must use 3 vertices.\n" ); exit ( 1 ); } } else if ( m == 3 ) { if ( element_order == 4 ) { } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "FEM_TO_XML - Fatal error!\n" ); fprintf ( stderr, " 3D elements must use 4vertices.\n" ); exit ( 1 ); } } printf ( "\n" ); printf ( " Read the header of \"%s\".\n", element_filename ); printf ( "\n" ); printf ( " Element order = %d\n", element_order ); printf ( " Number of elements = %d\n", element_num ); element_node = i4mat_data_read ( element_filename, element_order, element_num ); printf ( "\n" ); printf ( " Read the data in \"%s\".\n", element_filename ); i4mat_transpose_print_some ( element_order, element_num, element_node, 1, 1, element_order, 10, " Initial portion of element data:" ); /* Write the XML version of the data. */ if ( m == 1 ) { xml_mesh1d_write ( xml_filename, m, node_num, node_x, element_order, element_num, element_node ); } else if ( m == 2 ) { xml_mesh2d_write ( xml_filename, m, node_num, node_x, element_order, element_num, element_node ); } else if ( m == 3 ) { xml_mesh3d_write ( xml_filename, m, node_num, node_x, element_order, element_num, element_node ); } printf ( "\n" ); printf ( " Created XML file \"%s\".\n", xml_filename ); /* Free memory. */ free ( element_node ); free ( node_x ); /* Terminate. */ printf ( "\n" ); printf ( "FEM_TO_XML:\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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILE_COLUMN_COUNT - Fatal error!\n" ); fprintf ( stderr, " Could not open the input file: \"%s\"\n", input_filename ); exit ( 1 ); } /* 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILE_COLUMN_COUNT - Warning!\n" ); fprintf ( stderr, " The file does not seem to contain any data.\n" ); exit ( 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "FILE_ROW_COUNT - Fatal error!\n" ); fprintf ( stderr, " Could not open the input file: \"%s\"\n", input_filename ); exit ( 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 *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 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4MAT_DATA_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open the input file: \"%s\"\n", input_filename ); exit ( 1 ); } 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4MAT_HEADER_READ - Fatal error!\n" ); fprintf ( stderr, " FILE_COLUMN_COUNT failed.\n" ); exit ( 1 ); } *n = file_row_count ( input_filename ); if ( *n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4MAT_HEADER_READ - Fatal error!\n" ); fprintf ( stderr, " FILE_ROW_COUNT failed.\n" ); exit ( 1 ); } 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 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. Licensing: This code is distributed under the MIT license. Modified: 15 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 i; const int i4_huge = 2147483647; int j; int node_max; int node_min; int t; node_min = + i4_huge; node_max = - i4_huge; for ( j = 0; j < element_num; j++ ) { for ( i = 0; i < element_order; i++ ) { t = element_node[i+j*element_order]; if ( t < node_min ) { node_min = t; } if ( node_max < t ) { node_max = t; } } } 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 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 ( j = 0; j < element_num; j++ ) { for ( i = 0; i < element_order; i++ ) { element_node[i+j*element_order] = element_node[i+j*element_order] - 1; } } } 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 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8MAT_DATA_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open the input file: \"%s\"\n", input_filename ); exit ( 1 ); } 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 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8MAT_HEADER_READ - Fatal error!\n" ); fprintf ( stderr, " FILE_COLUMN_COUNT failed.\n" ); exit ( 1 ); } *n = file_row_count ( input_filename ); if ( *n <= 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "R8MAT_HEADER_READ - Fatal error!\n" ); fprintf ( stderr, " FILE_ROW_COUNT failed.\n" ); exit ( 1 ); } return; } /******************************************************************************/ 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 } /******************************************************************************/ 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 != ' ' && *t != '\n' ) { 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 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 } /******************************************************************************/ void xml_mesh1d_write ( char *xml_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: XML_MESH1D_WRITE writes a 1D mesh as a DOLFIN XML file. Licensing: This code is distributed under the MIT license. Modified: 05 October 2014 Author: John Burkardt Reference: Anders Logg, Kent-Andre Mardal, Garth Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Lecture Notes in Computational Science and Engineering, Springer, 2011, ISBN13: 978-364223098 Parameters: Input, char *XML_FILENAME, the name of the XML file to create. Input, int M, the spatial dimension. Input, int NODE_NUM, the number of nodes. Input, double NODE_X[M*NODE_NUM], the node coordinates. Input, int ELEMENT_ORDER, the order of the elements. Input, int ELEMENT_NUM, the number of elements. Input, int ELEMENT_NODE[ELEMENT_ORDER*ELEMENT_NUM], the nodes that make up each element. */ { int element; FILE *xml; int node; /* Force 0-based indexing. */ mesh_base_zero ( node_num, element_order, element_num, element_node ); /* Open the file. */ xml = fopen ( xml_filename, "wt" ); /* Write the data. */ fprintf ( xml, "\n" ); fprintf ( xml, "\n" ); fprintf ( xml, "\n" ); fprintf ( xml, " \n", m ); fprintf ( xml, " \n", node_num ); for ( node = 0; node < node_num; node++ ) { fprintf ( xml, " \n", node, node_x[0+node*m] ); } fprintf ( xml, " \n" ); fprintf ( xml, " \n", element_num ); for ( element = 0; element < element_num; element++ ) { fprintf ( xml, " \n", element, element_node[0+element*element_order], element_node[1+element*element_order] ); } fprintf ( xml, " \n" ); fprintf ( xml, " \n" ); fprintf ( xml, "\n" ); fclose ( xml ); return; } /******************************************************************************/ void xml_mesh2d_write ( char *xml_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: XML_MESH2D_WRITE writes a 2D mesh as a DOLFIN XML file. Licensing: This code is distributed under the MIT license. Modified: 26 May 2014 Author: John Burkardt Reference: Anders Logg, Kent-Andre Mardal, Garth Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Lecture Notes in Computational Science and Engineering, Springer, 2011, ISBN13: 978-364223098 Parameters: Input, char *XML_FILENAME, the name of the XML file to create. Input, int M, the spatial dimension. Input, int NODE_NUM, the number of nodes. Input, double NODE_X[M*NODE_NUM], the node coordinates. Input, int ELEMENT_ORDER, the order of the elements. Input, int ELEMENT_NUM, the number of elements. Input, int ELEMENT_NODE[ELEMENT_ORDER*ELEMENT_NUM], the nodes that make up each element. */ { int element; FILE *xml; int node; /* Force 0-based indexing. */ mesh_base_zero ( node_num, element_order, element_num, element_node ); /* Open the file. */ xml = fopen ( xml_filename, "wt" ); /* Write the data. */ fprintf ( xml, "\n" ); fprintf ( xml, "\n" ); fprintf ( xml, "\n" ); fprintf ( xml, " \n", m ); fprintf ( xml, " \n", node_num ); for ( node = 0; node < node_num; node++ ) { fprintf ( xml, " \n", node, node_x[0+node*m], node_x[1+node*m] ); } fprintf ( xml, " \n" ); fprintf ( xml, " \n", element_num ); for ( element = 0; element < element_num; element++ ) { fprintf ( xml, " \n", element, element_node[0+element*element_order], element_node[1+element*element_order], element_node[2+element*element_order] ); } fprintf ( xml, " \n" ); fprintf ( xml, " \n" ); fprintf ( xml, "\n" ); fclose ( xml ); return; } /******************************************************************************/ void xml_mesh3d_write ( char *xml_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: XML_MESH3D_WRITE writes a 3D mesh as a DOLFIN XML file. Licensing: This code is distributed under the MIT license. Modified: 05 October 2014 Author: John Burkardt Reference: Anders Logg, Kent-Andre Mardal, Garth Wells, Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Lecture Notes in Computational Science and Engineering, Springer, 2011, ISBN13: 978-364223098 Parameters: Input, char *XML_FILENAME, the name of the XML file to create. Input, int M, the spatial dimension. Input, int NODE_NUM, the number of nodes. Input, double NODE_X[M*NODE_NUM], the node coordinates. Input, int ELEMENT_ORDER, the order of the elements. Input, int ELEMENT_NUM, the number of elements. Input, int ELEMENT_NODE[ELEMENT_ORDER*ELEMENT_NUM], the nodes that make up each element. */ { int element; FILE *xml; int node; /* Force 0-based indexing. */ mesh_base_zero ( node_num, element_order, element_num, element_node ); /* Open the file. */ xml = fopen ( xml_filename, "wt" ); /* Write the data. */ fprintf ( xml, "\n" ); fprintf ( xml, "\n" ); fprintf ( xml, "\n" ); fprintf ( xml, " \n", m ); fprintf ( xml, " \n", node_num ); for ( node = 0; node < node_num; node++ ) { fprintf ( xml, " \n", node, node_x[0+node*m], node_x[1+node*m], node_x[2+node*m] ); } fprintf ( xml, " \n" ); fprintf ( xml, " \n", element_num ); for ( element = 0; element < element_num; element++ ) { fprintf ( xml, " \n", element, element_node[0+element*element_order], element_node[1+element*element_order], element_node[2+element*element_order], element_node[3+element*element_order] ); } fprintf ( xml, " \n" ); fprintf ( xml, " \n" ); fprintf ( xml, "\n" ); fclose ( xml ); return; }