program main !*****************************************************************************80 ! !! MAIN is the main program for TRIANGULATION_Q2L. ! ! Discussion: ! ! TRIANGULATION_Q2L makes a linear triangulation from a quadratic one. ! ! A quadratic triangulation is assumed to consist of 6-node triangles, ! as in the following: ! ! 11-12-13-14-15 ! |\ |\ | ! | \ | \ | ! 6 7 8 9 10 ! | \ | \ | ! | \| \| ! 1--2--3--4--5 ! ! This routine rearranges information so as to define the 3-node ! triangulation: ! ! 11-12-13-14-15 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 6--7--8--9-10 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 1--2--3--4--5 ! ! Usage: ! ! triangulation_q2l prefix ! ! where 'prefix' is the common filename prefix: ! ! * prefix_nodes.txt contains the node coordinates (not needed), ! * prefix_elements.txt contains the element definitions. ! * prefix_q2l_elements.txt will contain the linearized element definitions. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 October 2009 ! ! Author: ! ! John Burkardt ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer arg_num integer iarg integer iargc character ( len = 255 ) :: element_filename = ' ' character ( len = 255 ) :: triangle_q2l_filename = ' ' character ( len = 255 ) prefix integer, allocatable, dimension ( :, : ) :: triangle_node1 integer, allocatable, dimension ( :, : ) :: triangle_node2 integer triangle_num1 integer triangle_num2 integer triangle_order1 integer triangle_order2 call timestamp ( ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TRIANGULATION_Q2L' write ( *, '(a)' ) ' FORTRAN90 version' write ( *, '(a)' ) ' Read a "quadratic" triangulation and' write ( *, '(a)' ) ' write out a "linear" triangulation.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Read a triangulation dataset of TRIANGLE_NUM1' write ( *, '(a)' ) ' triangles using 6 nodes.' write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Create a 3 node triangulation by breaking' write ( *, '(a)' ) ' every 6 node triangle into 4 smaller ones.' write ( *, '(a)' ) ' Write the new linear triangulation to a file.' ! ! Get the number of command line arguments. ! arg_num = iargc ( ) ! ! Argument 1 is the common filename prefix. ! if ( 1 <= arg_num ) then iarg = 1 call getarg ( iarg, prefix ) else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TRIANGULATION_Q2l:' write ( *, '(a)' ) ' Please enter the filename prefix.' read ( *, '(a)' ) prefix end if ! ! Create the filenames. ! element_filename = trim ( prefix ) // '_elements.txt' triangle_q2l_filename = trim ( prefix ) // '_q2l_elements.txt' ! ! Read the element data. ! call i4mat_header_read ( element_filename, triangle_order1, & triangle_num1 ) if ( triangle_order1 /= 6 ) then write ( *, * ) ' ' write ( *, '(a)' ) 'TRIANGULATION_Q2L - Fatal error!' write ( *, '(a)' ) ' Input data is not for a 6-node triangulation.' stop end if write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Read the header of "' & // trim ( element_filename ) //'".' write ( *, '(a)' ) ' ' write ( *, '(a,i8)' ) ' Triangle order = ', triangle_order1 write ( *, '(a,i8)' ) ' Number of triangles TRIANGLE_NUM1 = ', triangle_num1 allocate ( triangle_node1(1:triangle_order1,1:triangle_num1) ) call i4mat_data_read ( element_filename, triangle_order1, & triangle_num1, triangle_node1 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Read the data in "' & // trim ( element_filename ) //'".' call i4mat_transpose_print_some ( triangle_order1, triangle_num1, & triangle_node1, & 1, 1, 6, 10, ' 6 by 10 portion TRIANGLE_NODE1:' ) ! ! Set the number of linear triangles: ! triangle_num2 = 4 * triangle_num1 triangle_order2 = 3 ! ! Allocate space. ! allocate ( triangle_node2(1:triangle_order2,1:triangle_num2) ) ! ! Convert the data. ! call triangulation_order6_to_order3 ( triangle_num1, triangle_node1, & triangle_num2, triangle_node2 ) call i4mat_transpose_print ( triangle_order2, triangle_num2, & triangle_node2, ' TRIANGLE_NODE2' ) ! ! Write out the triangle data for the quadratic mesh. ! call i4mat_write ( triangle_q2l_filename, triangle_order2, & triangle_num2, triangle_node2 ) write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' Wrote the linearized element data to "' & // trim ( triangle_q2l_filename ) //'".' ! ! Free memory. ! deallocate ( triangle_node1 ) deallocate ( triangle_node2 ) ! ! Terminate. ! write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'TRIANGULATION_Q2L' write ( *, '(a)' ) ' Normal end of execution.' write ( *, '(a)' ) ' ' call timestamp ( ) stop 0 end subroutine ch_cap ( c ) !*****************************************************************************80 ! !! CH_CAP capitalizes a single character. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 19 July 1998 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input/output, character C, the character to capitalize. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character c integer itemp itemp = ichar ( c ) if ( 97 <= itemp .and. itemp <= 122 ) then c = char ( itemp - 32 ) end if return end function ch_eqi ( c1, c2 ) !*****************************************************************************80 ! !! CH_EQI is a case insensitive comparison of two characters for equality. ! ! Example: ! ! CH_EQI ( 'A', 'a' ) is .TRUE. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 July 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C1, C2, the characters to compare. ! ! Output, logical CH_EQI, the result of the comparison. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) logical ch_eqi character c1 character c1_cap character c2 character c2_cap c1_cap = c1 c2_cap = c2 call ch_cap ( c1_cap ) call ch_cap ( c2_cap ) if ( c1_cap == c2_cap ) then ch_eqi = .true. else ch_eqi = .false. end if return end subroutine ch_to_digit ( c, digit ) !*****************************************************************************80 ! !! CH_TO_DIGIT returns the value of a base 10 digit. ! ! Example: ! ! C DIGIT ! --- ----- ! '0' 0 ! '1' 1 ! ... ... ! '9' 9 ! ' ' 0 ! 'X' -1 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 August 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character C, the decimal digit, '0' through '9' or blank ! are legal. ! ! Output, integer DIGIT, the corresponding value. If C was ! 'illegal', then DIGIT is -1. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character c integer digit if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then digit = ichar ( c ) - 48 else if ( c == ' ' ) then digit = 0 else digit = -1 end if return end subroutine file_column_count ( input_filename, column_num ) !*****************************************************************************80 ! !! 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: ! ! 21 June 2001 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the file. ! ! Output, integer COLUMN_NUM, the number of columns in the file. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer column_num logical got_one character ( len = * ) input_filename integer input_unit integer ios character ( len = 256 ) line ! ! Open the file. ! call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & form = 'formatted', access = 'sequential', iostat = ios ) if ( ios /= 0 ) then column_num = -1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Fatal error!' write ( *, '(a)' ) ' Could not open the file:' write ( *, '(a)' ) ' ' // trim ( input_filename ) return end if ! ! Read one line, but skip blank lines and comment lines. ! got_one = .false. do read ( input_unit, '(a)', iostat = ios ) line if ( ios /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if if ( line(1:1) == '#' ) then cycle end if got_one = .true. exit end do if ( .not. got_one ) then rewind ( input_unit ) do read ( input_unit, '(a)', iostat = ios ) line if ( ios /= 0 ) then exit end if if ( len_trim ( line ) == 0 ) then cycle end if got_one = .true. exit end do end if close ( unit = input_unit ) if ( .not. got_one ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_COLUMN_COUNT - Warning!' write ( *, '(a)' ) ' The file does not seem to contain any data.' column_num = -1 return end if call s_word_count ( line, column_num ) return end subroutine file_row_count ( input_filename, row_num ) !*****************************************************************************80 ! !! 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: ! ! 06 March 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output, integer ROW_NUM, the number of rows found. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer bad_num integer comment_num integer ierror character ( len = * ) input_filename integer input_unit integer ios character ( len = 100 ) line integer record_num integer row_num call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = ios ) if ( ios /= 0 ) then row_num = -1; ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'FILE_ROW_COUNT - Fatal error!' write ( *, '(a)' ) ' Could not open the input file: ' // & trim ( input_filename ) stop end if comment_num = 0 row_num = 0 record_num = 0 bad_num = 0 do read ( input_unit, '(a)', iostat = ios ) line if ( ios /= 0 ) then ierror = record_num exit end if record_num = record_num + 1 if ( line(1:1) == '#' ) then comment_num = comment_num + 1 cycle end if if ( len_trim ( line ) == 0 ) then comment_num = comment_num + 1 cycle end if row_num = row_num + 1 end do close ( unit = input_unit ) return end subroutine get_unit ( iunit ) !*****************************************************************************80 ! !! GET_UNIT returns a free FORTRAN unit number. ! ! Discussion: ! ! A "free" FORTRAN unit number is an integer between 1 and 99 which ! is not currently associated with an I/O device. A free FORTRAN unit ! number is needed in order to open a file with the OPEN command. ! ! If IUNIT = 0, then no free FORTRAN unit could be found, although ! all 99 units were checked (except for units 5, 6 and 9, which ! are commonly reserved for console I/O). ! ! Otherwise, IUNIT is an integer between 1 and 99, representing a ! free FORTRAN unit. Note that GET_UNIT assumes that units 5 and 6 ! are special, and will never return those values. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 September 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Output, integer IUNIT, the free unit number. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer i integer ios integer iunit logical lopen iunit = 0 do i = 1, 99 if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then inquire ( unit = i, opened = lopen, iostat = ios ) if ( ios == 0 ) then if ( .not. lopen ) then iunit = i return end if end if end if end do return end subroutine i4mat_data_read ( input_filename, m, n, table ) !*****************************************************************************80 ! !! I4MAT_DATA_READ reads data from an I4MAT file. ! ! Discussion: ! ! The file may contain more than N points, but this routine ! will return after reading N points. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 27 January 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Output, integer TABLE(M,N), the table data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer ierror character ( len = * ) input_filename integer input_status integer input_unit integer j character ( len = 255 ) line integer table(m,n) integer x(m) ierror = 0 call get_unit ( input_unit ) open ( unit = input_unit, file = input_filename, status = 'old', & iostat = input_status ) if ( input_status /= 0 ) then ierror = 1 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_DATA_READ - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the input file "' // & trim ( input_filename ) // '" on unit ', input_unit stop end if j = 0 do while ( j < n ) read ( input_unit, '(a)', iostat = input_status ) line if ( input_status /= 0 ) then ierror = 2 write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_DATA_READ - Fatal error!' write ( *, '(a)' ) ' Error while reading lines of data.' write ( *, '(a,i8)' ) ' Number of values expected per line M = ', m write ( *, '(a,i8)' ) ' Number of data lines read, J = ', j write ( *, '(a,i8)' ) ' Number of data lines needed, N = ', n stop end if if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then cycle end if call s_to_i4vec ( line, m, x, ierror ) if ( ierror /= 0 ) then cycle end if j = j + 1 table(1:m,j) = x(1:m) end do close ( unit = input_unit ) return end subroutine i4mat_header_read ( input_filename, m, n ) !*****************************************************************************80 ! !! I4MAT_HEADER_READ reads the header from an I4MAT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 04 June 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) INPUT_FILENAME, the name of the input file. ! ! Output, integer M, spatial dimension. ! ! Output, integer N, the number of points. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character ( len = * ) input_filename integer m integer n call file_column_count ( input_filename, m ) if ( m <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data columns in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if call file_row_count ( input_filename, n ) if ( n <= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_HEADER_READ - Fatal error!' write ( *, '(a)' ) ' There was some kind of I/O problem while trying' write ( *, '(a)' ) ' to count the number of data rows in' write ( *, '(a)' ) ' the file "' // trim ( input_filename ) // '".' stop end if return end subroutine i4mat_transpose_print ( m, n, a, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT prints an I4MAT, transposed. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 December 2004 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer a(m,n) character ( len = * ) title call i4mat_transpose_print_some ( m, n, a, 1, 1, m, n, title ) return end subroutine i4mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title ) !*****************************************************************************80 ! !! I4MAT_TRANSPOSE_PRINT_SOME prints some of the transpose of an I4MAT. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 09 February 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer M, N, the number of rows and columns. ! ! Input, integer A(M,N), an M by N matrix to be printed. ! ! Input, integer ILO, JLO, the first row and column to print. ! ! Input, integer IHI, JHI, the last row and column to print. ! ! Input, character ( len = * ) TITLE, an optional title. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer, parameter :: incx = 10 integer m integer n integer a(m,n) character ( len = 7 ) ctemp(incx) integer i integer i2 integer i2hi integer i2lo integer ihi integer ilo integer inc integer j integer j2hi integer j2lo integer jhi integer jlo character ( len = * ) title if ( 0 < len_trim ( title ) ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) trim ( title ) end if do i2lo = max ( ilo, 1 ), min ( ihi, m ), incx i2hi = i2lo + incx - 1 i2hi = min ( i2hi, m ) i2hi = min ( i2hi, ihi ) inc = i2hi + 1 - i2lo write ( *, '(a)' ) ' ' do i = i2lo, i2hi i2 = i + 1 - i2lo write ( ctemp(i2), '(i7)') i end do write ( *, '('' Row '',10a7)' ) ctemp(1:inc) write ( *, '(a)' ) ' Col' write ( *, '(a)' ) ' ' j2lo = max ( jlo, 1 ) j2hi = min ( jhi, n ) do j = j2lo, j2hi do i2 = 1, inc i = i2lo - 1 + i2 write ( ctemp(i2), '(i7)' ) a(i,j) end do write ( *, '(i5,1x,10a7)' ) j, ( ctemp(i), i = 1, inc ) end do end do return end subroutine i4mat_write ( output_filename, m, n, table ) !*****************************************************************************80 ! !! I4MAT_WRITE writes an I4MAT file. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 31 May 2009 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) OUTPUT_FILENAME, the output file name. ! ! Input, integer M, the spatial dimension. ! ! Input, integer N, the number of points. ! ! Input, integer TABLE(M,N), the table data. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer m integer n integer j character ( len = * ) output_filename integer output_status integer output_unit character ( len = 30 ) string integer table(m,n) ! ! Open the file. ! call get_unit ( output_unit ) open ( unit = output_unit, file = output_filename, & status = 'replace', iostat = output_status ) if ( output_status /= 0 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'I4MAT_WRITE - Fatal error!' write ( *, '(a,i8)' ) ' Could not open the output file "' // & trim ( output_filename ) // '" on unit ', output_unit output_unit = -1 stop end if ! ! Create a format string. ! write ( string, '(a1,i8,a4)' ) '(', m, 'i10)' ! ! Write the data. ! do j = 1, n write ( output_unit, string ) table(1:m,j) end do ! ! Close the file. ! close ( unit = output_unit ) return end subroutine mesh_base_one ( node_num, element_order, element_num, element_node ) !*****************************************************************************80 ! !! 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. ! Since a FORTRAN90 program will naturally assume a 1-based indexing, it is ! necessary to check a given element definition and, if it is actually ! 0-based, to convert it. ! ! This function attempts to detect 9-based node indexing and correct it. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 02 October 2009 ! ! 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. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer element_num integer element_order integer element_node(element_order,element_num) integer node_max integer node_min integer node_num node_min = node_num + 1 node_max = -1 node_min = minval ( element_node(1:element_order,1:element_num) ) node_max = maxval ( element_node(1:element_order,1:element_num) ) if ( node_min == 0 .and. node_max == node_num - 1 ) then write ( *, '(a)' ) ' ' write ( *, '(a)' )'MESH_BASE_ONE:' write ( *, '(a)' )' The element indexing appears to be 0-based!' write ( *, '(a)' )' This will be converted to 1-based.' element_node(1:element_order,1:element_num) = & element_node(1:element_order,1:element_num) + 1 else if ( node_min == 1 .and. node_max == node_num ) then write ( *, '(a)' ) ' ' write ( *, '(a)' )'MESH_BASE_ONE:' write ( *, '(a)' )' The element indexing appears to be 1-based!' write ( *, '(a)' )' No conversion is necessary.' else write ( *, '(a)' ) ' ' write ( *, '(a)' ) 'MESH_BASE_ONE - Warning!' write ( *, '(a)' ) ' The element indexing is not of a recognized type.' write ( *, '(a,i8)' ) ' NODE_MIN = ', node_min write ( *, '(a,i8)' ) ' NODE_MAX = ', node_max write ( *, '(a,i8)' ) ' NODE_NUM = ', node_num end if return end subroutine s_to_i4 ( s, ival, ierror, length ) !*****************************************************************************80 ! !! S_TO_I4 reads an I4 from a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 28 June 2000 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, a string to be examined. ! ! Output, integer IVAL, the integer value read from the string. ! If the string is blank, then IVAL will be returned 0. ! ! Output, integer IERROR, an error flag. ! 0, no error. ! 1, an error occurred. ! ! Output, integer LENGTH, the number of characters of S used. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character c integer i integer ierror integer isgn integer istate integer ival integer length character ( len = * ) s ierror = 0 istate = 0 isgn = 1 ival = 0 do i = 1, len_trim ( s ) c = s(i:i) ! ! Haven't read anything. ! if ( istate == 0 ) then if ( c == ' ' ) then else if ( c == '-' ) then istate = 1 isgn = -1 else if ( c == '+' ) then istate = 1 isgn = + 1 else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then istate = 2 ival = ichar ( c ) - ichar ( '0' ) else ierror = 1 return end if ! ! Have read the sign, expecting digits. ! else if ( istate == 1 ) then if ( c == ' ' ) then else if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then istate = 2 ival = ichar ( c ) - ichar ( '0' ) else ierror = 1 return end if ! ! Have read at least one digit, expecting more. ! else if ( istate == 2 ) then if ( lle ( '0', c ) .and. lle ( c, '9' ) ) then ival = 10 * ival + ichar ( c ) - ichar ( '0' ) else ival = isgn * ival length = i - 1 return end if end if end do ! ! If we read all the characters in the string, see if we're OK. ! if ( istate == 2 ) then ival = isgn * ival length = len_trim ( s ) else ierror = 1 length = 0 end if return end subroutine s_to_i4vec ( s, n, ivec, ierror ) !*****************************************************************************80 ! !! S_TO_I4VEC reads an I4VEC from a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 08 October 2003 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be read. ! ! Input, integer N, the number of values expected. ! ! Output, integer IVEC(N), the values read from the string. ! ! Output, integer IERROR, error flag. ! 0, no errors occurred. ! -K, could not read data for entries -K through N. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer n integer i integer ierror integer ilo integer ivec(n) integer length character ( len = * ) s i = 0 ilo = 1 do while ( i < n ) i = i + 1 call s_to_i4 ( s(ilo:), ivec(i), ierror, length ) if ( ierror /= 0 ) then ierror = -i exit end if ilo = ilo + length end do return end subroutine s_word_count ( s, nword ) !*****************************************************************************80 ! !! S_WORD_COUNT counts the number of "words" in a string. ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 14 April 1999 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, character ( len = * ) S, the string to be examined. ! ! Output, integer NWORD, the number of "words" in the string. ! Words are presumed to be separated by one or more blanks. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) logical blank integer i integer lens integer nword character ( len = * ) s nword = 0 lens = len ( s ) if ( lens <= 0 ) then return end if blank = .true. do i = 1, lens if ( s(i:i) == ' ' ) then blank = .true. else if ( blank ) then nword = nword + 1 blank = .false. end if end do return end subroutine timestamp ( ) !*****************************************************************************80 ! !! TIMESTAMP prints the current YMDHMS date as a time stamp. ! ! Example: ! ! 31 May 2001 9:45:54.872 AM ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 18 May 2013 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! None ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) character ( len = 8 ) ampm integer d integer h integer m integer mm character ( len = 9 ), parameter, dimension(12) :: month = (/ & 'January ', 'February ', 'March ', 'April ', & 'May ', 'June ', 'July ', 'August ', & 'September', 'October ', 'November ', 'December ' /) integer n integer s integer values(8) integer y call date_and_time ( values = values ) y = values(1) m = values(2) d = values(3) h = values(5) n = values(6) s = values(7) mm = values(8) if ( h < 12 ) then ampm = 'AM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Noon' else ampm = 'PM' end if else h = h - 12 if ( h < 12 ) then ampm = 'PM' else if ( h == 12 ) then if ( n == 0 .and. s == 0 ) then ampm = 'Midnight' else ampm = 'AM' end if end if end if write ( *, '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) & d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, trim ( ampm ) return end subroutine triangulation_order6_to_order3 ( tri_num1, triangle1, tri_num2, & triangle2 ) !*****************************************************************************80 ! !! TRIANGULATION_ORDER6_TO_ORDER3 linearizes a quadratic triangulation. ! ! Discussion: ! ! A quadratic triangulation is assumed to consist of 6-node triangles, ! as in the following: ! ! 11-12-13-14-15 ! |\ |\ | ! | \ | \ | ! 6 7 8 9 10 ! | \ | \ | ! | \| \| ! 1--2--3--4--5 ! ! This routine rearranges information so as to define the 3-node ! triangulation: ! ! 11-12-13-14-15 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 6--7--8--9-10 ! |\ |\ |\ |\ | ! | \| \| \| \| ! 1--2--3--4--5 ! ! Licensing: ! ! This code is distributed under the MIT license. ! ! Modified: ! ! 22 March 2005 ! ! Author: ! ! John Burkardt ! ! Parameters: ! ! Input, integer TRI_NUM1, the number of triangles in the ! quadratic triangulation. ! ! Input, integer TRIANGLE1(6,TRI_NUM1), the quadratic ! triangulation. ! ! Input, integer TRI_NUM2, the number of triangles in the linear ! triangulation. TRI_NUM2 = 4 * TRI_NUM1. ! ! Output, integer TRIANGLE2(3,TRI_NUM2), the linear ! triangulation. ! implicit none integer, parameter :: rk = kind ( 1.0D+00 ) integer tri_num1 integer tri_num2 integer n1 integer n2 integer n3 integer n4 integer n5 integer n6 integer tri1 integer tri2 integer triangle1(6,tri_num1) integer triangle2(3,tri_num2) tri2 = 0 do tri1 = 1, tri_num1 n1 = triangle1(1,tri1) n2 = triangle1(2,tri1) n3 = triangle1(3,tri1) n4 = triangle1(4,tri1) n5 = triangle1(5,tri1) n6 = triangle1(6,tri1) write ( *, * ) tri1, n1, n2, n3, n4, n5, n6 tri2 = tri2 + 1 triangle2(1:3,tri2) = (/ n1, n4, n6 /) tri2 = tri2 + 1 triangle2(1:3,tri2) = (/ n2, n5, n4 /) tri2 = tri2 + 1 triangle2(1:3,tri2) = (/ n3, n6, n5 /) tri2 = tri2 + 1 triangle2(1:3,tri2) = (/ n4, n5, n6 /) end do return end