diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 9407ffc7c..ed5e762b7 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -492,7 +492,7 @@ mergeIntoMaster: removeData: stage: clean before_script: - - echo 'Do nothing' + - echo "Removing data and lock of pipeline $CI_PIPELINE_ID" script: - rm -rf $TESTROOT/GitLabCI_Pipeline_$CI_PIPELINE_ID - sed -i "/$CI_PIPELINE_ID/d" $TESTROOT/GitLabCI.queue # in case pipeline was manually (web GUI) restarted and releaseLock was performed already @@ -504,7 +504,7 @@ removeData: removeLock: stage: releaseLock before_script: - - echo 'Do nothing' + - echo "Removing lock of pipeline $CI_PIPELINE_ID" when: always script: sed -i "/$CI_PIPELINE_ID/d" $TESTROOT/GitLabCI.queue except: diff --git a/PRIVATE b/PRIVATE index 3bc603489..92ca3e83b 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 3bc60348981a68bc22a5ebaafe4173b7513ac264 +Subproject commit 92ca3e83b6093c1af277cfc06a504e4bb09fe8bc diff --git a/VERSION b/VERSION index 5dac27a44..2d2d2952b 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha-136-gc0ea8e3f +v3.0.0-alpha-147-gf0806a9e diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 7f3e64aef..d23bfd22e 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -168,7 +168,7 @@ class VTK: def _write(writer): """Wrapper for parallel writing.""" writer.Write() - def to_file(self,fname,parallel=True): + def to_file(self,fname,parallel=True,compress=True): """ Write to file. @@ -192,7 +192,10 @@ class VTK: if ext and ext != '.'+default_ext: raise ValueError(f'Given extension {ext} does not match default .{default_ext}') writer.SetFileName(str(Path(fname).with_suffix('.'+default_ext))) - writer.SetCompressorTypeToZLib() + if compress: + writer.SetCompressorTypeToZLib() + else: + writer.SetCompressorTypeToNone() writer.SetDataModeToBinary() writer.SetInputData(self.vtk_data) diff --git a/src/IO.f90 b/src/IO.f90 index 23f70eb8c..f434539b0 100644 --- a/src/IO.f90 +++ b/src/IO.f90 @@ -585,6 +585,8 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg) msg = 'incomplete information in grid mesh header' case (843) msg = 'microstructure count mismatch' + case (844) + msg = 'invalid VTR file' case (846) msg = 'rotation for load case rotation ill-defined (R:RT != I)' case (891) diff --git a/src/base64.f90 b/src/base64.f90 index eb2a1e929..3d7a51987 100644 --- a/src/base64.f90 +++ b/src/base64.f90 @@ -16,7 +16,7 @@ module base64 public :: & base64_init, & base64_to_bytes, & - base64_nBase64, & + base64_nChar, & base64_nByte contains @@ -37,14 +37,14 @@ end subroutine base64_init !-------------------------------------------------------------------------------------------------- !> @brief Calculate number of Base64 characters required for storage of N bytes. !-------------------------------------------------------------------------------------------------- -pure function base64_nBase64(nByte) +pure function base64_nChar(nByte) - integer(pLongInt), intent(in) :: nByte - integer(pLongInt) :: base64_nBase64 + integer(pI64), intent(in) :: nByte + integer(pI64) :: base64_nChar - base64_nBase64 = 4_pLongInt * (nByte/3_pLongInt + merge(1_pLongInt,0_pLongInt,mod(nByte,3_pLongInt) /= 0_pLongInt)) + base64_nChar = 4_pI64 * (nByte/3_pI64 + merge(1_pI64,0_pI64,mod(nByte,3_pI64) /= 0_pI64)) -end function base64_nBase64 +end function base64_nChar !-------------------------------------------------------------------------------------------------- @@ -52,10 +52,10 @@ end function base64_nBase64 !-------------------------------------------------------------------------------------------------- pure function base64_nByte(nBase64) - integer(pLongInt), intent(in) :: nBase64 - integer(pLongInt) :: base64_nByte + integer(pI64), intent(in) :: nBase64 + integer(pI64) :: base64_nByte - base64_nByte = 3_pLongInt * (nBase64/4_pLongInt) + base64_nByte = 3_pI64 * (nBase64/4_pI64) end function base64_nByte @@ -66,36 +66,36 @@ end function base64_nByte function base64_to_bytes(base64_str,s,e) result(bytes) character(len=*), intent(in) :: base64_str !< Base64 string representation - integer(pLongInt), intent(in), optional :: & + integer(pI64), intent(in), optional :: & s, & !< start (in bytes) e !< end (in bytes) - integer(pLongInt) :: s_bytes, e_bytes, s_str, e_str + integer(pI64) :: s_bytes, e_bytes, s_str, e_str integer(C_SIGNED_CHAR), dimension(:), allocatable :: bytes - if(.not. valid_base64(base64_str)) call IO_error(114,ext_msg='invalid character') + if(.not. validBase64(base64_str)) call IO_error(114,ext_msg='invalid character') if(present(s)) then - if(s<1_pLongInt) call IO_error(114, ext_msg='s out of range') - s_str = ((s-1_pLongInt)/3_pLongInt)*4_pLongInt + 1_pLongInt - s_bytes = mod(s-1_pLongInt,3_pLongInt) + 1_pLongInt + if(s<1_pI64) call IO_error(114, ext_msg='s out of range') + s_str = ((s-1_pI64)/3_pI64)*4_pI64 + 1_pI64 + s_bytes = mod(s-1_pI64,3_pI64) + 1_pI64 else - s_str = 1_pLongInt - s_bytes = 1_pLongInt + s_str = 1_pI64 + s_bytes = 1_pI64 endif if(present(e)) then - if(e>base64_nByte(len(base64_str,kind=pLongInt))) call IO_error(114, ext_msg='e out of range') - e_str = ((e-1_pLongInt)/3_pLongInt)*4_pLongInt + 4_pLongInt + if(e>base64_nByte(len(base64_str,kind=pI64))) call IO_error(114, ext_msg='e out of range') + e_str = ((e-1_pI64)/3_pI64)*4_pI64 + 4_pI64 e_bytes = e - base64_nByte(s_str) else - e_str = len(base64_str,kind=pLongInt) - e_bytes = base64_nByte(len(base64_str,kind=pLongInt)) - base64_nByte(s_str) - if(base64_str(e_str-0_pLongInt:e_str-0_pLongInt) == '=') e_bytes = e_bytes - 1_pLongInt - if(base64_str(e_str-1_pLongInt:e_str-1_pLongInt) == '=') e_bytes = e_bytes - 1_pLongInt + e_str = len(base64_str,kind=pI64) + e_bytes = base64_nByte(len(base64_str,kind=pI64)) - base64_nByte(s_str) + if(base64_str(e_str-0_pI64:e_str-0_pI64) == '=') e_bytes = e_bytes - 1_pI64 + if(base64_str(e_str-1_pI64:e_str-1_pI64) == '=') e_bytes = e_bytes - 1_pI64 endif - bytes = decode_base64(base64_str(s_str:e_str)) + bytes = decodeBase64(base64_str(s_str:e_str)) bytes = bytes(s_bytes:e_bytes) end function base64_to_bytes @@ -104,21 +104,21 @@ end function base64_to_bytes !-------------------------------------------------------------------------------------------------- !> @brief Convert a Base64 ASCII string into its byte-wise binary representation. !-------------------------------------------------------------------------------------------------- -pure function decode_base64(base64_str) result(bytes) +pure function decodeBase64(base64_str) result(bytes) character(len=*), intent(in) :: base64_str !< Base64 string representation - integer(C_SIGNED_CHAR), dimension(base64_nByte(len(base64_str,pLongInt))) :: bytes + integer(C_SIGNED_CHAR), dimension(base64_nByte(len(base64_str,pI64))) :: bytes integer(C_SIGNED_CHAR), dimension(0:3) :: charPos - integer(pLongInt) :: c, b, p + integer(pI64) :: c, b, p - c = 1_pLongInt - b = 1_pLongInt + c = 1_pI64 + b = 1_pI64 - do while(c < len(base64_str,kind=pLongInt)) - do p=0_pLongInt,3_pLongInt - if(c+p<=len(base64_str,kind=pLongInt)) then + do while(c < len(base64_str,kind=pI64)) + do p=0_pI64,3_pI64 + if(c+p<=len(base64_str,kind=pI64)) then charPos(p) = int(index(base64_encoding,base64_str(c+p:c+p))-1,C_SIGNED_CHAR) else charPos(p) = 0_C_SIGNED_CHAR @@ -131,31 +131,31 @@ pure function decode_base64(base64_str) result(bytes) call mvbits(charPos(2),2,4,bytes(b+1),0) call mvbits(charPos(2),0,2,bytes(b+2),6) call mvbits(charPos(3),0,6,bytes(b+2),0) - b = b+3_pLongInt - c = c+4_pLongInt + b = b+3_pI64 + c = c+4_pI64 enddo -end function decode_base64 +end function decodeBase64 !-------------------------------------------------------------------------------------------------- !> @brief Test for valid Base64 encoded string. !> @details Input string must be properly padded. !-------------------------------------------------------------------------------------------------- -pure logical function valid_base64(base64_str) +pure logical function validBase64(base64_str) character(len=*), intent(in) :: base64_str !< Base64 string representation - integer(pLongInt) :: l + integer(pI64) :: l - l = len(base64_str,pLongInt) - valid_base64 = .true. + l = len(base64_str,pI64) + validBase64 = .true. - if(mod(l,4_pLongInt)/=0_pLongInt .or. l < 4_pInt) valid_base64 = .false. - if(verify(base64_str(:l-2_pLongInt),base64_encoding, kind=pLongInt) /= 0_pLongInt) valid_base64 = .false. - if(verify(base64_str(l-1_pLongInt:),base64_encoding//'=',kind=pLongInt) /= 0_pLongInt) valid_base64 = .false. + if(mod(l,4_pI64)/=0_pI64 .or. l < 4_pInt) validBase64 = .false. + if(verify(base64_str(:l-2_pI64),base64_encoding, kind=pI64) /= 0_pI64) validBase64 = .false. + if(verify(base64_str(l-1_pI64:),base64_encoding//'=',kind=pI64) /= 0_pI64) validBase64 = .false. -end function valid_base64 +end function validBase64 !-------------------------------------------------------------------------------------------------- @@ -167,58 +167,58 @@ subroutine selfTest character(len=*), parameter :: zero_to_three = 'AAECAw==' ! https://en.wikipedia.org/wiki/Base64#Output_padding - if(base64_nBase64(20_pLongInt) /= 28_pLongInt) call IO_error(0,ext_msg='base64_nBase64/20/28') - if(base64_nBase64(19_pLongInt) /= 28_pLongInt) call IO_error(0,ext_msg='base64_nBase64/19/28') - if(base64_nBase64(18_pLongInt) /= 24_pLongInt) call IO_error(0,ext_msg='base64_nBase64/18/24') - if(base64_nBase64(17_pLongInt) /= 24_pLongInt) call IO_error(0,ext_msg='base64_nBase64/17/24') - if(base64_nBase64(16_pLongInt) /= 24_pLongInt) call IO_error(0,ext_msg='base64_nBase64/16/24') + if(base64_nChar(20_pI64) /= 28_pI64) call IO_error(0,ext_msg='base64_nChar/20/28') + if(base64_nChar(19_pI64) /= 28_pI64) call IO_error(0,ext_msg='base64_nChar/19/28') + if(base64_nChar(18_pI64) /= 24_pI64) call IO_error(0,ext_msg='base64_nChar/18/24') + if(base64_nChar(17_pI64) /= 24_pI64) call IO_error(0,ext_msg='base64_nChar/17/24') + if(base64_nChar(16_pI64) /= 24_pI64) call IO_error(0,ext_msg='base64_nChar/16/24') - if(base64_nByte(4_pLongInt) /= 3_pLongInt) call IO_error(0,ext_msg='base64_nByte/4/3') - if(base64_nByte(8_pLongInt) /= 6_pLongInt) call IO_error(0,ext_msg='base64_nByte/8/6') + if(base64_nByte(4_pI64) /= 3_pI64) call IO_error(0,ext_msg='base64_nByte/4/3') + if(base64_nByte(8_pI64) /= 6_pI64) call IO_error(0,ext_msg='base64_nByte/8/6') bytes = base64_to_bytes(zero_to_three) if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) call IO_error(0,ext_msg='base64_to_bytes//') - bytes = base64_to_bytes(zero_to_three,e=1_pLongInt) + bytes = base64_to_bytes(zero_to_three,e=1_pI64) if(any(bytes /= int([0],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes//1') - bytes = base64_to_bytes(zero_to_three,e=2_pLongInt) + bytes = base64_to_bytes(zero_to_three,e=2_pI64) if(any(bytes /= int([0,1],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes//2') - bytes = base64_to_bytes(zero_to_three,e=3_pLongInt) + bytes = base64_to_bytes(zero_to_three,e=3_pI64) if(any(bytes /= int([0,1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 3) call IO_error(0,ext_msg='base64_to_bytes//3') - bytes = base64_to_bytes(zero_to_three,e=4_pLongInt) + bytes = base64_to_bytes(zero_to_three,e=4_pI64) if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) call IO_error(0,ext_msg='base64_to_bytes//4') - bytes = base64_to_bytes(zero_to_three,s=1_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=1_pI64) if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) call IO_error(0,ext_msg='base64_to_bytes/1/') - bytes = base64_to_bytes(zero_to_three,s=2_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=2_pI64) if(any(bytes /= int([1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 3) call IO_error(0,ext_msg='base64_to_bytes/2/') - bytes = base64_to_bytes(zero_to_three,s=3_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=3_pI64) if(any(bytes /= int([2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes/3/') - bytes = base64_to_bytes(zero_to_three,s=4_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=4_pI64) if(any(bytes /= int([3],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/4/') - bytes = base64_to_bytes(zero_to_three,s=1_pLongInt,e=1_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=1_pI64,e=1_pI64) if(any(bytes /= int([0],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/1/1') - bytes = base64_to_bytes(zero_to_three,s=2_pLongInt,e=2_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=2_pI64,e=2_pI64) if(any(bytes /= int([1],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/2/2') - bytes = base64_to_bytes(zero_to_three,s=3_pLongInt,e=3_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=3_pI64,e=3_pI64) if(any(bytes /= int([2],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/3/3') - bytes = base64_to_bytes(zero_to_three,s=4_pLongInt,e=4_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=4_pI64,e=4_pI64) if(any(bytes /= int([3],C_SIGNED_CHAR)) .or. size(bytes) /= 1) call IO_error(0,ext_msg='base64_to_bytes/4/4') - bytes = base64_to_bytes(zero_to_three,s=1_pLongInt,e=2_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=1_pI64,e=2_pI64) if(any(bytes /= int([0,1],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes/1/2') - bytes = base64_to_bytes(zero_to_three,s=2_pLongInt,e=3_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=2_pI64,e=3_pI64) if(any(bytes /= int([1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes/2/3') - bytes = base64_to_bytes(zero_to_three,s=3_pLongInt,e=4_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=3_pI64,e=4_pI64) if(any(bytes /= int([2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 2) call IO_error(0,ext_msg='base64_to_bytes/3/4') - bytes = base64_to_bytes(zero_to_three,s=1_pLongInt,e=3_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=1_pI64,e=3_pI64) if(any(bytes /= int([0,1,2],C_SIGNED_CHAR)) .or. size(bytes) /= 3) call IO_error(0,ext_msg='base64_to_bytes/1/3') - bytes = base64_to_bytes(zero_to_three,s=2_pLongInt,e=4_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=2_pI64,e=4_pI64) if(any(bytes /= int([1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 3) call IO_error(0,ext_msg='base64_to_bytes/2/4') - bytes = base64_to_bytes(zero_to_three,s=1_pLongInt,e=4_pLongInt) + bytes = base64_to_bytes(zero_to_three,s=1_pI64,e=4_pI64) if(any(bytes /= int([0,1,2,3],C_SIGNED_CHAR)) .or. size(bytes) /= 4) call IO_error(0,ext_msg='base64_to_bytes/1/4') end subroutine selfTest diff --git a/src/grid/discretization_grid.f90 b/src/grid/discretization_grid.f90 index 995f9e104..7c089cc88 100644 --- a/src/grid/discretization_grid.f90 +++ b/src/grid/discretization_grid.f90 @@ -10,6 +10,8 @@ module discretization_grid use prec use system_routines + use base64 + use zlib use DAMASK_interface use IO use config @@ -64,7 +66,11 @@ subroutine discretization_grid_init(restart) write(6,'(/,a)') ' <<<+- discretization_grid init -+>>>'; flush(6) - call readGeom(grid,geomSize,origin,microstructureAt) + if(index(geometryFile,'.vtr') /= 0) then + call readVTR(grid,geomSize,origin,microstructureAt) + else + call readGeom(grid,geomSize,origin,microstructureAt) + endif !-------------------------------------------------------------------------------------------------- ! grid solver specific quantities @@ -139,10 +145,10 @@ end subroutine discretization_grid_init subroutine readGeom(grid,geomSize,origin,microstructure) integer, dimension(3), intent(out) :: & - grid ! grid (for all processes!) + grid ! grid (across all processes!) real(pReal), dimension(3), intent(out) :: & - geomSize, & ! size (for all processes!) - origin ! origin (for all processes!) + geomSize, & ! size (across all processes!) + origin ! origin (across all processes!) integer, dimension(:), intent(out), allocatable :: & microstructure @@ -150,7 +156,6 @@ subroutine readGeom(grid,geomSize,origin,microstructure) character(len=65536) :: line integer, allocatable, dimension(:) :: chunkPos integer :: & - h =- 1, & headerLength = -1, & !< length of header (in lines) fileLength, & !< length of the geom file (in characters) fileUnit, & @@ -189,7 +194,7 @@ subroutine readGeom(grid,geomSize,origin,microstructure) endif !-------------------------------------------------------------------------------------------------- -! read and interprete header +! read and interpret header origin = 0.0_pReal l = 0 do while (l < headerLength .and. startPos < len(rawData)) @@ -294,6 +299,360 @@ subroutine readGeom(grid,geomSize,origin,microstructure) end subroutine readGeom +!-------------------------------------------------------------------------------------------------- +!> @brief Parse vtk rectilinear grid (.vtr) +!> @details https://vtk.org/Wiki/VTK_XML_Formats +!-------------------------------------------------------------------------------------------------- +subroutine readVTR(grid,geomSize,origin,microstructure) + + integer, dimension(3), intent(out) :: & + grid ! grid (across all processes!) + real(pReal), dimension(3), intent(out) :: & + geomSize, & ! size (across all processes!) + origin ! origin (across all processes!) + integer, dimension(:), intent(out), allocatable :: & + microstructure + + character(len=:), allocatable :: fileContent, dataType, headerType + logical :: inFile,inGrid,gotCoordinates,gotCellData,compressed + integer :: fileUnit, myStat, coord + integer(pI64) :: & + fileLength, & !< length of the geom file (in characters) + startPos, endPos, & + s + + grid = -1 + geomSize = -1.0_pReal + +!-------------------------------------------------------------------------------------------------- +! read raw data as stream + inquire(file = trim(geometryFile), size=fileLength) + open(newunit=fileUnit, file=trim(geometryFile), access='stream',& + status='old', position='rewind', action='read',iostat=myStat) + if(myStat /= 0) call IO_error(100,ext_msg=trim(geometryFile)) + allocate(character(len=fileLength)::fileContent) + read(fileUnit) fileContent + close(fileUnit) + + inFile = .false. + inGrid = .false. + gotCoordinates = .false. + gotCelldata = .false. + +!-------------------------------------------------------------------------------------------------- +! interpret 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) then + gotCellData = .true. + startPos = endPos + 2_pI64 + do while (index(fileContent(startPos:endPos),'',kind=pI64) == 0_pI64) + endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64 + if(index(fileContent(startPos:endPos),'',kind=pI64) /= 0_pI64) then + gotCoordinates = .true. + startPos = endPos + 2_pI64 + + coord = 0 + do while (startPos',kind=pI64) /= 0_pI64) exit + startPos = endPos + 2_pI64 + enddo + endif + endif + endif + + if(gotCellData .and. gotCoordinates) exit + startPos = endPos + 2_pI64 + + end do + + if(.not. allocated(microstructure)) call IO_error(error_ID = 844, ext_msg='materialpoint not found') + if(size(microstructure) /= product(grid)) call IO_error(error_ID = 844, ext_msg='size(materialpoint)') + if(any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='size') + if(any(grid<1)) call IO_error(error_ID = 844, ext_msg='grid') + + contains + + !------------------------------------------------------------------------------------------------ + !> @brief determine size and origin from coordinates + !------------------------------------------------------------------------------------------------ + subroutine gridSizeOrigin(base64_str,headerType,compressed,dataType,direction) + + character(len=*), intent(in) :: base64_str, & ! base64 encoded string of 1D coordinates + headerType, & ! header type (UInt32 or Uint64) + dataType ! data type (Int32, Int64, Float32, Float64) + logical, intent(in) :: compressed ! indicate whether data is zlib compressed + integer, intent(in) :: direction ! direction (1=x,2=y,3=z) + + real(pReal), dimension(:), allocatable :: coords,delta + + coords = as_pReal(base64_str,headerType,compressed,dataType) + + delta = coords(2:) - coords(:size(coords)-1) + if(any(delta<0.0_pReal) .or. dNeq(maxval(delta),minval(delta))) & + call IO_error(error_ID = 844, ext_msg = 'grid spacing') + + grid(direction) = size(coords)-1 + origin(direction) = coords(1) + geomSize(direction) = coords(size(coords)) - coords(1) + + 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(bytes_to_C_INT32_T(asBytes(base64_str,headerType,compressed))) + case('Int64') + as_Int = int(bytes_to_C_INT64_T(asBytes(base64_str,headerType,compressed))) + case('Float32') + as_Int = int(bytes_to_C_FLOAT (asBytes(base64_str,headerType,compressed))) + case('Float64') + as_Int = int(bytes_to_C_DOUBLE (asBytes(base64_str,headerType,compressed))) + case default + call IO_error(844_pInt,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(bytes_to_C_INT32_T(asBytes(base64_str,headerType,compressed)),pReal) + case('Int64') + as_pReal = real(bytes_to_C_INT64_T(asBytes(base64_str,headerType,compressed)),pReal) + case('Float32') + as_pReal = real(bytes_to_C_FLOAT (asBytes(base64_str,headerType,compressed)),pReal) + case('Float64') + as_pReal = real(bytes_to_C_DOUBLE (asBytes(base64_str,headerType,compressed)),pReal) + case default + call IO_error(844_pInt,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) + endif + + 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(bytes_to_C_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(bytes_to_C_INT32_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64) + elseif(headerType == 'UInt64') then + temp = int(bytes_to_C_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(bytes_to_C_INT64_T(base64_to_bytes(base64_str(:base64_nChar(headerLen)))),pI64) + endif + + 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(0)) + e = 0_pI64 + do b = 1, nBlock + s = e + 1_pI64 + e = s + size_deflated(b) - 1_pI64 + bytes = [bytes,zlib_inflate(bytes_inflated(s:e),size_inflated(b))] + enddo + + 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(bytes_to_C_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)) + enddo + elseif(headerType == 'UInt64') then + do while(s+base64_nChar(8_pI64)<(len(base64_str,pI64))) + nByte = int(bytes_to_C_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)) + enddo + endif + + 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) + endif + endif + + end function + + + !------------------------------------------------------------------------------------------------ + !> @brief check for supported file format + !------------------------------------------------------------------------------------------------ + pure function fileFormatOk(line) + + character(len=*),intent(in) :: line + logical :: fileFormatOk + + fileFormatOk = getXMLValue(line,'type') == 'RectilinearGrid' .and. & + getXMLValue(line,'byte_order') == 'LittleEndian' .and. & + getXMLValue(line,'compressor') /= 'vtkLZ4DataCompressor' .and. & + getXMLValue(line,'compressor') /= 'vtkLZMADataCompressor' + + end function fileFormatOk + +end subroutine readVTR + + !--------------------------------------------------------------------------------------------------- !> @brief Calculate undeformed position of IPs/cell centers (pretend to be an element) !--------------------------------------------------------------------------------------------------- diff --git a/src/prec.f90 b/src/prec.f90 index 8bb3332f3..d3ec108fe 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -15,12 +15,13 @@ module prec ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds integer, parameter :: pReal = IEEE_selected_real_kind(15,307) !< number with 15 significant digits, up to 1e+-307 (typically 64 bit) + integer, parameter :: pI32 = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit) + integer, parameter :: pI64 = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) #if(INT==8) - integer, parameter :: pInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) + integer, parameter :: pInt = pI64 #else - integer, parameter :: pInt = selected_int_kind(9) !< number with at least up to +-1e9 (typically 32 bit) + integer, parameter :: pInt = pI32 #endif - integer, parameter :: pLongInt = selected_int_kind(18) !< number with at least up to +-1e18 (typically 64 bit) integer, parameter :: pStringLen = 256 !< default string length integer, parameter :: pPathLen = 4096 !< maximum length of a path name on linux @@ -237,7 +238,7 @@ end function cNeq pure function bytes_to_C_FLOAT(bytes) integer(C_SIGNED_CHAR), dimension(:), intent(in) :: bytes !< byte-wise representation of a C_FLOAT array - real(C_FLOAT), dimension(size(bytes,kind=pLongInt)/(storage_size(0._C_FLOAT,pLongInt)/8_pLongInt)) :: & + real(C_FLOAT), dimension(size(bytes,kind=pI64)/(storage_size(0._C_FLOAT,pI64)/8_pI64)) :: & bytes_to_C_FLOAT bytes_to_C_FLOAT = transfer(bytes,bytes_to_C_FLOAT,size(bytes_to_C_FLOAT)) @@ -251,7 +252,7 @@ end function bytes_to_C_FLOAT pure function bytes_to_C_DOUBLE(bytes) integer(C_SIGNED_CHAR), dimension(:), intent(in) :: bytes !< byte-wise representation of a C_DOUBLE array - real(C_DOUBLE), dimension(size(bytes,kind=pLongInt)/(storage_size(0._C_DOUBLE,pLongInt)/8_pLongInt)) :: & + real(C_DOUBLE), dimension(size(bytes,kind=pI64)/(storage_size(0._C_DOUBLE,pI64)/8_pI64)) :: & bytes_to_C_DOUBLE bytes_to_C_DOUBLE = transfer(bytes,bytes_to_C_DOUBLE,size(bytes_to_C_DOUBLE)) @@ -265,7 +266,7 @@ end function bytes_to_C_DOUBLE pure function bytes_to_C_INT32_T(bytes) integer(C_SIGNED_CHAR), dimension(:), intent(in) :: bytes !< byte-wise representation of a C_INT32_T array - integer(C_INT32_T), dimension(size(bytes,kind=pLongInt)/(storage_size(0_C_INT32_T,pLongInt)/8_pLongInt)) :: & + integer(C_INT32_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT32_T,pI64)/8_pI64)) :: & bytes_to_C_INT32_T bytes_to_C_INT32_T = transfer(bytes,bytes_to_C_INT32_T,size(bytes_to_C_INT32_T)) @@ -279,7 +280,7 @@ end function bytes_to_C_INT32_T pure function bytes_to_C_INT64_T(bytes) integer(C_SIGNED_CHAR), dimension(:), intent(in) :: bytes !< byte-wise representation of a C_INT64_T array - integer(C_INT64_T), dimension(size(bytes,kind=pLongInt)/(storage_size(0_C_INT64_T,pLongInt)/8_pLongInt)) :: & + integer(C_INT64_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT64_T,pI64)/8_pI64)) :: & bytes_to_C_INT64_T bytes_to_C_INT64_T = transfer(bytes,bytes_to_C_INT64_T,size(bytes_to_C_INT64_T)) diff --git a/src/zlib.f90 b/src/zlib.f90 index 670131334..21428255c 100644 --- a/src/zlib.f90 +++ b/src/zlib.f90 @@ -33,7 +33,7 @@ contains function zlib_inflate(deflated,size_inflated) integer(C_SIGNED_CHAR), dimension(:), intent(in) :: deflated - integer(pLongInt), intent(in) :: size_inflated + integer(pI64), intent(in) :: size_inflated integer(C_SIGNED_CHAR), dimension(size_inflated) :: zlib_inflate