From 7ee440c1b1f9e2fe0864d283bb166927e0499cb2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 27 Feb 2022 17:54:11 +0100 Subject: [PATCH] separating functionality for more flexibility --- src/IO.f90 | 2 +- src/VTK.f90 | 224 ++++++++++++++++++++----------- src/grid/discretization_grid.f90 | 7 +- 3 files changed, 151 insertions(+), 82 deletions(-) diff --git a/src/IO.f90 b/src/IO.f90 index 441a0bf3d..591a8c733 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -440,7 +440,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) case (155) msg = 'material index out of bounds' case (180) - msg = 'missing/invalid material definition via State Variable 2' + msg = 'missing/invalid material definition' case (190) msg = 'unknown element type:' case (191) diff --git a/src/VTK.f90 b/src/VTK.f90 index cf5a47a20..b3ea65c67 100644 --- a/src/VTK.f90 +++ b/src/VTK.f90 @@ -11,24 +11,117 @@ module VTK implicit none private - public :: VTK_readVTI + public :: & + VTK_readVTIdataset_int, & + VTK_readVTI_cellsSizeOrigin contains !-------------------------------------------------------------------------------------------------- -!> @brief Parse vtk image data (.vti) +!> @brief Read integer dataset from a VTK image data (*.vti) file. !> @details https://vtk.org/Wiki/VTK_XML_Formats !-------------------------------------------------------------------------------------------------- -subroutine VTK_readVTI(cells,geomSize,origin,material, & - fileContent) +function VTK_readVTIdataset_int(fileContent,label) result(dataset) + + character(len=*), intent(in) :: & + label, & + fileContent + integer, dimension(:), allocatable :: & + dataset + + character(len=:), allocatable :: dataType, headerType, base64_str + logical :: compressed + + + call VTK_readVTIdataset_raw(base64_str,dataType,headerType,compressed, & + fileContent,label) + dataset = as_Int(base64_str,headerType,compressed,dataType) + + if (.not. allocated(dataset)) call IO_error(error_ID = 844, ext_msg='dataset "'//label//'" not found') + +end function VTK_readVTIdataset_int + + +!-------------------------------------------------------------------------------------------------- +!> @brief Read dataset as raw data (base64 string) from a VTK image data (*.vti) file. +!> @details https://vtk.org/Wiki/VTK_XML_Formats +!-------------------------------------------------------------------------------------------------- +subroutine VTK_readVTIdataset_raw(base64_str,dataType,headerType,compressed, & + fileContent,label) + + character(len=*), intent(in) :: & + label, & + fileContent + character(len=:), allocatable, intent(out) :: dataType, headerType, base64_str + logical, intent(out) :: compressed + + logical :: inFile,inImage,gotCellData + integer(pI64) :: & + startPos, endPos, & + s + + + inFile = .false. + inImage = .false. + startPos = 1_pI64 + outer: 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),'',kind=pI64) == 0_pI64) + if (index(fileContent(startPos:endPos),' @brief Read cells, size, and origin of an VTK image data (*.vti) file. +!> @details https://vtk.org/Wiki/VTK_XML_Formats +!-------------------------------------------------------------------------------------------------- +subroutine VTK_readVTI_cellsSizeOrigin(cells,geomSize,origin, & + fileContent) integer, dimension(3), intent(out) :: & cells ! # of cells (across all processes!) real(pReal), dimension(3), intent(out) :: & - geomSize, & ! physical size (across all processes!) + geomSize, & ! size (across all processes!) origin ! origin (across all processes!) - integer, dimension(:), intent(out), allocatable :: & - material character(len=*), intent(in) :: & fileContent @@ -42,14 +135,10 @@ subroutine VTK_readVTI(cells,geomSize,origin,material, & cells = -1 geomSize = -1.0_pReal - inFile = .false. - inImage = .false. - gotCelldata = .false. - -!-------------------------------------------------------------------------------------------------- -! parse XML file + inFile = .false. + inImage = .false. startPos = 1_pI64 - do while (startPos < len(fileContent,kind=pI64)) + outer: 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 @@ -65,49 +154,24 @@ subroutine VTK_readVTI(cells,geomSize,origin,material, & if (index(fileContent(startPos:endPos),'',kind=pI64) == 0_pI64) - if (index(fileContent(startPos:endPos),' @brief determine size and origin from coordinates -!------------------------------------------------------------------------------------------------ +!-------------------------------------------------------------------------------------------------- +!> @brief Determine size and origin from coordinates. +!-------------------------------------------------------------------------------------------------- subroutine cellsSizeOrigin(c,s,o,header) integer, dimension(3), intent(out) :: c @@ -120,7 +184,7 @@ subroutine cellsSizeOrigin(c,s,o,header) 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 + 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') @@ -138,15 +202,15 @@ subroutine cellsSizeOrigin(c,s,o,header) end subroutine -!------------------------------------------------------------------------------------------------ -!> @brief Interpret Base64 string in vtk XML file as integer of default kind -!------------------------------------------------------------------------------------------------ +!-------------------------------------------------------------------------------------------------- +!> @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 + 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 @@ -167,15 +231,15 @@ function as_Int(base64_str,headerType,compressed,dataType) end function as_Int -!------------------------------------------------------------------------------------------------ -!> @brief Interpret Base64 string in vtk XML file as integer of pReal kind -!------------------------------------------------------------------------------------------------ +!-------------------------------------------------------------------------------------------------- +!> @brief Interpret Base64 string in vtk XML file as real of kind pReal. +!-------------------------------------------------------------------------------------------------- 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 + 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 @@ -196,14 +260,14 @@ function as_pReal(base64_str,headerType,compressed,dataType) end function as_pReal -!------------------------------------------------------------------------------------------------ -!> @brief Interpret Base64 string in vtk XML file as bytes -!------------------------------------------------------------------------------------------------ +!-------------------------------------------------------------------------------------------------- +!> @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 + 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 @@ -217,19 +281,19 @@ function asBytes(base64_str,headerType,compressed) result(bytes) end function asBytes -!------------------------------------------------------------------------------------------------ -!> @brief Interpret compressed Base64 string in vtk XML file as bytes +!-------------------------------------------------------------------------------------------------- +!> @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) + character(len=*), intent(in) :: base64_str, & ! base64 encoded string + headerType ! header type (UInt32 or Uint64) integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes_inflated @@ -265,15 +329,15 @@ function asBytes_compressed(base64_str,headerType) result(bytes) end function asBytes_compressed -!------------------------------------------------------------------------------------------------ -!> @brief Interprete uncompressed Base64 string in vtk XML file as bytes +!-------------------------------------------------------------------------------------------------- +!> @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) + character(len=*), intent(in) :: base64_str, & ! base64 encoded string + headerType ! header type (UInt32 or Uint64) integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes integer(pI64) :: s @@ -300,9 +364,9 @@ function asBytes_uncompressed(base64_str,headerType) result(bytes) end function asBytes_uncompressed -!------------------------------------------------------------------------------------------------ +!-------------------------------------------------------------------------------------------------- !> @brief Get XML string value for given key. -!------------------------------------------------------------------------------------------------ +!-------------------------------------------------------------------------------------------------- pure function getXMLValue(line,key) character(len=*), intent(in) :: line, key @@ -337,9 +401,9 @@ pure function getXMLValue(line,key) end function -!------------------------------------------------------------------------------------------------ +!-------------------------------------------------------------------------------------------------- !> @brief Check for supported file format variants. -!------------------------------------------------------------------------------------------------ +!-------------------------------------------------------------------------------------------------- pure function fileFormatOk(line) character(len=*),intent(in) :: line diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 4601d411b..d00b217b0 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -76,7 +76,12 @@ subroutine discretization_grid_init(restart) if (worldrank == 0) then fileContent = IO_read(interface_geomFile) - call VTK_readVTI(cells,geomSize,origin,materialAt_global,fileContent) + call VTK_readVTI_cellsSizeOrigin(cells,geomSize,origin,fileContent) + materialAt_global = VTK_readVTIdataset_int(fileContent,'material') + 1 + if (any(materialAt_global < 1)) & + call IO_error(180,ext_msg='material ID < 1') + if (size(materialAt_global) /= product(cells)) & + call IO_error(180,ext_msg='mismatch of material IDs and # cells') fname = interface_geomFile if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) call results_openJobFile(parallel=.false.)