keep functionality separated

This commit is contained in:
Martin Diehl 2022-02-27 15:43:34 +01:00
parent 6d50066a04
commit 1cae6c4533
2 changed files with 352 additions and 336 deletions

350
src/VTK.f90 Normal file
View File

@ -0,0 +1,350 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Read file data from the Visualization toolkit
!--------------------------------------------------------------------------------------------------
module VTK
use prec
use zlib
use base64
use IO
implicit none
private
public :: VTK_readVTI
contains
!--------------------------------------------------------------------------------------------------
!> @brief Parse vtk image data (.vti)
!> @details https://vtk.org/Wiki/VTK_XML_Formats
!--------------------------------------------------------------------------------------------------
subroutine VTK_readVTI(cells,geomSize,origin,material, &
fileContent)
integer, dimension(3), intent(out) :: &
cells ! # of cells (across all processes!)
real(pReal), dimension(3), intent(out) :: &
geomSize, & ! physical size (across all processes!)
origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: &
material
character(len=*), intent(in) :: &
fileContent
character(len=:), allocatable :: dataType, headerType
logical :: inFile,inImage,gotCellData,compressed
integer(pI64) :: &
startPos, endPos, &
s
cells = -1
geomSize = -1.0_pReal
inFile = .false.
inImage = .false.
gotCelldata = .false.
!--------------------------------------------------------------------------------------------------
! parse XML file
startPos = 1_pI64
do while (startPos < len(fileContent,kind=pI64))
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
if (endPos < startPos) endPos = len(fileContent,kind=pI64) ! end of file without new line
if (.not. inFile) then
if (index(fileContent(startPos:endPos),'<VTKFile',kind=pI64) /= 0_pI64) then
inFile = .true.
if (.not. fileFormatOk(fileContent(startPos:endPos))) call IO_error(error_ID = 844, ext_msg='file format')
headerType = merge('UInt64','UInt32',getXMLValue(fileContent(startPos:endPos),'header_type')=='UInt64')
compressed = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor'
end if
else
if (.not. inImage) then
if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
inImage = .true.
call cellsSizeOrigin(cells,geomSize,origin,fileContent(startPos:endPos))
end if
else
if (index(fileContent(startPos:endPos),'<CellData',kind=pI64) /= 0_pI64) then
gotCellData = .true.
do while (index(fileContent(startPos:endPos),'</CellData>',kind=pI64) == 0_pI64)
if (index(fileContent(startPos:endPos),'<DataArray',kind=pI64) /= 0_pI64 .and. &
getXMLValue(fileContent(startPos:endPos),'Name') == 'material' ) then
if (getXMLValue(fileContent(startPos:endPos),'format') /= 'binary') &
call IO_error(error_ID = 844, ext_msg='format (material)')
dataType = getXMLValue(fileContent(startPos:endPos),'type')
startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace)
material = as_Int(fileContent(s:endPos),headerType,compressed,dataType)
exit
end if
startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
end do
end if
end if
end if
if (gotCellData) exit
startPos = endPos + 2_pI64
end do
if (.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if (size(material) /= product(cells)) call IO_error(error_ID = 844, ext_msg='size(material)')
if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if (any(cells<1)) call IO_error(error_ID = 844, ext_msg='cells')
material = material + 1
if (any(material<1)) call IO_error(error_ID = 844, ext_msg='material ID < 0')
contains
!------------------------------------------------------------------------------------------------
!> @brief determine size and origin from coordinates
!------------------------------------------------------------------------------------------------
subroutine cellsSizeOrigin(c,s,o,header)
integer, dimension(3), intent(out) :: c
real(pReal), dimension(3), intent(out) :: s,o
character(len=*), intent(in) :: header
character(len=:), allocatable :: temp
real(pReal), dimension(3) :: delta
integer :: i
temp = getXMLValue(header,'Direction')
if (temp /= '1 0 0 0 1 0 0 0 1' .and. temp /= '') & ! https://discourse.vtk.org/t/vti-specification/6526
call IO_error(error_ID = 844, ext_msg = 'coordinate order')
temp = getXMLValue(header,'WholeExtent')
if (any([(IO_intValue(temp,IO_stringPos(temp),i),i=1,5,2)] /= 0)) &
call IO_error(error_ID = 844, ext_msg = 'coordinate start')
c = [(IO_intValue(temp,IO_stringPos(temp),i),i=2,6,2)]
temp = getXMLValue(header,'Spacing')
delta = [(IO_floatValue(temp,IO_stringPos(temp),i),i=1,3)]
s = delta * real(c,pReal)
temp = getXMLValue(header,'Origin')
o = [(IO_floatValue(temp,IO_stringPos(temp),i),i=1,3)]
end subroutine
!------------------------------------------------------------------------------------------------
!> @brief Interpret Base64 string in vtk XML file as integer of default kind
!------------------------------------------------------------------------------------------------
function as_Int(base64_str,headerType,compressed,dataType)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType, & ! header type (UInt32 or Uint64)
dataType ! data type (Int32, Int64, Float32, Float64)
logical, intent(in) :: compressed ! indicate whether data is zlib compressed
integer, dimension(:), allocatable :: as_Int
select case(dataType)
case('Int32')
as_Int = int(prec_bytesToC_INT32_T(asBytes(base64_str,headerType,compressed)))
case('Int64')
as_Int = int(prec_bytesToC_INT64_T(asBytes(base64_str,headerType,compressed)))
case('Float32')
as_Int = int(prec_bytesToC_FLOAT (asBytes(base64_str,headerType,compressed)))
case('Float64')
as_Int = int(prec_bytesToC_DOUBLE (asBytes(base64_str,headerType,compressed)))
case default
call IO_error(844,ext_msg='unknown data type: '//trim(dataType))
end select
end function as_Int
!------------------------------------------------------------------------------------------------
!> @brief Interpret Base64 string in vtk XML file as integer of pReal kind
!------------------------------------------------------------------------------------------------
function as_pReal(base64_str,headerType,compressed,dataType)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType, & ! header type (UInt32 or Uint64)
dataType ! data type (Int32, Int64, Float32, Float64)
logical, intent(in) :: compressed ! indicate whether data is zlib compressed
real(pReal), dimension(:), allocatable :: as_pReal
select case(dataType)
case('Int32')
as_pReal = real(prec_bytesToC_INT32_T(asBytes(base64_str,headerType,compressed)),pReal)
case('Int64')
as_pReal = real(prec_bytesToC_INT64_T(asBytes(base64_str,headerType,compressed)),pReal)
case('Float32')
as_pReal = real(prec_bytesToC_FLOAT (asBytes(base64_str,headerType,compressed)),pReal)
case('Float64')
as_pReal = real(prec_bytesToC_DOUBLE (asBytes(base64_str,headerType,compressed)),pReal)
case default
call IO_error(844,ext_msg='unknown data type: '//trim(dataType))
end select
end function as_pReal
!------------------------------------------------------------------------------------------------
!> @brief Interpret Base64 string in vtk XML file as bytes
!------------------------------------------------------------------------------------------------
function asBytes(base64_str,headerType,compressed) result(bytes)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType ! header type (UInt32 or Uint64)
logical, intent(in) :: compressed ! indicate whether data is zlib compressed
integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes
if (compressed) then
bytes = asBytes_compressed(base64_str,headerType)
else
bytes = asBytes_uncompressed(base64_str,headerType)
end if
end function asBytes
!------------------------------------------------------------------------------------------------
!> @brief Interpret compressed Base64 string in vtk XML file as bytes
!> @details A compressed Base64 string consists of a header block and a data block
! [#blocks/#u-size/#p-size/#c-size-1/#c-size-2/.../#c-size-#blocks][DATA-1/DATA-2...]
! #blocks = Number of blocks
! #u-size = Block size before compression
! #p-size = Size of last partial block (zero if it not needed)
! #c-size-i = Size in bytes of block i after compression
!------------------------------------------------------------------------------------------------
function asBytes_compressed(base64_str,headerType) result(bytes)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType ! header type (UInt32 or Uint64)
integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes, bytes_inflated
integer(pI64), dimension(:), allocatable :: temp, size_inflated, size_deflated
integer(pI64) :: headerLen, nBlock, b,s,e
if (headerType == 'UInt32') then
temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(4_pI64)))),pI64)
nBlock = int(temp(1),pI64)
headerLen = 4_pI64 * (3_pI64 + nBlock)
temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
else if (headerType == 'UInt64') then
temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(8_pI64)))),pI64)
nBlock = int(temp(1),pI64)
headerLen = 8_pI64 * (3_pI64 + nBlock)
temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
end if
allocate(size_inflated(nBlock),source=temp(2))
size_inflated(nBlock) = merge(temp(3),temp(2),temp(3)/=0_pI64)
size_deflated = temp(4:)
bytes_inflated = base64_to_bytes(base64_str(base64_nChar(headerLen)+1_pI64:))
allocate(bytes(sum(size_inflated)))
e = 0_pI64
do b = 1, nBlock
s = e + 1_pI64
e = s + size_deflated(b) - 1_pI64
bytes(sum(size_inflated(:b-1))+1_pI64:sum(size_inflated(:b))) = zlib_inflate(bytes_inflated(s:e),size_inflated(b))
end do
end function asBytes_compressed
!------------------------------------------------------------------------------------------------
!> @brief Interprete uncompressed Base64 string in vtk XML file as bytes
!> @details An uncompressed Base64 string consists of N headers blocks and a N data blocks
![#bytes-1/DATA-1][#bytes-2/DATA-2]...
!------------------------------------------------------------------------------------------------
function asBytes_uncompressed(base64_str,headerType) result(bytes)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType ! header type (UInt32 or Uint64)
integer(pI64) :: s
integer(pI64), dimension(1) :: nByte
integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes
allocate(bytes(0))
s=0_pI64
if (headerType == 'UInt32') then
do while(s+base64_nChar(4_pI64)<(len(base64_str,pI64)))
nByte = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64)))),pI64)
bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64+nByte(1))),5_pI64)]
s = s + base64_nChar(4_pI64+nByte(1))
end do
else if (headerType == 'UInt64') then
do while(s+base64_nChar(8_pI64)<(len(base64_str,pI64)))
nByte = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64)))),pI64)
bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64+nByte(1))),9_pI64)]
s = s + base64_nChar(8_pI64+nByte(1))
end do
end if
end function asBytes_uncompressed
!------------------------------------------------------------------------------------------------
!> @brief Get XML string value for given key
!------------------------------------------------------------------------------------------------
pure function getXMLValue(line,key)
character(len=*), intent(in) :: line, key
character(len=:), allocatable :: getXMLValue
integer :: s,e
#ifdef __INTEL_COMPILER
character :: q
#endif
s = index(line," "//key,back=.true.)
if (s==0) then
getXMLValue = ''
else
e = s + 1 + scan(line(s+1:),"'"//'"')
if (scan(line(s:e-2),'=') == 0) then
getXMLValue = ''
else
s = e
! https://community.intel.com/t5/Intel-Fortran-Compiler/ICE-for-merge-with-strings/m-p/1207204#M151657
#ifdef __INTEL_COMPILER
q = line(s-1:s-1)
e = s + index(line(s:),q) - 1
#else
e = s + index(line(s:),merge("'",'"',line(s-1:s-1)=="'")) - 1
#endif
getXMLValue = line(s:e-1)
end if
end if
end function
!------------------------------------------------------------------------------------------------
!> @brief check for supported file format
!------------------------------------------------------------------------------------------------
pure function fileFormatOk(line)
character(len=*),intent(in) :: line
logical :: fileFormatOk
fileFormatOk = getXMLValue(line,'type') == 'ImageData' .and. &
getXMLValue(line,'byte_order') == 'LittleEndian' .and. &
getXMLValue(line,'compressor') /= 'vtkLZ4DataCompressor' .and. &
getXMLValue(line,'compressor') /= 'vtkLZMADataCompressor'
end function fileFormatOk
end subroutine VTK_readVTI
end module VTK

View File

@ -14,8 +14,7 @@ module discretization_grid
use prec use prec
use parallelization use parallelization
use system_routines use system_routines
use base64 use VTK
use zlib
use DAMASK_interface use DAMASK_interface
use IO use IO
use config use config
@ -77,7 +76,7 @@ subroutine discretization_grid_init(restart)
if (worldrank == 0) then if (worldrank == 0) then
fileContent = IO_read(interface_geomFile) fileContent = IO_read(interface_geomFile)
call readVTI(cells,geomSize,origin,materialAt_global,fileContent) call VTK_readVTI(cells,geomSize,origin,materialAt_global,fileContent)
fname = interface_geomFile fname = interface_geomFile
if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
call results_openJobFile(parallel=.false.) call results_openJobFile(parallel=.false.)
@ -166,339 +165,6 @@ subroutine discretization_grid_init(restart)
end subroutine discretization_grid_init end subroutine discretization_grid_init
!--------------------------------------------------------------------------------------------------
!> @brief Parse vtk image data (.vti)
!> @details https://vtk.org/Wiki/VTK_XML_Formats
!--------------------------------------------------------------------------------------------------
subroutine readVTI(cells,geomSize,origin,material, &
fileContent)
integer, dimension(3), intent(out) :: &
cells ! cells (across all processes!)
real(pReal), dimension(3), intent(out) :: &
geomSize, & ! size (across all processes!)
origin ! origin (across all processes!)
integer, dimension(:), intent(out), allocatable :: &
material
character(len=*), intent(in) :: &
fileContent
character(len=:), allocatable :: dataType, headerType
logical :: inFile,inImage,gotCellData,compressed
integer(pI64) :: &
startPos, endPos, &
s
cells = -1
geomSize = -1.0_pReal
inFile = .false.
inImage = .false.
gotCelldata = .false.
!--------------------------------------------------------------------------------------------------
! parse XML file
startPos = 1_pI64
do while (startPos < len(fileContent,kind=pI64))
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
if (endPos < startPos) endPos = len(fileContent,kind=pI64) ! end of file without new line
if (.not. inFile) then
if (index(fileContent(startPos:endPos),'<VTKFile',kind=pI64) /= 0_pI64) then
inFile = .true.
if (.not. fileFormatOk(fileContent(startPos:endPos))) call IO_error(error_ID = 844, ext_msg='file format')
headerType = merge('UInt64','UInt32',getXMLValue(fileContent(startPos:endPos),'header_type')=='UInt64')
compressed = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor'
end if
else
if (.not. inImage) then
if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
inImage = .true.
call cellsSizeOrigin(cells,geomSize,origin,fileContent(startPos:endPos))
end if
else
if (index(fileContent(startPos:endPos),'<CellData',kind=pI64) /= 0_pI64) then
gotCellData = .true.
do while (index(fileContent(startPos:endPos),'</CellData>',kind=pI64) == 0_pI64)
if (index(fileContent(startPos:endPos),'<DataArray',kind=pI64) /= 0_pI64 .and. &
getXMLValue(fileContent(startPos:endPos),'Name') == 'material' ) then
if (getXMLValue(fileContent(startPos:endPos),'format') /= 'binary') &
call IO_error(error_ID = 844, ext_msg='format (material)')
dataType = getXMLValue(fileContent(startPos:endPos),'type')
startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace)
material = as_Int(fileContent(s:endPos),headerType,compressed,dataType)
exit
end if
startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
end do
end if
end if
end if
if (gotCellData) exit
startPos = endPos + 2_pI64
end do
if (.not. allocated(material)) call IO_error(error_ID = 844, ext_msg='material data not found')
if (size(material) /= product(cells)) call IO_error(error_ID = 844, ext_msg='size(material)')
if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size')
if (any(cells<1)) call IO_error(error_ID = 844, ext_msg='cells')
material = material + 1
if (any(material<1)) call IO_error(error_ID = 844, ext_msg='material ID < 0')
contains
!------------------------------------------------------------------------------------------------
!> @brief determine size and origin from coordinates
!------------------------------------------------------------------------------------------------
subroutine cellsSizeOrigin(c,s,o,header)
integer, dimension(3), intent(out) :: c
real(pReal), dimension(3), intent(out) :: s,o
character(len=*), intent(in) :: header
character(len=:), allocatable :: temp
real(pReal), dimension(:), allocatable :: delta
integer :: i
temp = getXMLValue(header,'Direction')
if (temp /= '1 0 0 0 1 0 0 0 1' .and. temp /= '') & ! https://discourse.vtk.org/t/vti-specification/6526
call IO_error(error_ID = 844, ext_msg = 'coordinate order')
temp = getXMLValue(header,'WholeExtent')
if (any([(IO_intValue(temp,IO_stringPos(temp),i),i=1,5,2)] /= 0)) &
call IO_error(error_ID = 844, ext_msg = 'coordinate start')
c = [(IO_intValue(temp,IO_stringPos(temp),i),i=2,6,2)]
temp = getXMLValue(header,'Spacing')
delta = [(IO_floatValue(temp,IO_stringPos(temp),i),i=1,3)]
s = delta * real(c,pReal)
temp = getXMLValue(header,'Origin')
o = [(IO_floatValue(temp,IO_stringPos(temp),i),i=1,3)]
end subroutine
!------------------------------------------------------------------------------------------------
!> @brief Interpret Base64 string in vtk XML file as integer of default kind
!------------------------------------------------------------------------------------------------
function as_Int(base64_str,headerType,compressed,dataType)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType, & ! header type (UInt32 or Uint64)
dataType ! data type (Int32, Int64, Float32, Float64)
logical, intent(in) :: compressed ! indicate whether data is zlib compressed
integer, dimension(:), allocatable :: as_Int
select case(dataType)
case('Int32')
as_Int = int(prec_bytesToC_INT32_T(asBytes(base64_str,headerType,compressed)))
case('Int64')
as_Int = int(prec_bytesToC_INT64_T(asBytes(base64_str,headerType,compressed)))
case('Float32')
as_Int = int(prec_bytesToC_FLOAT (asBytes(base64_str,headerType,compressed)))
case('Float64')
as_Int = int(prec_bytesToC_DOUBLE (asBytes(base64_str,headerType,compressed)))
case default
call IO_error(844,ext_msg='unknown data type: '//trim(dataType))
end select
end function as_Int
!------------------------------------------------------------------------------------------------
!> @brief Interpret Base64 string in vtk XML file as integer of pReal kind
!------------------------------------------------------------------------------------------------
function as_pReal(base64_str,headerType,compressed,dataType)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType, & ! header type (UInt32 or Uint64)
dataType ! data type (Int32, Int64, Float32, Float64)
logical, intent(in) :: compressed ! indicate whether data is zlib compressed
real(pReal), dimension(:), allocatable :: as_pReal
select case(dataType)
case('Int32')
as_pReal = real(prec_bytesToC_INT32_T(asBytes(base64_str,headerType,compressed)),pReal)
case('Int64')
as_pReal = real(prec_bytesToC_INT64_T(asBytes(base64_str,headerType,compressed)),pReal)
case('Float32')
as_pReal = real(prec_bytesToC_FLOAT (asBytes(base64_str,headerType,compressed)),pReal)
case('Float64')
as_pReal = real(prec_bytesToC_DOUBLE (asBytes(base64_str,headerType,compressed)),pReal)
case default
call IO_error(844,ext_msg='unknown data type: '//trim(dataType))
end select
end function as_pReal
!------------------------------------------------------------------------------------------------
!> @brief Interpret Base64 string in vtk XML file as bytes
!------------------------------------------------------------------------------------------------
function asBytes(base64_str,headerType,compressed) result(bytes)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType ! header type (UInt32 or Uint64)
logical, intent(in) :: compressed ! indicate whether data is zlib compressed
integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes
if (compressed) then
bytes = asBytes_compressed(base64_str,headerType)
else
bytes = asBytes_uncompressed(base64_str,headerType)
end if
end function asBytes
!------------------------------------------------------------------------------------------------
!> @brief Interpret compressed Base64 string in vtk XML file as bytes
!> @details A compressed Base64 string consists of a header block and a data block
! [#blocks/#u-size/#p-size/#c-size-1/#c-size-2/.../#c-size-#blocks][DATA-1/DATA-2...]
! #blocks = Number of blocks
! #u-size = Block size before compression
! #p-size = Size of last partial block (zero if it not needed)
! #c-size-i = Size in bytes of block i after compression
!------------------------------------------------------------------------------------------------
function asBytes_compressed(base64_str,headerType) result(bytes)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType ! header type (UInt32 or Uint64)
integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes, bytes_inflated
integer(pI64), dimension(:), allocatable :: temp, size_inflated, size_deflated
integer(pI64) :: headerLen, nBlock, b,s,e
if (headerType == 'UInt32') then
temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(4_pI64)))),pI64)
nBlock = int(temp(1),pI64)
headerLen = 4_pI64 * (3_pI64 + nBlock)
temp = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
else if (headerType == 'UInt64') then
temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(8_pI64)))),pI64)
nBlock = int(temp(1),pI64)
headerLen = 8_pI64 * (3_pI64 + nBlock)
temp = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64)
end if
allocate(size_inflated(nBlock),source=temp(2))
size_inflated(nBlock) = merge(temp(3),temp(2),temp(3)/=0_pI64)
size_deflated = temp(4:)
bytes_inflated = base64_to_bytes(base64_str(base64_nChar(headerLen)+1_pI64:))
allocate(bytes(sum(size_inflated)))
e = 0_pI64
do b = 1, nBlock
s = e + 1_pI64
e = s + size_deflated(b) - 1_pI64
bytes(sum(size_inflated(:b-1))+1_pI64:sum(size_inflated(:b))) = zlib_inflate(bytes_inflated(s:e),size_inflated(b))
end do
end function asBytes_compressed
!------------------------------------------------------------------------------------------------
!> @brief Interprete uncompressed Base64 string in vtk XML file as bytes
!> @details An uncompressed Base64 string consists of N headers blocks and a N data blocks
![#bytes-1/DATA-1][#bytes-2/DATA-2]...
!------------------------------------------------------------------------------------------------
function asBytes_uncompressed(base64_str,headerType) result(bytes)
character(len=*), intent(in) :: base64_str, & ! base64 encoded string
headerType ! header type (UInt32 or Uint64)
integer(pI64) :: s
integer(pI64), dimension(1) :: nByte
integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes
allocate(bytes(0))
s=0_pI64
if (headerType == 'UInt32') then
do while(s+base64_nChar(4_pI64)<(len(base64_str,pI64)))
nByte = int(prec_bytesToC_INT32_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64)))),pI64)
bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(4_pI64+nByte(1))),5_pI64)]
s = s + base64_nChar(4_pI64+nByte(1))
end do
else if (headerType == 'UInt64') then
do while(s+base64_nChar(8_pI64)<(len(base64_str,pI64)))
nByte = int(prec_bytesToC_INT64_T(base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64)))),pI64)
bytes = [bytes,base64_to_bytes(base64_str(s+1_pI64:s+base64_nChar(8_pI64+nByte(1))),9_pI64)]
s = s + base64_nChar(8_pI64+nByte(1))
end do
end if
end function asBytes_uncompressed
!------------------------------------------------------------------------------------------------
!> @brief Get XML string value for given key
!------------------------------------------------------------------------------------------------
pure function getXMLValue(line,key)
character(len=*), intent(in) :: line, key
character(len=:), allocatable :: getXMLValue
integer :: s,e
#ifdef __INTEL_COMPILER
character :: q
#endif
s = index(line," "//key,back=.true.)
if (s==0) then
getXMLValue = ''
else
e = s + 1 + scan(line(s+1:),"'"//'"')
if (scan(line(s:e-2),'=') == 0) then
getXMLValue = ''
else
s = e
! https://community.intel.com/t5/Intel-Fortran-Compiler/ICE-for-merge-with-strings/m-p/1207204#M151657
#ifdef __INTEL_COMPILER
q = line(s-1:s-1)
e = s + index(line(s:),q) - 1
#else
e = s + index(line(s:),merge("'",'"',line(s-1:s-1)=="'")) - 1
#endif
getXMLValue = line(s:e-1)
end if
end if
end function
!------------------------------------------------------------------------------------------------
!> @brief check for supported file format
!------------------------------------------------------------------------------------------------
pure function fileFormatOk(line)
character(len=*),intent(in) :: line
logical :: fileFormatOk
fileFormatOk = getXMLValue(line,'type') == 'ImageData' .and. &
getXMLValue(line,'byte_order') == 'LittleEndian' .and. &
getXMLValue(line,'compressor') /= 'vtkLZ4DataCompressor' .and. &
getXMLValue(line,'compressor') /= 'vtkLZMADataCompressor'
end function fileFormatOk
end subroutine readVTI
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
!> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element) !> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element)
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------