program main
!*****************************************************************************80
!
!! MAIN is the main program for FEM_TO_XML.
!
! Discussion:
!
! FEM_TO_XML converts a 1D, 2D or 3D mesh to DOLFIN 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 connectivities.
! * 'prefix'.xml will contain the DOLFIN XML version of the data.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 03 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-3642230981
!
implicit none
integer ( kind = 4 ) arg_num
integer ( kind = 4 ) element
character ( len = 255 ) :: element_filename = ''
integer ( kind = 4 ), allocatable, dimension ( :, : ) :: element_node
integer ( kind = 4 ) element_num
integer ( kind = 4 ) element_order
integer ( kind = 4 ) i
integer ( kind = 4 ) iarg
integer ( kind = 4 ) iargc
integer ( kind = 4 ) m
character ( len = 255 ) :: node_filename = ''
integer ( kind = 4 ) node_num
real ( kind = 8 ), allocatable, dimension ( :, : ) :: node_x
character ( len = 255 ) prefix
character ( len = 255 ) :: xml_filename = ''
call timestamp ( )
write ( *, '(a)' ) ''
write ( *, '(a)' ) 'FEM_TO_XML'
write ( *, '(a)' ) ' FORTRAN90 version'
write ( *, '(a)' ) ' Convert a 1D, 2D or 3D mesh from FEM to DOLFIN XML format.'
write ( *, '(a)' ) ''
write ( *, '(a)' ) ' Read "prefix"_nodes.txt, node coordinates.'
write ( *, '(a)' ) &
' Read "prefix"_elements.txt, the element node connectivity.'
write ( *, '(a)' ) ''
write ( *, '(a)' ) &
' Create "prefix".xml, a corresponding DOLFIN XML mesh 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)' ) 'FEM_TO_XML:'
write ( *, '(a)' ) ' Please enter the filename prefix.'
read ( *, '(a)' ) prefix
end if
!
! Create the filenames.
!
node_filename = trim ( prefix ) // '_nodes.txt'
element_filename = trim ( prefix ) // '_elements.txt'
xml_filename = trim ( prefix ) // '.xml'
!
! Read the node data.
!
call r8mat_header_read ( node_filename, m, node_num )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Read the header of "' &
// trim ( node_filename ) //'".'
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Spatial dimension = ', m
write ( *, '(a,i8)' ) ' Number of points = ', node_num
allocate ( node_x(1:m,1:node_num) )
call r8mat_data_read ( node_filename, m, node_num, node_x )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Read the data in "' // trim ( node_filename ) //'".'
call r8mat_transpose_print_some ( m, node_num, node_x, 1, 1, &
m, 5, ' Initial portion of node coordinate data:' )
!
! Read the element data.
!
call i4mat_header_read ( element_filename, element_order, element_num )
if ( m == 1 ) then
if ( element_order == 2 ) then
else
write ( *, * ) ' '
write ( *, '(a)' ) 'FEM_TO_XML - Fatal error!'
write ( *, '(a)' ) ' 1D mesh elements must use 2 vertices.'
stop 1
end if
else if ( m == 2 ) then
if ( element_order == 3 ) then
elseif ( element_order == 6 ) then
else
write ( *, * ) ' '
write ( *, '(a)' ) 'FEM_TO_XML - Fatal error!'
write ( *, '(a)' ) ' 2D mesh elements must use 3 vertices.'
stop 1
end if
else if ( m == 3 ) then
if ( element_order == 4 ) then
else
write ( *, * ) ' '
write ( *, '(a)' ) 'FEM_TO_XML - Fatal error!'
write ( *, '(a)' ) ' 3D mesh elements must use 4 vertices.'
stop 1
end if
end if
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Read the header of "' &
// trim ( element_filename ) //'".'
write ( *, '(a)' ) ' '
write ( *, '(a,i8)' ) ' Element order = ', element_order
write ( *, '(a,i8)' ) ' Number of elements = ', element_num
allocate ( element_node(1:element_order,1:element_num) )
call i4mat_data_read ( element_filename, element_order, element_num, &
element_node )
write ( *, '(a)' ) ' '
write ( *, '(a)' ) ' Read the data in "' &
// trim ( element_filename ) //'".'
call 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 ) then
call xml_mesh1d_write ( xml_filename, m, node_num, node_x, &
element_order, element_num, element_node )
else if ( m == 2 ) then
call xml_mesh2d_write ( xml_filename, m, node_num, node_x, &
element_order, element_num, element_node )
else if ( m == 3 ) then
call xml_mesh3d_write ( xml_filename, m, node_num, node_x, &
element_order, element_num, element_node )
end if
write ( *, '(a)' ) ''
write ( *, '(a)' ) ' Created XML file "' // trim ( xml_filename ) // '".'
!
! Free memory.
!
deallocate ( element_node )
deallocate ( node_x )
!
! Terminate.
!
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'FEM_TO_XML'
write ( *, '(a)' ) ' Normal end of execution.'
write ( *, '(a)' ) ' '
call timestamp ( )
stop
end
subroutine ch_cap ( c )
!*****************************************************************************80
!
!! CH_CAP capitalizes a single character.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 19 July 1998
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input/output, character C, the character to capitalize.
!
implicit none
character c
integer ( kind = 4 ) 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 GNU LGPL 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
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 GNU LGPL license.
!
! Modified:
!
! 04 August 1999
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character C, the decimal digit, '0' through '9' or blank
! are legal.
!
! Output, integer ( kind = 4 ) DIGIT, the corresponding value. If C was
! 'illegal', then DIGIT is -1.
!
implicit none
character c
integer ( kind = 4 ) 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 GNU LGPL license.
!
! Modified:
!
! 21 June 2001
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) INPUT_FILENAME, the name of the file.
!
! Output, integer ( kind = 4 ) COLUMN_NUM, the number of columns in the file.
!
implicit none
integer ( kind = 4 ) column_num
logical got_one
character ( len = * ) input_filename
integer ( kind = 4 ) input_unit
integer ( kind = 4 ) ios
character ( len = 255 ) 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 )
stop 1
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 GNU LGPL license.
!
! Modified:
!
! 06 March 2003
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) INPUT_FILENAME, the name of the input file.
!
! Output, integer ( kind = 4 ) ROW_NUM, the number of rows found.
!
implicit none
integer ( kind = 4 ) bad_num
integer ( kind = 4 ) comment_num
integer ( kind = 4 ) ierror
character ( len = * ) input_filename
integer ( kind = 4 ) input_unit
integer ( kind = 4 ) ios
character ( len = 255 ) line
integer ( kind = 4 ) record_num
integer ( kind = 4 ) 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 1
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 GNU LGPL license.
!
! Modified:
!
! 18 September 2005
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Output, integer ( kind = 4 ) IUNIT, the free unit number.
!
implicit none
integer ( kind = 4 ) i
integer ( kind = 4 ) ios
integer ( kind = 4 ) 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 GNU LGPL license.
!
! Modified:
!
! 27 January 2005
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) INPUT_FILENAME, the name of the input file.
!
! Input, integer ( kind = 4 ) M, the spatial dimension.
!
! Input, integer ( kind = 4 ) N, the number of points.
!
! Output, integer ( kind = 4 ) TABLE(M,N), the table data.
!
implicit none
integer ( kind = 4 ) m
integer ( kind = 4 ) n
integer ( kind = 4 ) ierror
character ( len = * ) input_filename
integer ( kind = 4 ) input_status
integer ( kind = 4 ) input_unit
integer ( kind = 4 ) j
character ( len = 255 ) line
integer ( kind = 4 ) table(m,n)
integer ( kind = 4 ) 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 1
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 1
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 GNU LGPL license.
!
! Modified:
!
! 04 June 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) INPUT_FILENAME, the name of the input file.
!
! Output, integer ( kind = 4 ) M, spatial dimension.
!
! Output, integer ( kind = 4 ) N, the number of points.
!
implicit none
character ( len = * ) input_filename
integer ( kind = 4 ) m
integer ( kind = 4 ) 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 1
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 1
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 GNU LGPL license.
!
! Modified:
!
! 28 December 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer ( kind = 4 ) M, N, the number of rows and columns.
!
! Input, integer ( kind = 4 ) A(M,N), an M by N matrix to be printed.
!
! Input, character ( len = * ) TITLE, a title.
!
implicit none
integer ( kind = 4 ) m
integer ( kind = 4 ) n
integer ( kind = 4 ) 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 GNU LGPL license.
!
! Modified:
!
! 09 February 2005
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer ( kind = 4 ) M, N, the number of rows and columns.
!
! Input, integer ( kind = 4 ) A(M,N), an M by N matrix to be printed.
!
! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print.
!
! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print.
!
! Input, character ( len = * ) TITLE, a title.
!
implicit none
integer ( kind = 4 ), parameter :: incx = 10
integer ( kind = 4 ) m
integer ( kind = 4 ) n
integer ( kind = 4 ) a(m,n)
character ( len = 7 ) ctemp(incx)
integer ( kind = 4 ) i
integer ( kind = 4 ) i2
integer ( kind = 4 ) i2hi
integer ( kind = 4 ) i2lo
integer ( kind = 4 ) ihi
integer ( kind = 4 ) ilo
integer ( kind = 4 ) inc
integer ( kind = 4 ) j
integer ( kind = 4 ) j2hi
integer ( kind = 4 ) j2lo
integer ( kind = 4 ) jhi
integer ( kind = 4 ) jlo
character ( len = * ) title
write ( *, '(a)' ) ' '
write ( *, '(a)' ) trim ( title )
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,a1,10a7)' ) j, ':', ( ctemp(i), i = 1, inc )
end do
end do
return
end
subroutine mesh_base_zero ( node_num, element_order, element_num, element_node )
!*****************************************************************************80
!
!! MESH_BASE_ZERO ensures that the element definition is zero-based.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 02 October 2009
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes.
!
! Input, integer ( kind = 4 ) ELEMENT_ORDER, the order of the elements.
!
! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements.
!
! Input/output, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM),
! the element definitions.
!
implicit none
integer ( kind = 4 ) element_num
integer ( kind = 4 ) element_order
integer ( kind = 4 ) element
integer ( kind = 4 ) element_node(element_order,element_num)
integer ( kind = 4 ) node
integer ( kind = 4 ) node_max
integer ( kind = 4 ) node_min
integer ( kind = 4 ) node_num
integer ( kind = 4 ) order
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_ZERO:'
write ( *, '(a)' )' The element indexing appears to be 0-based!'
write ( *, '(a)' )' No conversion is necessary.'
else if ( node_min == 1 .and. node_max == node_num ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' )'MESH_BASE_ZERO:'
write ( *, '(a)' )' The element indexing appears to be 1-based!'
write ( *, '(a)' )' This will be converted to 0-based.'
element_node(1:element_order,1:element_num) = &
element_node(1:element_order,1:element_num) - 1
else
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'MESH_BASE_ZERO - 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 r8mat_data_read ( input_filename, m, n, table )
!*****************************************************************************80
!
!! R8MAT_DATA_READ reads data from an R8MAT file.
!
! Discussion:
!
! The file may contain more than N points, but this routine will
! return after reading N of them.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 18 October 2008
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) INPUT_FILENAME, the name of the input file.
!
! Input, integer ( kind = 4 ) M, the spatial dimension.
!
! Input, integer ( kind = 4 ) N, the number of points.
!
! Output, real ( kind = 8 ) TABLE(M,N), the table data.
!
implicit none
integer ( kind = 4 ) m
integer ( kind = 4 ) n
integer ( kind = 4 ) ierror
character ( len = * ) input_filename
integer ( kind = 4 ) input_status
integer ( kind = 4 ) input_unit
integer ( kind = 4 ) j
character ( len = 255 ) line
real ( kind = 8 ) table(m,n)
real ( kind = 8 ) 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)' ) 'R8MAT_DATA_READ - Fatal error!'
write ( *, '(a,i8)' ) ' Could not open the input file "' // &
trim ( input_filename ) // '" on unit ', input_unit
stop 1
end if
j = 0
do while ( j < n )
read ( input_unit, '(a)', iostat = input_status ) line
if ( input_status /= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'R8MAT_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 1
end if
if ( line(1:1) == '#' .or. len_trim ( line ) == 0 ) then
cycle
end if
call s_to_r8vec ( 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 r8mat_header_read ( input_filename, m, n )
!*****************************************************************************80
!
!! R8MAT_HEADER_READ reads the header from an R8MAT file.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 07 September 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) INPUT_FILENAME, the name of the input file.
!
! Output, integer ( kind = 4 ) M, spatial dimension.
!
! Output, integer ( kind = 4 ) N, the number of points.
!
implicit none
character ( len = * ) input_filename
integer ( kind = 4 ) m
integer ( kind = 4 ) n
call file_column_count ( input_filename, m )
if ( m <= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'R8MAT_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 1
end if
call file_row_count ( input_filename, n )
if ( n <= 0 ) then
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'R8MAT_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 1
end if
return
end
subroutine r8mat_transpose_print ( m, n, a, title )
!*****************************************************************************80
!
!! R8MAT_TRANSPOSE_PRINT prints an R8MAT, transposed.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 14 June 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer ( kind = 4 ) M, N, the number of rows and columns.
!
! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed.
!
! Input, character ( len = * ) TITLE, a title.
!
implicit none
integer ( kind = 4 ) m
integer ( kind = 4 ) n
real ( kind = 8 ) a(m,n)
character ( len = * ) title
call r8mat_transpose_print_some ( m, n, a, 1, 1, m, n, title )
return
end
subroutine r8mat_transpose_print_some ( m, n, a, ilo, jlo, ihi, jhi, title )
!*****************************************************************************80
!
!! R8MAT_TRANSPOSE_PRINT_SOME prints some of an R8MAT, transposed.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 14 June 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, integer ( kind = 4 ) M, N, the number of rows and columns.
!
! Input, real ( kind = 8 ) A(M,N), an M by N matrix to be printed.
!
! Input, integer ( kind = 4 ) ILO, JLO, the first row and column to print.
!
! Input, integer ( kind = 4 ) IHI, JHI, the last row and column to print.
!
! Input, character ( len = * ) TITLE, a title.
!
implicit none
integer ( kind = 4 ), parameter :: incx = 5
integer ( kind = 4 ) m
integer ( kind = 4 ) n
real ( kind = 8 ) a(m,n)
character ( len = 14 ) ctemp(incx)
integer ( kind = 4 ) i
integer ( kind = 4 ) i2
integer ( kind = 4 ) i2hi
integer ( kind = 4 ) i2lo
integer ( kind = 4 ) ihi
integer ( kind = 4 ) ilo
integer ( kind = 4 ) inc
integer ( kind = 4 ) j
integer ( kind = 4 ) j2hi
integer ( kind = 4 ) j2lo
integer ( kind = 4 ) jhi
integer ( kind = 4 ) jlo
character ( len = * ) title
write ( *, '(a)' ) ' '
write ( *, '(a)' ) trim ( title )
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,7x)') i
end do
write ( *, '('' Row '',5a14)' ) ctemp(1:inc)
write ( *, '(a)' ) ' Col'
j2lo = max ( jlo, 1 )
j2hi = min ( jhi, n )
do j = j2lo, j2hi
do i2 = 1, inc
i = i2lo - 1 + i2
write ( ctemp(i2), '(g14.6)' ) a(i,j)
end do
write ( *, '(i5,a1,5a14)' ) j, ':', ( ctemp(i), i = 1, inc )
end do
end do
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 GNU LGPL license.
!
! Modified:
!
! 28 June 2000
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) S, a string to be examined.
!
! Output, integer ( kind = 4 ) IVAL, the value read from the string.
! If the string is blank, then IVAL will be returned 0.
!
! Output, integer ( kind = 4 ) IERROR, an error flag.
! 0, no error.
! 1, an error occurred.
!
! Output, integer ( kind = 4 ) LENGTH, the number of characters of S
! used to make IVAL.
!
implicit none
character c
integer ( kind = 4 ) i
integer ( kind = 4 ) ierror
integer ( kind = 4 ) isgn
integer ( kind = 4 ) istate
integer ( kind = 4 ) ival
integer ( kind = 4 ) 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 GNU LGPL license.
!
! Modified:
!
! 08 October 2003
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) S, the string to be read.
!
! Input, integer ( kind = 4 ) N, the number of values expected.
!
! Output, integer ( kind = 4 ) IVEC(N), the values read from the string.
!
! Output, integer ( kind = 4 ) IERROR, error flag.
! 0, no errors occurred.
! -K, could not read data for entries -K through N.
!
implicit none
integer ( kind = 4 ) n
integer ( kind = 4 ) i
integer ( kind = 4 ) ierror
integer ( kind = 4 ) ilo
integer ( kind = 4 ) ivec(n)
integer ( kind = 4 ) 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_to_r8 ( s, dval, ierror, length )
!*****************************************************************************80
!
!! S_TO_R8 reads an R8 from a string.
!
! Discussion:
!
! The 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 number.
!
! Legal input is:
!
! 1 blanks,
! 2 '+' or '-' sign,
! 2.5 blanks
! 3 integer ( kind = 4 ) part,
! 4 decimal point,
! 5 fraction part,
! 6 'E' or 'e' or 'D' or 'd', exponent marker,
! 7 exponent sign,
! 8 exponent integer ( kind = 4 ) part,
! 9 exponent decimal point,
! 10 exponent fraction part,
! 11 blanks,
! 12 final comma or semicolon,
!
! with most quantities optional.
!
! Example:
!
! S DVAL
!
! '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 GNU LGPL license.
!
! Modified:
!
! 07 September 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) 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, real ( kind = 8 ) DVAL, the value read from the string.
!
! Output, integer ( kind = 4 ) IERROR, error flag.
! 0, no errors occurred.
! 1, 2, 6 or 7, the input number was garbled. The
! value of IERROR is the last type of input successfully
! read. For instance, 1 means initial blanks, 2 means
! a plus or minus sign, and so on.
!
! Output, integer ( kind = 4 ) LENGTH, the number of characters read
! to form the number, including any terminating
! characters such as a trailing comma or blanks.
!
implicit none
logical ch_eqi
character c
real ( kind = 8 ) dval
integer ( kind = 4 ) ierror
integer ( kind = 4 ) ihave
integer ( kind = 4 ) isgn
integer ( kind = 4 ) iterm
integer ( kind = 4 ) jbot
integer ( kind = 4 ) jsgn
integer ( kind = 4 ) jtop
integer ( kind = 4 ) length
integer ( kind = 4 ) nchar
integer ( kind = 4 ) ndig
real ( kind = 8 ) rbot
real ( kind = 8 ) rexp
real ( kind = 8 ) rtop
character ( len = * ) s
nchar = len_trim ( s )
ierror = 0
dval = 0.0D+00
length = -1
isgn = 1
rtop = 0
rbot = 1
jsgn = 1
jtop = 0
jbot = 1
ihave = 1
iterm = 0
do
length = length + 1
if ( nchar < length+1 ) then
exit
end if
c = s(length+1:length+1)
!
! Blank character.
!
if ( c == ' ' ) then
if ( ihave == 2 ) then
else if ( ihave == 6 .or. ihave == 7 ) then
iterm = 1
else if ( 1 < ihave ) then
ihave = 11
end if
!
! Comma.
!
else if ( c == ',' .or. c == ';' ) then
if ( ihave /= 1 ) then
iterm = 1
ihave = 12
length = length + 1
end if
!
! Minus sign.
!
else if ( c == '-' ) then
if ( ihave == 1 ) then
ihave = 2
isgn = -1
else if ( ihave == 6 ) then
ihave = 7
jsgn = -1
else
iterm = 1
end if
!
! Plus sign.
!
else if ( c == '+' ) then
if ( ihave == 1 ) then
ihave = 2
else if ( ihave == 6 ) then
ihave = 7
else
iterm = 1
end if
!
! Decimal point.
!
else if ( c == '.' ) then
if ( ihave < 4 ) then
ihave = 4
else if ( 6 <= ihave .and. ihave <= 8 ) then
ihave = 9
else
iterm = 1
end if
!
! Scientific notation exponent marker.
!
else if ( ch_eqi ( c, 'E' ) .or. ch_eqi ( c, 'D' ) ) then
if ( ihave < 6 ) then
ihave = 6
else
iterm = 1
end if
!
! Digit.
!
else if ( ihave < 11 .and. lle ( '0', c ) .and. lle ( c, '9' ) ) then
if ( ihave <= 2 ) then
ihave = 3
else if ( ihave == 4 ) then
ihave = 5
else if ( ihave == 6 .or. ihave == 7 ) then
ihave = 8
else if ( ihave == 9 ) then
ihave = 10
end if
call ch_to_digit ( c, ndig )
if ( ihave == 3 ) then
rtop = 10.0D+00 * rtop + real ( ndig, kind = 8 )
else if ( ihave == 5 ) then
rtop = 10.0D+00 * rtop + real ( ndig, kind = 8 )
rbot = 10.0D+00 * rbot
else if ( ihave == 8 ) then
jtop = 10 * jtop + ndig
else if ( ihave == 10 ) then
jtop = 10 * jtop + ndig
jbot = 10 * jbot
end if
!
! Anything else is regarded as a terminator.
!
else
iterm = 1
end if
!
! If we haven't seen a terminator, and we haven't examined the
! entire string, go get the next character.
!
if ( iterm == 1 ) then
exit
end if
end do
!
! If we haven't seen a terminator, and we have examined the
! entire string, then we're done, and LENGTH is equal to NCHAR.
!
if ( iterm /= 1 .and. length+1 == nchar ) then
length = nchar
end if
!
! 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 .or. ihave == 2 .or. ihave == 6 .or. ihave == 7 ) then
ierror = ihave
write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'S_TO_R8 - Serious error!'
write ( *, '(a)' ) ' Illegal or nonnumeric input:'
write ( *, '(a)' ) ' ' // trim ( s )
return
end if
!
! Number seems OK. Form it.
!
if ( jtop == 0 ) then
rexp = 1.0D+00
else
if ( jbot == 1 ) then
rexp = 10.0D+00 ** ( jsgn * jtop )
else
rexp = 10.0D+00 ** ( real ( jsgn * jtop, kind = 8 ) &
/ real ( jbot, kind = 8 ) )
end if
end if
dval = real ( isgn, kind = 8 ) * rexp * rtop / rbot
return
end
subroutine s_to_r8vec ( s, n, rvec, ierror )
!*****************************************************************************80
!
!! S_TO_R8VEC reads an R8VEC from a string.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 07 September 2004
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) S, the string to be read.
!
! Input, integer ( kind = 4 ) N, the number of values expected.
!
! Output, real ( kind = 8 ) RVEC(N), the values read from the string.
!
! Output, integer ( kind = 4 ) IERROR, error flag.
! 0, no errors occurred.
! -K, could not read data for entries -K through N.
!
implicit none
integer ( kind = 4 ) n
integer ( kind = 4 ) i
integer ( kind = 4 ) ierror
integer ( kind = 4 ) ilo
integer ( kind = 4 ) lchar
real ( kind = 8 ) rvec(n)
character ( len = * ) s
i = 0
ilo = 1
do while ( i < n )
i = i + 1
call s_to_r8 ( s(ilo:), rvec(i), ierror, lchar )
if ( ierror /= 0 ) then
ierror = -i
exit
end if
ilo = ilo + lchar
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 GNU LGPL license.
!
! Modified:
!
! 14 April 1999
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! Input, character ( len = * ) S, the string to be examined.
!
! Output, integer ( kind = 4 ) NWORD, the number of "words" in the string.
! Words are presumed to be separated by one or more blanks.
!
implicit none
logical blank
integer ( kind = 4 ) i
integer ( kind = 4 ) lens
integer ( kind = 4 ) 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 GNU LGPL license.
!
! Modified:
!
! 18 May 2013
!
! Author:
!
! John Burkardt
!
! Parameters:
!
! None
!
implicit none
character ( len = 8 ) ampm
integer ( kind = 4 ) d
integer ( kind = 4 ) h
integer ( kind = 4 ) m
integer ( kind = 4 ) mm
character ( len = 9 ), parameter, dimension(12) :: month = (/ &
'January ', 'February ', 'March ', 'April ', &
'May ', 'June ', 'July ', 'August ', &
'September', 'October ', 'November ', 'December ' /)
integer ( kind = 4 ) n
integer ( kind = 4 ) s
integer ( kind = 4 ) values(8)
integer ( kind = 4 ) 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 xml_mesh1d_write ( xml_filename, m, node_num, node_x, &
element_order, element_num, element_node )
!*****************************************************************************80
!
!! XML_MESH1D_WRITE writes 1D triangulation data as a DOLFIN XML mesh file.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 03 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-3642230981
!
! Parameters:
!
! Input, character * ( * ) XML_FILENAME, the name of the XML file
! to create.
!
! Input, integer ( kind = 4 ) M, the spatial dimension.
!
! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes.
!
! Input, real ( kind = 8 ) NODE_XY(M,NODE_NUM), the node coordinates.
!
! Input, integer ( kind = 4 ) ELEMENT_ORDER, the order of the elements.
!
! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements.
!
! Input, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM),
! the nodes that make up each element.
!
implicit none
integer ( kind = 4 ) element_num
integer ( kind = 4 ) element_order
integer ( kind = 4 ) m
integer ( kind = 4 ) node_num
integer ( kind = 4 ) element
integer ( kind = 4 ) element_node(element_order,element_num)
integer ( kind = 4 ) element_type
character * ( * ) xml_filename
integer ( kind = 4 ) xml_unit
integer ( kind = 4 ) node
real ( kind = 8 ) node_x(m,node_num)
character ( len = 80 ) s1
character ( len = 80 ) s2
character ( len = 80 ) s3
!
! Force 0-based indexing.
!
call mesh_base_zero ( node_num, element_order, element_num, element_node )
!
! Open the file.
!
call get_unit ( xml_unit )
open ( unit = xml_unit, file = xml_filename, status = 'replace' )
!
! Write the data.
!
write ( xml_unit, '(a)' ) ''
write ( xml_unit, '(a)' ) ''
write ( xml_unit, '(a)' ) ''
write ( s1, '(i1) ' ) m
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
write ( s1, '(i6)' ) node_num
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
do node = 1, node_num
write ( s1, '(i6)' ) node - 1
s1 = adjustl ( s1 )
write ( s2, '(g14.6)' ) node_x(1,node)
s2 = adjustl ( s2 )
write ( xml_unit, '(a)' ) ' '
end do
write ( xml_unit, '(a)' ) ' '
write ( s1, '(i6)' ) element_num
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
do element = 1, element_num
write ( s1, '(i6)' ) element - 1
s1 = adjustl ( s1 )
write ( s2, '(i6)' ) element_node(1,element)
s2 = adjustl ( s2 )
write ( s3, '(i6)' ) element_node(2,element)
s3 = adjustl ( s3 )
write ( xml_unit, '(a)' ) ' '
end do
write ( xml_unit, '(a)' ) ' '
write ( xml_unit, '(a)' ) ' '
write ( xml_unit, '(a)' ) ''
close ( unit = xml_unit )
return
end
subroutine xml_mesh2d_write ( xml_filename, m, node_num, node_x, &
element_order, element_num, element_node )
!*****************************************************************************80
!
!! XML_MESH2D_WRITE writes 2D mesh data as a DOLFIN XML mesh file.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 16 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-3642230981
!
! Parameters:
!
! Input, character * ( * ) XML_FILENAME, the name of the XML file
! to create.
!
! Input, integer ( kind = 4 ) M, the spatial dimension.
!
! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes.
!
! Input, real ( kind = 8 ) NODE_XY(M,NODE_NUM), the node coordinates.
!
! Input, integer ( kind = 4 ) ELEMENT_ORDER, the order of the elements.
!
! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements.
!
! Input, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM),
! the nodes that make up each element.
!
implicit none
integer ( kind = 4 ) element_num
integer ( kind = 4 ) element_order
integer ( kind = 4 ) m
integer ( kind = 4 ) node_num
integer ( kind = 4 ) element
integer ( kind = 4 ) element_node(element_order,element_num)
integer ( kind = 4 ) element_type
character * ( * ) xml_filename
integer ( kind = 4 ) xml_unit
integer ( kind = 4 ) node
real ( kind = 8 ) node_x(m,node_num)
character ( len = 80 ) s1
character ( len = 80 ) s2
character ( len = 80 ) s3
character ( len = 80 ) s4
!
! Force 0-based indexing.
!
call mesh_base_zero ( node_num, element_order, element_num, element_node )
!
! Open the file.
!
call get_unit ( xml_unit )
open ( unit = xml_unit, file = xml_filename, status = 'replace' )
!
! Write the data.
!
write ( xml_unit, '(a)' ) ''
write ( xml_unit, '(a)' ) ''
write ( xml_unit, '(a)' ) ''
write ( s1, '(i1) ' ) m
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
write ( s1, '(i6)' ) node_num
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
do node = 1, node_num
write ( s1, '(i6)' ) node - 1
s1 = adjustl ( s1 )
write ( s2, '(g14.6)' ) node_x(1,node)
s2 = adjustl ( s2 )
write ( s3, '(g14.6)' ) node_x(2,node)
s3 = adjustl ( s3 )
write ( xml_unit, '(a)' ) ' '
end do
write ( xml_unit, '(a)' ) ' '
write ( s1, '(i6)' ) element_num
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
do element = 1, element_num
write ( s1, '(i6)' ) element - 1
s1 = adjustl ( s1 )
write ( s2, '(i6)' ) element_node(1,element)
s2 = adjustl ( s2 )
write ( s3, '(i6)' ) element_node(2,element)
s3 = adjustl ( s3 )
write ( s4, '(i6)' ) element_node(3,element)
s4 = adjustl ( s4 )
write ( xml_unit, '(a)' ) &
' '
end do
write ( xml_unit, '(a)' ) ' '
write ( xml_unit, '(a)' ) ' '
write ( xml_unit, '(a)' ) ''
close ( unit = xml_unit )
return
end
subroutine xml_mesh3d_write ( xml_filename, m, node_num, node_x, &
element_order, element_num, element_node )
!*****************************************************************************80
!
!! XML_MESH3D_WRITE writes 3D mesh data as a DOLFIN XML mesh file.
!
! Licensing:
!
! This code is distributed under the GNU LGPL license.
!
! Modified:
!
! 16 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-3642230981
!
! Parameters:
!
! Input, character * ( * ) XML_FILENAME, the name of the XML file
! to create.
!
! Input, integer ( kind = 4 ) M, the spatial dimension.
!
! Input, integer ( kind = 4 ) NODE_NUM, the number of nodes.
!
! Input, real ( kind = 8 ) NODE_XY(M,NODE_NUM), the node coordinates.
!
! Input, integer ( kind = 4 ) ELEMENT_ORDER, the order of the elements.
!
! Input, integer ( kind = 4 ) ELEMENT_NUM, the number of elements.
!
! Input, integer ( kind = 4 ) ELEMENT_NODE(ELEMENT_ORDER,ELEMENT_NUM),
! the nodes that make up each element.
!
implicit none
integer ( kind = 4 ) element_num
integer ( kind = 4 ) element_order
integer ( kind = 4 ) m
integer ( kind = 4 ) node_num
integer ( kind = 4 ) element
integer ( kind = 4 ) element_node(element_order,element_num)
integer ( kind = 4 ) element_type
character * ( * ) xml_filename
integer ( kind = 4 ) xml_unit
integer ( kind = 4 ) node
real ( kind = 8 ) node_x(m,node_num)
character ( len = 80 ) s1
character ( len = 80 ) s2
character ( len = 80 ) s3
character ( len = 80 ) s4
character ( len = 80 ) s5
!
! Force 0-based indexing.
!
call mesh_base_zero ( node_num, element_order, element_num, element_node )
!
! Open the file.
!
call get_unit ( xml_unit )
open ( unit = xml_unit, file = xml_filename, status = 'replace' )
!
! Write the data.
!
write ( xml_unit, '(a)' ) ''
write ( xml_unit, '(a)' ) ''
write ( xml_unit, '(a)' ) ''
write ( s1, '(i1) ' ) m
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
write ( s1, '(i6)' ) node_num
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
do node = 1, node_num
write ( s1, '(i6)' ) node - 1
s1 = adjustl ( s1 )
write ( s2, '(g14.6)' ) node_x(1,node)
s2 = adjustl ( s2 )
write ( s3, '(g14.6)' ) node_x(2,node)
s3 = adjustl ( s3 )
write ( s4, '(g14.6)' ) node_x(3,node)
s4 = adjustl ( s4 )
write ( xml_unit, '(a)' ) ' '
end do
write ( xml_unit, '(a)' ) ' '
write ( s1, '(i6)' ) element_num
s1 = adjustl ( s1 )
write ( xml_unit, '(a)' ) ' '
do element = 1, element_num
write ( s1, '(i6)' ) element - 1
s1 = adjustl ( s1 )
write ( s2, '(i6)' ) element_node(1,element)
s2 = adjustl ( s2 )
write ( s3, '(i6)' ) element_node(2,element)
s3 = adjustl ( s3 )
write ( s4, '(i6)' ) element_node(3,element)
s4 = adjustl ( s4 )
write ( s5, '(i6)' ) element_node(4,element)
s5 = adjustl ( s5 )
write ( xml_unit, '(a)' ) &
' '
end do
write ( xml_unit, '(a)' ) ' '
write ( xml_unit, '(a)' ) ' '
write ( xml_unit, '(a)' ) ''
close ( unit = xml_unit )
return
end