From 1cae6c45333821b6d15200b8edbd5d2b797d8a02 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 27 Feb 2022 15:43:34 +0100 Subject: [PATCH 01/11] keep functionality separated --- src/VTK.f90 | 350 +++++++++++++++++++++++++++++++ src/grid/discretization_grid.f90 | 338 +---------------------------- 2 files changed, 352 insertions(+), 336 deletions(-) create mode 100644 src/VTK.f90 diff --git a/src/VTK.f90 b/src/VTK.f90 new file mode 100644 index 000000000..98735e344 --- /dev/null +++ b/src/VTK.f90 @@ -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),'',kind=pI64) == 0_pI64) + if (index(fileContent(startPos:endPos),' @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 diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index ced5f8eff..4601d411b 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -14,8 +14,7 @@ module discretization_grid use prec use parallelization use system_routines - use base64 - use zlib + use VTK use DAMASK_interface use IO use config @@ -77,7 +76,7 @@ subroutine discretization_grid_init(restart) if (worldrank == 0) then 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 if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) call results_openJobFile(parallel=.false.) @@ -166,339 +165,6 @@ subroutine discretization_grid_init(restart) 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),'',kind=pI64) == 0_pI64) - if (index(fileContent(startPos:endPos),' @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) !--------------------------------------------------------------------------------------------------- From efa30d1f3ae95a227b3291ad93727af71c7d9b7a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 27 Feb 2022 17:27:47 +0100 Subject: [PATCH 02/11] functions don't need acces to variables of function better keep them out --- src/VTK.f90 | 492 ++++++++++++++++++++++++++-------------------------- 1 file changed, 249 insertions(+), 243 deletions(-) diff --git a/src/VTK.f90 b/src/VTK.f90 index 98735e344..cf5a47a20 100644 --- a/src/VTK.f90 +++ b/src/VTK.f90 @@ -102,249 +102,255 @@ subroutine VTK_readVTI(cells,geomSize,origin,material, & 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 + +!------------------------------------------------------------------------------------------------ +!> @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 + + integer(C_SIGNED_CHAR), dimension(:), allocatable :: 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(C_SIGNED_CHAR), dimension(:), allocatable :: bytes + + integer(pI64) :: s + integer(pI64), dimension(1) :: nByte + + + 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 variants. +!------------------------------------------------------------------------------------------------ +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 module VTK From 7ee440c1b1f9e2fe0864d283bb166927e0499cb2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 27 Feb 2022 17:54:11 +0100 Subject: [PATCH 03/11] 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.) From f4f656cec479eac7ca2fada34f89a5d00d28c01f Mon Sep 17 00:00:00 2001 From: Sharan Date: Tue, 1 Mar 2022 20:39:48 +0100 Subject: [PATCH 04/11] point to current development in PRIVATE --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 4c8116ba3..39d26425f 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 4c8116ba3b9e9fbb325a580705028e8310139117 +Subproject commit 39d26425fadc312706b48acf5f465ac44746146c From f7b18981c93ed7ff464380e9edb5c8eae8aa44ab Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 4 Mar 2022 20:17:08 +0100 Subject: [PATCH 05/11] shorter names (w/o loss of clarity) --- src/{VTK.f90 => VTI.f90} | 28 ++++++++++++++-------------- src/grid/discretization_grid.f90 | 6 +++--- 2 files changed, 17 insertions(+), 17 deletions(-) rename src/{VTK.f90 => VTI.f90} (96%) diff --git a/src/VTK.f90 b/src/VTI.f90 similarity index 96% rename from src/VTK.f90 rename to src/VTI.f90 index b3ea65c67..763add989 100644 --- a/src/VTK.f90 +++ b/src/VTI.f90 @@ -1,8 +1,8 @@ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, KU Leuven -!> @brief Read file data from the Visualization toolkit +!> @brief Read data from image files of the visualization toolkit. !-------------------------------------------------------------------------------------------------- -module VTK +module VTI use prec use zlib use base64 @@ -12,8 +12,8 @@ module VTK private public :: & - VTK_readVTIdataset_int, & - VTK_readVTI_cellsSizeOrigin + VTI_readDataset_int, & + VTI_readCellsSizeOrigin contains @@ -21,7 +21,7 @@ contains !> @brief Read integer dataset from a VTK image data (*.vti) file. !> @details https://vtk.org/Wiki/VTK_XML_Formats !-------------------------------------------------------------------------------------------------- -function VTK_readVTIdataset_int(fileContent,label) result(dataset) +function VTI_readDataset_int(fileContent,label) result(dataset) character(len=*), intent(in) :: & label, & @@ -33,20 +33,20 @@ function VTK_readVTIdataset_int(fileContent,label) result(dataset) logical :: compressed - call VTK_readVTIdataset_raw(base64_str,dataType,headerType,compressed, & - fileContent,label) + call VTI_readDataset_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 +end function VTI_readDataset_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, & +subroutine VTI_readDataset_raw(base64_str,dataType,headerType,compressed, & fileContent,label) character(len=*), intent(in) :: & @@ -107,15 +107,15 @@ subroutine VTK_readVTIdataset_raw(base64_str,dataType,headerType,compressed, & end do outer -end subroutine VTK_readVTIdataset_raw +end subroutine VTI_readDataset_raw !-------------------------------------------------------------------------------------------------- !> @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) +subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, & + fileContent) integer, dimension(3), intent(out) :: & cells ! # of cells (across all processes!) @@ -166,7 +166,7 @@ subroutine VTK_readVTI_cellsSizeOrigin(cells,geomSize,origin, & 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') -end subroutine VTK_readVTI_cellsSizeOrigin +end subroutine VTI_readCellsSizeOrigin !-------------------------------------------------------------------------------------------------- @@ -417,4 +417,4 @@ pure function fileFormatOk(line) end function fileFormatOk -end module VTK +end module VTI diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index d00b217b0..691578fd0 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -14,7 +14,7 @@ module discretization_grid use prec use parallelization use system_routines - use VTK + use VTI use DAMASK_interface use IO use config @@ -76,8 +76,8 @@ subroutine discretization_grid_init(restart) if (worldrank == 0) then fileContent = IO_read(interface_geomFile) - call VTK_readVTI_cellsSizeOrigin(cells,geomSize,origin,fileContent) - materialAt_global = VTK_readVTIdataset_int(fileContent,'material') + 1 + call VTI_readCellsSizeOrigin(cells,geomSize,origin,fileContent) + materialAt_global = VTI_readDataset_int(fileContent,'material') + 1 if (any(materialAt_global < 1)) & call IO_error(180,ext_msg='material ID < 1') if (size(materialAt_global) /= product(cells)) & From d8de8f8b39ba410d65adf2a2db5a9b8615a29cb0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Mar 2022 22:13:37 +0100 Subject: [PATCH 06/11] initial conditions are of type real --- src/VTI.f90 | 44 +++++++++++++++++++++++++++++++++++--------- 1 file changed, 35 insertions(+), 9 deletions(-) diff --git a/src/VTI.f90 b/src/VTI.f90 index 763add989..d8b8f7778 100644 --- a/src/VTI.f90 +++ b/src/VTI.f90 @@ -13,6 +13,7 @@ module VTI public :: & VTI_readDataset_int, & + VTI_readDataset_real, & VTI_readCellsSizeOrigin contains @@ -23,10 +24,10 @@ contains !-------------------------------------------------------------------------------------------------- function VTI_readDataset_int(fileContent,label) result(dataset) - character(len=*), intent(in) :: & + character(len=*), intent(in) :: & label, & fileContent - integer, dimension(:), allocatable :: & + integer, dimension(:), allocatable :: & dataset character(len=:), allocatable :: dataType, headerType, base64_str @@ -42,6 +43,31 @@ function VTI_readDataset_int(fileContent,label) result(dataset) end function VTI_readDataset_int +!-------------------------------------------------------------------------------------------------- +!> @brief Read real dataset from a VTK image data (*.vti) file. +!> @details https://vtk.org/Wiki/VTK_XML_Formats +!-------------------------------------------------------------------------------------------------- +function VTI_readDataset_real(fileContent,label) result(dataset) + + character(len=*), intent(in) :: & + label, & + fileContent + real(pReal), dimension(:), allocatable :: & + dataset + + character(len=:), allocatable :: dataType, headerType, base64_str + logical :: compressed + + + call VTI_readDataset_raw(base64_str,dataType,headerType,compressed, & + fileContent,label) + dataset = as_real(base64_str,headerType,compressed,dataType) + + if (.not. allocated(dataset)) call IO_error(error_ID = 844, ext_msg='dataset "'//label//'" not found') + +end function VTI_readDataset_real + + !-------------------------------------------------------------------------------------------------- !> @brief Read dataset as raw data (base64 string) from a VTK image data (*.vti) file. !> @details https://vtk.org/Wiki/VTK_XML_Formats @@ -234,30 +260,30 @@ end function as_Int !-------------------------------------------------------------------------------------------------- !> @brief Interpret Base64 string in vtk XML file as real of kind pReal. !-------------------------------------------------------------------------------------------------- -function as_pReal(base64_str,headerType,compressed,dataType) +function as_real(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 + real(pReal), dimension(:), allocatable :: as_real select case(dataType) case('Int32') - as_pReal = real(prec_bytesToC_INT32_T(asBytes(base64_str,headerType,compressed)),pReal) + as_real = 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) + as_real = 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) + as_real = real(prec_bytesToC_FLOAT (asBytes(base64_str,headerType,compressed)),pReal) case('Float64') - as_pReal = real(prec_bytesToC_DOUBLE (asBytes(base64_str,headerType,compressed)),pReal) + as_real = 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 +end function as_real !-------------------------------------------------------------------------------------------------- From f54849f495bf38ac95c61975bf670ce11e65acb0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 5 Mar 2022 22:50:33 +0100 Subject: [PATCH 07/11] enhance grid to store initial conditions --- python/damask/_grid.py | 17 +++++++++++------ python/damask/_vtk.py | 29 ++++++++++++++++++++--------- 2 files changed, 31 insertions(+), 15 deletions(-) diff --git a/python/damask/_grid.py b/python/damask/_grid.py index ddc717b08..442afc6ef 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -34,7 +34,8 @@ class Grid: material: np.ndarray, size: FloatSequence, origin: FloatSequence = np.zeros(3), - comments: Union[str, Sequence[str]] = None): + comments: Union[str, Sequence[str]] = None, + initial_conditions = None): """ New geometry definition for grid solvers. @@ -55,7 +56,7 @@ class Grid: self.size = size # type: ignore self.origin = origin # type: ignore self.comments = [] if comments is None else comments # type: ignore - + self.ic = initial_conditions if initial_conditions is not None else {} def __repr__(self) -> str: """Give short human-readable summary.""" @@ -68,7 +69,7 @@ class Grid: f'origin: {util.srepr(self.origin," ")} m', f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else f' (min: {mat_min}, max: {mat_max})') - ]) + ]+(['initial_conditions:']+[f' - {f}' for f in self.ic.keys()] if self.ic else [])) def __copy__(self) -> 'Grid': @@ -187,11 +188,13 @@ class Grid: cells = np.array(v.vtk_data.GetDimensions())-1 bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T comments = v.comments + ic = {label:v.get(label).reshape(cells,order='F') for label in set(v.labels['Cell Data']) - {'material'}} return Grid(material = v.get('material').reshape(cells,order='F'), size = bbox[1] - bbox[0], origin = bbox[0], - comments = comments) + comments = comments, + initial_conditions = ic) @typing. no_type_check @@ -646,7 +649,9 @@ class Grid: """ v = VTK.from_image_data(self.cells,self.size,self.origin)\ .add(self.material.flatten(order='F'),'material') - v.comments += self.comments + for label,field in self.fields.items(): + v = v.add(field.flatten(order='F'),label) + v.comments = self.comments v.save(fname,parallel=False,compress=compress) @@ -683,7 +688,7 @@ class Grid: def show(self, - colormap: Colormap = Colormap.from_predefined('cividis')) -> None: + colormap: Colormap = None) -> None: """ Show on screen. diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 959cfe8ab..dcf0ca787 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -503,11 +503,21 @@ class VTK: def show(self, label: str = None, - colormap: Colormap = Colormap.from_predefined('cividis')): + colormap: Colormap = None): """ Render. - See http://compilatrix.com/article/vtk-1 for further ideas. + Parameters + ---------- + label : str, optional + Label of the dataset to show. + colormap : damask.Colormap, optional + Colormap for visualization of dataset. + + Notes + ----- + See http://compilatrix.com/article/vtk-1 for further ideas. + """ try: import wx @@ -525,8 +535,9 @@ class VTK: height = 768 lut = vtk.vtkLookupTable() - lut.SetNumberOfTableValues(len(colormap.colors)) - for i,c in enumerate(colormap.colors): + colormap_ = Colormap.from_predefined('cividis') if colormap is None else colormap + lut.SetNumberOfTableValues(len(colormap_.colors)) + for i,c in enumerate(colormap_.colors): lut.SetTableValue(i,c if len(c)==4 else np.append(c,1.0)) lut.Build() @@ -545,11 +556,11 @@ class VTK: if label is None: ren.SetBackground(67/255,128/255,208/255) else: - colormap = vtk.vtkScalarBarActor() - colormap.SetLookupTable(lut) - colormap.SetTitle(label) - colormap.SetMaximumWidthInPixels(width//100) - ren.AddActor2D(colormap) + colormap_vtk = vtk.vtkScalarBarActor() + colormap_vtk.SetLookupTable(lut) + colormap_vtk.SetTitle(label) + colormap_vtk.SetMaximumWidthInPixels(width//100) + ren.AddActor2D(colormap_vtk) ren.SetBackground(0.3,0.3,0.3) window = vtk.vtkRenderWindow() From bc020a9580ab25beab575a000eb33f8de425398d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 7 Mar 2022 10:38:56 +0100 Subject: [PATCH 08/11] correct names --- python/damask/_grid.py | 4 ++-- python/damask/_vtk.py | 8 ++++---- 2 files changed, 6 insertions(+), 6 deletions(-) diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 442afc6ef..0fd14dbf5 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -649,8 +649,8 @@ class Grid: """ v = VTK.from_image_data(self.cells,self.size,self.origin)\ .add(self.material.flatten(order='F'),'material') - for label,field in self.fields.items(): - v = v.add(field.flatten(order='F'),label) + for label,data in self.ic.items(): + v = v.add(data.flatten(order='F'),label) v.comments = self.comments v.save(fname,parallel=False,compress=compress) diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index dcf0ca787..a2a6bd4fa 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -88,10 +88,10 @@ class VTK: @property def comments(self) -> List[str]: """Return the comments.""" - fielddata = self.vtk_data.GetFieldData() - for a in range(fielddata.GetNumberOfArrays()): - if fielddata.GetArrayName(a) == 'comments': - comments = fielddata.GetAbstractArray(a) + field_data = self.vtk_data.GetFieldData() + for a in range(field_data.GetNumberOfArrays()): + if field_data.GetArrayName(a) == 'comments': + comments = field_data.GetAbstractArray(a) return [comments.GetValue(i) for i in range(comments.GetNumberOfValues())] return [] From fe8e55e4706b63abcdeb9252b1f1c2280a2fa33a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 7 Mar 2022 11:28:17 +0100 Subject: [PATCH 09/11] KeyError more sensible here --- python/damask/_vtk.py | 2 +- python/tests/test_Grid.py | 2 +- python/tests/test_VTK.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index a2a6bd4fa..90f5e2a1c 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -498,7 +498,7 @@ class VTK: # string array return np.array([vtk_array.GetValue(i) for i in range(vtk_array.GetNumberOfValues())]).astype(str) except UnboundLocalError: - raise ValueError(f'array "{label}" not found') + raise KeyError(f'array "{label}" not found') def show(self, diff --git a/python/tests/test_Grid.py b/python/tests/test_Grid.py index c453b6c74..4ef62ff93 100644 --- a/python/tests/test_Grid.py +++ b/python/tests/test_Grid.py @@ -61,7 +61,7 @@ class TestGrid: def test_invalid_no_material(self,tmp_path): v = VTK.from_image_data(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0) v.save(tmp_path/'no_materialpoint.vti',parallel=False) - with pytest.raises(ValueError): + with pytest.raises(KeyError): Grid.load(tmp_path/'no_materialpoint.vti') def test_invalid_material_type(self): diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index fc2ca179b..684a05678 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -141,7 +141,7 @@ class TestVTK: def test_invalid_get(self,default): - with pytest.raises(ValueError): + with pytest.raises(KeyError): default.get('does_not_exist') def test_invalid_add_shape(self,default): From 73f01c07d007deb11402cde9495fe4163b8fe4f7 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 8 Mar 2022 09:31:08 -0500 Subject: [PATCH 10/11] clarified colormap default; accept string as colormap name --- python/damask/_grid.py | 6 +++--- python/damask/_vtk.py | 10 ++++++---- python/tests/test_Grid.py | 6 ++++-- python/tests/test_VTK.py | 7 ++++--- src/grid/discretization_grid.f90 | 2 +- 5 files changed, 18 insertions(+), 13 deletions(-) diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 0fd14dbf5..0077a0309 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -688,14 +688,14 @@ class Grid: def show(self, - colormap: Colormap = None) -> None: + colormap: Union[Colormap, str] = None) -> None: """ Show on screen. Parameters ---------- - colormap : damask.Colormap, optional - Colors used to map material IDs. Defaults to 'cividis'. + colormap : damask.Colormap or str, optional + Colormap for visualization of material IDs. Defaults to 'cividis'. """ VTK.from_image_data(self.cells,self.size,self.origin) \ diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 90f5e2a1c..13b345091 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -503,7 +503,7 @@ class VTK: def show(self, label: str = None, - colormap: Colormap = None): + colormap: Union[Colormap, str] = None): """ Render. @@ -511,8 +511,8 @@ class VTK: ---------- label : str, optional Label of the dataset to show. - colormap : damask.Colormap, optional - Colormap for visualization of dataset. + colormap : damask.Colormap or str, optional + Colormap for visualization of dataset. Defaults to 'cividis'. Notes ----- @@ -535,7 +535,9 @@ class VTK: height = 768 lut = vtk.vtkLookupTable() - colormap_ = Colormap.from_predefined('cividis') if colormap is None else colormap + colormap_ = Colormap.from_predefined('cividis') if colormap is None \ + else Colormap.from_predefined(colormap) if isinstance(colormap,str) \ + else colormap lut.SetNumberOfTableValues(len(colormap_.colors)) for i,c in enumerate(colormap_.colors): lut.SetTableValue(i,c if len(c)==4 else np.append(c,1.0)) diff --git a/python/tests/test_Grid.py b/python/tests/test_Grid.py index 4ef62ff93..417d92507 100644 --- a/python/tests/test_Grid.py +++ b/python/tests/test_Grid.py @@ -6,6 +6,7 @@ from damask import VTK from damask import Grid from damask import Table from damask import Rotation +from damask import Colormap from damask import util from damask import seeds from damask import grid_filters @@ -42,9 +43,10 @@ class TestGrid: print('patched datetime.datetime.now') - def test_show(sef,default,monkeypatch): + @pytest.mark.parametrize('cmap',[Colormap.from_predefined('cividis'),'cividis',None]) + def test_show(sef,default,cmap,monkeypatch): monkeypatch.delenv('DISPLAY',raising=False) - default.show() + default.show(cmap) def test_equal(self,default): assert default == default diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index 684a05678..2cf7e411f 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -10,6 +10,7 @@ import vtk from damask import VTK from damask import Table +from damask import Colormap @pytest.fixture def ref_path(ref_path_base): @@ -29,10 +30,10 @@ class TestVTK: def _patch_execution_stamp(self, patch_execution_stamp): print('patched damask.util.execution_stamp') - def test_show(sef,default,monkeypatch): + @pytest.mark.parametrize('cmap',[Colormap.from_predefined('cividis'),'cividis',None]) + def test_show(sef,default,cmap,monkeypatch): monkeypatch.delenv('DISPLAY',raising=False) - default.show() - + default.show(colormap=cmap) def test_imageData(self,tmp_path): cells = np.random.randint(5,10,3) diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 691578fd0..399bb180a 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -81,7 +81,7 @@ subroutine discretization_grid_init(restart) 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') + call IO_error(180,ext_msg='mismatch in # of material IDs and cells') fname = interface_geomFile if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:) call results_openJobFile(parallel=.false.) From c4a7c0096a7ae69c343bf40d767226eb57a3ab55 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 8 Mar 2022 22:43:54 +0100 Subject: [PATCH 11/11] give default directly only Colormap object caused problem (tab completion triggered '__repr__' which means showing colormap in maplotlib window --- python/damask/_grid.py | 2 +- python/damask/_vtk.py | 9 ++++----- python/tests/test_Grid.py | 2 +- python/tests/test_VTK.py | 2 +- 4 files changed, 7 insertions(+), 8 deletions(-) diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 0077a0309..b20b52366 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -688,7 +688,7 @@ class Grid: def show(self, - colormap: Union[Colormap, str] = None) -> None: + colormap: Union[Colormap, str] = 'cividis') -> None: """ Show on screen. diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 13b345091..f8dce1754 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -503,7 +503,7 @@ class VTK: def show(self, label: str = None, - colormap: Union[Colormap, str] = None): + colormap: Union[Colormap, str] = 'cividis'): """ Render. @@ -516,7 +516,7 @@ class VTK: Notes ----- - See http://compilatrix.com/article/vtk-1 for further ideas. + See http://compilatrix.com/article/vtk-1 for further ideas. """ try: @@ -535,9 +535,8 @@ class VTK: height = 768 lut = vtk.vtkLookupTable() - colormap_ = Colormap.from_predefined('cividis') if colormap is None \ - else Colormap.from_predefined(colormap) if isinstance(colormap,str) \ - else colormap + colormap_ = Colormap.from_predefined(colormap) if isinstance(colormap,str) else \ + colormap lut.SetNumberOfTableValues(len(colormap_.colors)) for i,c in enumerate(colormap_.colors): lut.SetTableValue(i,c if len(c)==4 else np.append(c,1.0)) diff --git a/python/tests/test_Grid.py b/python/tests/test_Grid.py index 417d92507..5721c6428 100644 --- a/python/tests/test_Grid.py +++ b/python/tests/test_Grid.py @@ -43,7 +43,7 @@ class TestGrid: print('patched datetime.datetime.now') - @pytest.mark.parametrize('cmap',[Colormap.from_predefined('cividis'),'cividis',None]) + @pytest.mark.parametrize('cmap',[Colormap.from_predefined('stress'),'viridis']) def test_show(sef,default,cmap,monkeypatch): monkeypatch.delenv('DISPLAY',raising=False) default.show(cmap) diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index 2cf7e411f..c776135e1 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -30,7 +30,7 @@ class TestVTK: def _patch_execution_stamp(self, patch_execution_stamp): print('patched damask.util.execution_stamp') - @pytest.mark.parametrize('cmap',[Colormap.from_predefined('cividis'),'cividis',None]) + @pytest.mark.parametrize('cmap',[Colormap.from_predefined('cividis'),'strain']) def test_show(sef,default,cmap,monkeypatch): monkeypatch.delenv('DISPLAY',raising=False) default.show(colormap=cmap)