# include # include # include # include # include # include "gmsh_io.h" /******************************************************************************/ 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; } /******************************************************************************/ void gmsh_data_read ( char *gmsh_filename, int node_dim, int node_num, double node_x[], int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: GMSH_DATA_READ reads data from a GMSH file. Licensing: This code is distributed under the MIT license. Modified: 20 October 2014 Author: John Burkardt Parameters: Input, character *GMSH_FILENAME, the GMSH filename. Input, int NODE_DIM, the spatial dimension. Input, int NODE_NUM, the number of nodes. Input, double NODE_X[NODE_DIM*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. */ { char *error; int i; int ierror; FILE *input; int j; int k; int length; int level; char text[255]; char* text_pointer; double x; input = fopen ( gmsh_filename, "rt" ); if ( ! input ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GMSH_DATA_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open input file \"%s\"\n", gmsh_filename ); exit ( 1 ); } level = 0; for ( ; ; ) { text_pointer = text; error = fgets ( text_pointer, 255, input ); if ( !error ) { break; } if ( level == 0 ) { if ( s_begin ( text_pointer, "$Nodes" ) ) { level = 1; j = 0; } } else if ( level == 1 ) { s_to_i4 ( text_pointer, &length, &ierror ); level = 2; } else if ( level == 2 ) { if ( s_begin ( text_pointer, "$EndNodes" ) ) { break; } else { s_to_i4 ( text_pointer, &length, &ierror ); text_pointer = text_pointer + length; for ( i = 0; i < node_dim; i++ ) { x = s_to_r8 ( text_pointer, &length, &ierror ); text_pointer = text_pointer + length; node_x[i+j*node_dim] = x; } j = j + 1; } } } /* Now read element information. */ level = 0; for ( ; ; ) { text_pointer = text; error = fgets ( text_pointer, 255, input ); if ( !error ) { break; } if ( level == 0 ) { if ( s_begin ( text_pointer, "$Elements" ) ) { level = 1; j = 0; } } else if ( level == 1 ) { s_to_i4 ( text_pointer, &length, &ierror ); level = 2; } else if ( level == 2 ) { if ( s_begin ( text_pointer, "$EndElements" ) ) { break; } else { for ( k = 1; k <= 5; k++ ) { s_to_i4 ( text_pointer, &length, &ierror ); text_pointer = text_pointer + length; } for ( i = 0; i < element_order; i++ ) { k = s_to_i4 ( text_pointer, &length, &ierror ); text_pointer = text_pointer + length; element_node[i+j*element_order] = k; } j = j + 1; } } } fclose ( input ); return; } /******************************************************************************/ void gmsh_size_read ( char *gmsh_filename, int *node_num, int *node_dim, int *element_num, int *element_order ) /******************************************************************************/ /* Purpose: GMSH_SIZE_READ reads sizes from a GMSH file. Licensing: This code is distributed under the MIT license. Modified: 19 October 2014 Author: John Burkardt Parameters: Input, character *GMSH_FILENAME, the GMSH filename. Output, int *NODE_NUM, the number of nodes. Output, int *NODE_DIM, the spatial dimension. Output, int *ELEMENT_NUM, the number of elements. Output, int *ELEMENT_ORDER, the order of the elements. */ { char *error; int ierror; FILE *input; int k; int length; int level; const double r8_big = 1.0E+30; char text[255]; char* text_pointer; double x; double x_max; double x_min; double y; double y_max; double y_min; double z; double z_max; double z_min; *node_num = 0; *node_dim = 0; x_max = - r8_big; x_min = + r8_big; y_max = - r8_big; y_min = + r8_big; z_max = - r8_big; z_min = + r8_big; input = fopen ( gmsh_filename, "rt" ); if ( ! input ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "GMSH_SIZE_READ - Fatal error!\n" ); fprintf ( stderr, " Could not open input file \"%s\"\n", gmsh_filename ); exit ( 1 ); } level = 0; for ( ; ; ) { text_pointer = text; error = fgets ( text_pointer, 255, input ); if ( !error ) { break; } if ( level == 0 ) { if ( s_begin ( text_pointer, "$Nodes" ) ) { level = 1; } } else if ( level == 1 ) { *node_num = s_to_i4 ( text_pointer, &length, &ierror ); level = 2; } else if ( level == 2 ) { if ( s_begin ( text_pointer, "$EndNodes" ) ) { break; } else { s_to_i4 ( text_pointer, &length, &ierror ); text_pointer = text_pointer + length; x = s_to_r8 ( text_pointer, &length, &ierror ); x_min = fmin ( x_min, x ); x_max = fmax ( x_max, x ); text_pointer = text_pointer + length; y = s_to_r8 ( text_pointer, &length, &ierror ); y_min = fmin ( y_min, y ); y_max = fmax ( y_max, y ); text_pointer = text_pointer + length; z = s_to_r8 ( text_pointer, &length, &ierror); text_pointer = text_pointer + length; z_min = fmin ( z_min, z ); z_max = fmax ( z_max, z ); } } } /* Make a very simple guess as to the dimensionality of the data. */ *node_dim = 3; if ( z_max == z_min ) { *node_dim = 2; if ( y_max == y_min ) { *node_dim = 1; } } /* Now read element information. */ level = 0; for ( ; ; ) { text_pointer = text; error = fgets ( text_pointer, 255, input ); if ( !error ) { break; } if ( level == 0 ) { if ( s_begin ( text_pointer, "$Elements" ) ) { level = 1; } } else if ( level == 1 ) { *element_num = s_to_i4 ( text_pointer, &length, &ierror ); level = 2; } else if ( level == 2 ) { if ( s_begin ( text_pointer, "$EndElements" ) ) { break; } else { k = 0; for ( ; ; ) { s_to_i4 ( text_pointer, &length, &ierror ); text_pointer = text_pointer + length; if ( ierror != 0 ) { break; } k = k + 1; } *element_order = k - 5; break; } } } fclose ( input ); return; } /******************************************************************************/ int *gmsh_mesh2d_element_data_example ( int element_num, int element_order ) /******************************************************************************/ /* Purpose: GMSH_MESH2D_ELEMENT_DATA_EXAMPLE: element information for the example. Licensing: This code is distributed under the MIT license. Modified: 18 October 2014 Author: John Burkardt Parameters: Input, int ELEMENT_NUM, the number of elements. Input, int ELEMENT_ORDER, the order of the elements. Output, int GMSH_MESH2D_ELEMENT_DATA_EXAMPLE[ELEMENT_ORDER*ELEMENT_NUM], the indices of the nodes that make up each element. */ { int *element_node; int element_node_save[3*24] = { 1, 2, 6, 7, 6, 2, 2, 3, 7, 8, 7, 3, 3, 4, 8, 9, 8, 4, 4, 5, 9, 10, 9, 5, 6, 7, 11, 12, 11, 7, 7, 8, 12, 13, 12, 8, 8, 9, 13, 14, 13, 9, 9, 10, 14, 15, 14, 10, 11, 12, 16, 17, 16, 12, 12, 13, 17, 18, 17, 13, 16, 17, 19, 20, 19, 17, 17, 18, 20, 21, 20, 18 }; element_node = i4mat_copy_new ( element_order, element_num, element_node_save ); return element_node; } /******************************************************************************/ void gmsh_mesh2d_element_size_example ( int *element_num, int *element_order ) /******************************************************************************/ /* Purpose: GMSH_MESH2D_ELEMENT_SIZE_EXAMPLE: element size information for the example. Licensing: This code is distributed under the MIT license. Modified: 18 October 2014 Author: John Burkardt Parameters: Output, int *ELEMENT_NUM, the number of elements. Output, int *ELEMENT_ORDER, the order of the elements. */ { *element_num = 24; *element_order = 3; return; } /******************************************************************************/ double *gmsh_mesh2d_node_data_example ( int node_num, int node_dim ) /******************************************************************************/ /* Purpose: GMSH_MESH2D_NODE_DATA_EXAMPLE returns the node information for the example. 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 NODE_DIM, the spatial dimension. Output, double GMSH_MESH2D_NODE_DATA_EXAMPLE[NODE_DIM*NODE_NUM], the nodal coordinates. */ { double *node_coord; double node_coord_save[2*21] = { 0.0, 0.0, 1.0, 0.0, 2.0, 0.0, 3.0, 0.0, 4.0, 0.0, 0.0, 1.0, 1.0, 1.0, 2.0, 1.0, 3.0, 1.0, 4.0, 1.0, 0.0, 2.0, 1.0, 2.0, 2.0, 2.0, 3.0, 2.0, 4.0, 2.0, 0.0, 3.0, 1.0, 3.0, 2.0, 3.0, 0.0, 4.0, 1.0, 4.0, 2.0, 4.0 }; node_coord = r8mat_copy_new ( 2, 21, node_coord_save ); return node_coord; } /******************************************************************************/ void gmsh_mesh2d_node_size_example ( int *node_num, int *node_dim ) /******************************************************************************/ /* Purpose: GMSH_MESH2D_NODE_SIZE_EXAMPLE: sizes of node information for the example. Licensing: This code is distributed under the MIT license. Modified: 18 October 2014 Author: John Burkardt Parameters: Output, int *NODE_NUM, the number of nodes. Output, int *NODE_DIM, the spatial dimension. */ { *node_num = 21; *node_dim = 2; return; } /******************************************************************************/ void gmsh_mesh1d_write ( char *gmsh_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: GMSH_MESH1D_WRITE writes 1D mesh data as a Gmsh mesh file. Licensing: This code is distributed under the MIT license. Modified: 06 October 2014 Author: John Burkardt Reference: Christophe Geuzaine, Jean-Francois Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering, Volume 79, Number 11, pages 1309-1331, 2009. Parameters: Input, char *GMSH_FILENAME, the name of the Gmsh file. Input, int M, the spatial dimension. Input, inte 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; int element_type; FILE *gmsh; int i; int node; int tag_num; int tag1; /* Force 1-based indexing. */ mesh_base_one ( node_num, element_order, element_num, element_node ); /* Open the file. */ gmsh = fopen ( gmsh_filename, "wt" ); /* Write the data. */ fprintf ( gmsh, "$MeshFormat\n" ); fprintf ( gmsh, "2.2 0 8\n" ); fprintf ( gmsh, "$EndMeshFormat\n" ); fprintf ( gmsh, "$Nodes\n" ); fprintf ( gmsh, "%d\n", node_num ); for ( node = 0; node < node_num; node++ ) { fprintf ( gmsh, " %d %g 0.0 0.0\n", node + 1, node_x[0+node*m] ); } fprintf ( gmsh, "$EndNodes\n" ); element_type = 1; tag_num = 2; tag1 = 0; fprintf ( gmsh, "$Elements\n" ); fprintf ( gmsh, "%d\n", element_num ); for ( element = 0; element < element_num; element++ ) { fprintf ( gmsh, " %d %d %d %d %d", element + 1, element_type, tag_num, tag1, element + 1 ); for ( i = 0; i < element_order; i++ ) { fprintf ( gmsh, " %d", element_node[i+element*element_order] ); } fprintf ( gmsh, "\n" ); } fprintf ( gmsh, "$EndElements\n" ); fclose ( gmsh ); return; } /******************************************************************************/ void gmsh_mesh2d_write ( char *gmsh_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: GMSH_MESH2D_WRITE writes 2D mesh data as a Gmsh mesh file. Licensing: This code is distributed under the MIT license. Modified: 06 October 2014 Author: John Burkardt Reference: Christophe Geuzaine, Jean-Francois Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering, Volume 79, Number 11, pages 1309-1331, 2009. Parameters: Input, char *GMSH_FILENAME, the name of the Gmsh file. Input, int M, the spatial dimension. Input, inte 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; int element_type; FILE *gmsh; int i; int node; int tag_num; int tag1; /* Force 1-based indexing. */ mesh_base_one ( node_num, element_order, element_num, element_node ); /* Open the file. */ gmsh = fopen ( gmsh_filename, "wt" ); /* Write the data. */ fprintf ( gmsh, "$MeshFormat\n" ); fprintf ( gmsh, "2.2 0 8\n" ); fprintf ( gmsh, "$EndMeshFormat\n" ); fprintf ( gmsh, "$Nodes\n" ); fprintf ( gmsh, "%d\n", node_num ); for ( node = 0; node < node_num; node++ ) { fprintf ( gmsh, " %d %g %g 0.0\n", node + 1, node_x[0+node*2], node_x[1+node*2] ); } fprintf ( gmsh, "$EndNodes\n" ); if ( element_order == 3 ) { element_type = 2; } else if ( element_order == 6 ) { element_type = 9; } tag_num = 2; tag1 = 0; fprintf ( gmsh, "$Elements\n" ); fprintf ( gmsh, "%d\n", element_num ); for ( element = 0; element < element_num; element++ ) { fprintf ( gmsh, " %d %d %d %d %d", element + 1, element_type, tag_num, tag1, element + 1 ); for ( i = 0; i < element_order; i++ ) { fprintf ( gmsh, " %d", element_node[i+element*element_order] ); } fprintf ( gmsh, "\n" ); } fprintf ( gmsh, "$EndElements\n" ); fclose ( gmsh ); return; } /******************************************************************************/ void gmsh_mesh3d_write ( char *gmsh_filename, int m, int node_num, double node_x[], int element_order, int element_num, int element_node[] ) /******************************************************************************/ /* Purpose: GMSH_MESH3D_WRITE writes 3D mesh data as a Gmsh mesh file. Discussion: The node ordering for the 20 node element is not standard. Assuming the vertices are A, B, C and D, Gmsh uses the following ordering: 1: a 2: b 3: c 4: d 5: (2*a +b )/3 6: ( a+2*b )/3 7: ( 2*b+ c )/3 8: ( b+2*c )/3 9: ( a +2*c )/3 10: (2*a + c )/3 11: (2*a + d)/3 12: ( a +2*d)/3 13: ( b +2*d)/3 14: ( 2*b + d)/3 15: ( + c+2*d)/3 16: ( +2*c+ d)/3 17: ( a+ b+ c )/3 18: ( a+ b + d)/3 19: ( b+ c+ d)/3 20: ( a+ c+ d)/3 Leo Rebholz used the following ordering: 1: a 2: b 3: c 4: d 5: (2*a +b )/3 6: (2*a + c )/3 7: ( a+2*b )/3 8: ( a +2*c )/3 9: ( a+ b+ c )/3 10: ( 2*b+ c )/3 11: ( b+2*c )/3 12: (2*a + d)/3 13: ( 2*b + d)/3 14: ( +2*c+ d)/3 15: ( a+ b + d)/3 16: ( b+ c+ d)/3 17: ( a+ c+ d)/3 18: ( a +2*d)/3 19: ( b +2*d)/3 20: ( + c+2*d)/3 Since the only 20 node data we have is from Leo, we will assume that all 20 node input data is in Leo's format, and needs to be converted to the Gmsh convention. Licensing: This code is distributed under the MIT license. Modified: 06 October 2014 Author: John Burkardt Reference: Christophe Geuzaine, Jean-Francois Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering, Volume 79, Number 11, pages 1309-1331, 2009. Parameters: Input, char *GMSH_FILENAME, the name of the Gmsh file. 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; int element_type; FILE *gmsh; int i; int i2; int leo_to_gmsh[20] = { 0, 1, 2, 3, 4, 6, 9, 10, 7, 5, 11, 17, 18, 12, 19, 13, 8, 14, 15, 16 }; int node; int tag_num; int tag1; /* Force 1-based indexing. */ mesh_base_one ( node_num, element_order, element_num, element_node ); /* Open the file. */ gmsh = fopen ( gmsh_filename, "wt" ); /* Write the data. */ fprintf ( gmsh, "$MeshFormat\n" ); fprintf ( gmsh, "2.2 0 8\n" ); fprintf ( gmsh, "$EndMeshFormat\n" ); fprintf ( gmsh, "$Nodes\n" ); fprintf ( gmsh, "%d\n", node_num ); for ( node = 0; node < node_num; node++ ) { fprintf ( gmsh, "%d %g %g %g\n", node + 1, node_x[0+node*3], node_x[1+node*3], node_x[2+node*3] ); } fprintf ( gmsh, "$EndNodes\n" ); if ( element_order == 4 ) { element_type = 4; } else if ( element_order == 10 ) { element_type = 11; } else if ( element_order == 20 ) { element_type = 29; } tag_num = 2; tag1 = 0; fprintf ( gmsh, "$Elements\n" ); fprintf ( gmsh, "%d\n", element_num ); for ( element = 0; element < element_num; element++ ) { fprintf ( gmsh, "%d %d %d %d %d", element + 1, element_type, tag_num, tag1, element + 1 ); if ( element_order == 20 ) { for ( i = 0; i < element_order; i++ ) { i2 = leo_to_gmsh[i]; fprintf ( gmsh, " %d", element_node[i2+element*element_order] ); } fprintf ( gmsh, "\n" ); } else { for ( i = 0; i < element_order; i++ ) { fprintf ( gmsh, " %d", element_node[i+element*element_order] ); } fprintf ( gmsh, "\n" ); } } fprintf ( gmsh, "$EndElements\n" ); fclose ( gmsh ); return; } /******************************************************************************/ int *i4mat_copy_new ( int m, int n, int a1[] ) /******************************************************************************/ /* Purpose: I4MAT_COPY_NEW copies an I4MAT to a "new" I4MAT. 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: 27 August 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, int A1[M*N], the matrix to be copied. Output, int I4MAT_COPY_NEW[M*N], the copy of A1. */ { int *a2; int i; int j; a2 = ( int * ) malloc ( m * n * sizeof ( int ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a2[i+j*m] = a1[i+j*m]; } } return a2; } /******************************************************************************/ 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_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; } /******************************************************************************/ double *r8mat_copy_new ( int m, int n, double a1[] ) /******************************************************************************/ /* Purpose: R8MAT_COPY_NEW copies one R8MAT to a "new" 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 July 2008 Author: John Burkardt Parameters: Input, int M, N, the number of rows and columns. Input, double A1[M*N], the matrix to be copied. Output, double R8MAT_COPY_NEW[M*N], the copy of A1. */ { double *a2; int i; int j; a2 = ( double * ) malloc ( m * n * sizeof ( double ) ); for ( j = 0; j < n; j++ ) { for ( i = 0; i < m; i++ ) { a2[i+j*m] = a1[i+j*m]; } } return a2; } /******************************************************************************/ 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_begin ( char *s1, char *s2 ) /******************************************************************************/ /* Purpose: S_BEGIN reports whether string 1 begins with string 2. Licensing: This code is distributed under the MIT license. Modified: 22 November 2010 Author: John Burkardt Parameters: Input, char *S1, char *S2, two strings. Output, int S_BEGIN, is true if S1 is the same as S2 up to the end of S2, and false otherwise. */ { int i; int n1; int n2; n1 = strlen ( s1 ); n2 = strlen ( s2 ); if ( n1 < n2 ) { return 0; } for ( i = 0; i < n2; i++ ) { if ( ch_cap ( s1[i] ) != ch_cap ( s2[i] ) ) { return 0; } } return 1; } /******************************************************************************/ 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; } /******************************************************************************/ 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; } /******************************************************************************/ 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; } /******************************************************************************/ 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 }