Merge remote-tracking branch 'origin/development' into polishing

This commit is contained in:
Martin Diehl 2022-03-09 15:22:22 +01:00
commit c2453c56f1
8 changed files with 508 additions and 371 deletions

View File

@ -1 +1 @@
v3.0.0-alpha6-71-g3bc1a5eda v3.0.0-alpha6-88-g2162442e0

View File

@ -34,7 +34,8 @@ class Grid:
material: np.ndarray, material: np.ndarray,
size: FloatSequence, size: FloatSequence,
origin: FloatSequence = np.zeros(3), 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. New geometry definition for grid solvers.
@ -55,10 +56,10 @@ class Grid:
self.size = size # type: ignore self.size = size # type: ignore
self.origin = origin # type: ignore self.origin = origin # type: ignore
self.comments = [] if comments is None else comments # 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: def __repr__(self) -> str:
"""Basic information on grid definition.""" """Give short human-readable summary."""
mat_min = np.nanmin(self.material) mat_min = np.nanmin(self.material)
mat_max = np.nanmax(self.material) mat_max = np.nanmax(self.material)
mat_N = self.N_materials mat_N = self.N_materials
@ -68,7 +69,7 @@ class Grid:
f'origin: {util.srepr(self.origin," ")} m', f'origin: {util.srepr(self.origin," ")} m',
f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else
f' (min: {mat_min}, max: {mat_max})') 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': def __copy__(self) -> 'Grid':
@ -187,11 +188,13 @@ class Grid:
cells = np.array(v.vtk_data.GetDimensions())-1 cells = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
comments = v.comments 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'), return Grid(material = v.get('material').reshape(cells,order='F'),
size = bbox[1] - bbox[0], size = bbox[1] - bbox[0],
origin = bbox[0], origin = bbox[0],
comments=comments) comments = comments,
initial_conditions = ic)
@typing. no_type_check @typing. no_type_check
@ -660,6 +663,8 @@ class Grid:
""" """
v = VTK.from_image_data(self.cells,self.size,self.origin)\ v = VTK.from_image_data(self.cells,self.size,self.origin)\
.add(self.material.flatten(order='F'),'material') .add(self.material.flatten(order='F'),'material')
for label,data in self.ic.items():
v = v.add(data.flatten(order='F'),label)
v.comments += self.comments v.comments += self.comments
v.save(fname,parallel=False,compress=compress) v.save(fname,parallel=False,compress=compress)
@ -697,14 +702,14 @@ class Grid:
def show(self, def show(self,
colormap: Colormap = Colormap.from_predefined('cividis')) -> None: colormap: Union[Colormap, str] = 'cividis') -> None:
""" """
Show on screen. Show on screen.
Parameters Parameters
---------- ----------
colormap : damask.Colormap, optional colormap : damask.Colormap or str, optional
Colors used to map material IDs. Defaults to 'cividis'. Colormap for visualization of material IDs. Defaults to 'cividis'.
""" """
VTK.from_image_data(self.cells,self.size,self.origin) \ VTK.from_image_data(self.cells,self.size,self.origin) \

View File

@ -88,10 +88,10 @@ class VTK:
@property @property
def comments(self) -> List[str]: def comments(self) -> List[str]:
"""Return the comments.""" """Return the comments."""
fielddata = self.vtk_data.GetFieldData() field_data = self.vtk_data.GetFieldData()
for a in range(fielddata.GetNumberOfArrays()): for a in range(field_data.GetNumberOfArrays()):
if fielddata.GetArrayName(a) == 'comments': if field_data.GetArrayName(a) == 'comments':
comments = fielddata.GetAbstractArray(a) comments = field_data.GetAbstractArray(a)
return [comments.GetValue(i) for i in range(comments.GetNumberOfValues())] return [comments.GetValue(i) for i in range(comments.GetNumberOfValues())]
return [] return []
@ -498,16 +498,26 @@ class VTK:
# string array # string array
return np.array([vtk_array.GetValue(i) for i in range(vtk_array.GetNumberOfValues())]).astype(str) return np.array([vtk_array.GetValue(i) for i in range(vtk_array.GetNumberOfValues())]).astype(str)
except UnboundLocalError: except UnboundLocalError:
raise ValueError(f'array "{label}" not found') raise KeyError(f'array "{label}" not found')
def show(self, def show(self,
label: str = None, label: str = None,
colormap: Colormap = Colormap.from_predefined('cividis')): colormap: Union[Colormap, str] = 'cividis'):
""" """
Render. Render.
Parameters
----------
label : str, optional
Label of the dataset to show.
colormap : damask.Colormap or str, optional
Colormap for visualization of dataset. Defaults to 'cividis'.
Notes
-----
See http://compilatrix.com/article/vtk-1 for further ideas. See http://compilatrix.com/article/vtk-1 for further ideas.
""" """
try: try:
import wx import wx
@ -525,8 +535,10 @@ class VTK:
height = 768 height = 768
lut = vtk.vtkLookupTable() lut = vtk.vtkLookupTable()
lut.SetNumberOfTableValues(len(colormap.colors)) colormap_ = Colormap.from_predefined(colormap) if isinstance(colormap,str) else \
for i,c in enumerate(colormap.colors): 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.SetTableValue(i,c if len(c)==4 else np.append(c,1.0))
lut.Build() lut.Build()
@ -545,11 +557,11 @@ class VTK:
if label is None: if label is None:
ren.SetBackground(67/255,128/255,208/255) ren.SetBackground(67/255,128/255,208/255)
else: else:
colormap = vtk.vtkScalarBarActor() colormap_vtk = vtk.vtkScalarBarActor()
colormap.SetLookupTable(lut) colormap_vtk.SetLookupTable(lut)
colormap.SetTitle(label) colormap_vtk.SetTitle(label)
colormap.SetMaximumWidthInPixels(width//100) colormap_vtk.SetMaximumWidthInPixels(width//100)
ren.AddActor2D(colormap) ren.AddActor2D(colormap_vtk)
ren.SetBackground(0.3,0.3,0.3) ren.SetBackground(0.3,0.3,0.3)
window = vtk.vtkRenderWindow() window = vtk.vtkRenderWindow()

View File

@ -6,6 +6,7 @@ from damask import VTK
from damask import Grid from damask import Grid
from damask import Table from damask import Table
from damask import Rotation from damask import Rotation
from damask import Colormap
from damask import util from damask import util
from damask import seeds from damask import seeds
from damask import grid_filters from damask import grid_filters
@ -50,9 +51,10 @@ class TestGrid:
print('patched datetime.datetime.now') print('patched datetime.datetime.now')
def test_show(sef,default,monkeypatch): @pytest.mark.parametrize('cmap',[Colormap.from_predefined('stress'),'viridis'])
def test_show(sef,default,cmap,monkeypatch):
monkeypatch.delenv('DISPLAY',raising=False) monkeypatch.delenv('DISPLAY',raising=False)
default.show() default.show(cmap)
def test_equal(self,default): def test_equal(self,default):
assert default == default assert default == default
@ -69,7 +71,7 @@ class TestGrid:
def test_invalid_no_material(self,tmp_path): 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 = 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) v.save(tmp_path/'no_materialpoint.vti',parallel=False)
with pytest.raises(ValueError): with pytest.raises(KeyError):
Grid.load(tmp_path/'no_materialpoint.vti') Grid.load(tmp_path/'no_materialpoint.vti')
def test_invalid_material_type(self): def test_invalid_material_type(self):

View File

@ -10,6 +10,7 @@ import vtk
from damask import VTK from damask import VTK
from damask import Table from damask import Table
from damask import Colormap
@pytest.fixture @pytest.fixture
def ref_path(ref_path_base): def ref_path(ref_path_base):
@ -29,10 +30,10 @@ class TestVTK:
def _patch_execution_stamp(self, patch_execution_stamp): def _patch_execution_stamp(self, patch_execution_stamp):
print('patched damask.util.execution_stamp') print('patched damask.util.execution_stamp')
def test_show(sef,default,monkeypatch): @pytest.mark.parametrize('cmap',[Colormap.from_predefined('cividis'),'strain'])
def test_show(sef,default,cmap,monkeypatch):
monkeypatch.delenv('DISPLAY',raising=False) monkeypatch.delenv('DISPLAY',raising=False)
default.show() default.show(colormap=cmap)
def test_imageData(self,tmp_path): def test_imageData(self,tmp_path):
cells = np.random.randint(5,10,3) cells = np.random.randint(5,10,3)
@ -141,7 +142,7 @@ class TestVTK:
def test_invalid_get(self,default): def test_invalid_get(self,default):
with pytest.raises(ValueError): with pytest.raises(KeyError):
default.get('does_not_exist') default.get('does_not_exist')
def test_invalid_add_shape(self,default): def test_invalid_add_shape(self,default):

View File

@ -440,7 +440,7 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
case (155) case (155)
msg = 'material index out of bounds' msg = 'material index out of bounds'
case (180) case (180)
msg = 'missing/invalid material definition via State Variable 2' msg = 'missing/invalid material definition'
case (190) case (190)
msg = 'unknown element type:' msg = 'unknown element type:'
case (191) case (191)

446
src/VTI.f90 Normal file
View File

@ -0,0 +1,446 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Read data from image files of the visualization toolkit.
!--------------------------------------------------------------------------------------------------
module VTI
use prec
use zlib
use base64
use IO
implicit none
private
public :: &
VTI_readDataset_int, &
VTI_readDataset_real, &
VTI_readCellsSizeOrigin
contains
!--------------------------------------------------------------------------------------------------
!> @brief Read integer dataset from a VTK image data (*.vti) file.
!> @details https://vtk.org/Wiki/VTK_XML_Formats
!--------------------------------------------------------------------------------------------------
function VTI_readDataset_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 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 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
!--------------------------------------------------------------------------------------------------
subroutine VTI_readDataset_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),'<VTKFile',kind=pI64) /= 0_pI64) then
inFile = .true.
if (.not. fileFormatOk(fileContent(startPos:endPos))) call IO_error(error_ID = 844, ext_msg='file format')
headerType = merge('UInt64','UInt32',getXMLValue(fileContent(startPos:endPos),'header_type')=='UInt64')
compressed = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor'
end if
else
if (.not. inImage) then
if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
inImage = .true.
end if
else
if (index(fileContent(startPos:endPos),'<CellData',kind=pI64) /= 0_pI64) then
do while (index(fileContent(startPos:endPos),'</CellData>',kind=pI64) == 0_pI64)
if (index(fileContent(startPos:endPos),'<DataArray',kind=pI64) /= 0_pI64 .and. &
getXMLValue(fileContent(startPos:endPos),'Name') == label ) then
if (getXMLValue(fileContent(startPos:endPos),'format') /= 'binary') &
call IO_error(error_ID = 844, ext_msg='format ('//label//')')
dataType = getXMLValue(fileContent(startPos:endPos),'type')
startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
s = startPos + verify(fileContent(startPos:endPos),IO_WHITESPACE,kind=pI64) -1_pI64 ! start (no leading whitespace)
base64_str = fileContent(s:endPos)
exit outer
end if
startPos = endPos + 2_pI64
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
end do
end if
end if
end if
startPos = endPos + 2_pI64
end do outer
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 VTI_readCellsSizeOrigin(cells,geomSize,origin, &
fileContent)
integer, dimension(3), intent(out) :: &
cells ! # of cells (across all processes!)
real(pReal), dimension(3), intent(out) :: &
geomSize, & ! size (across all processes!)
origin ! origin (across all processes!)
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.
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),'<VTKFile',kind=pI64) /= 0_pI64) then
inFile = .true.
if (.not. fileFormatOk(fileContent(startPos:endPos))) call IO_error(error_ID = 844, ext_msg='file format')
headerType = merge('UInt64','UInt32',getXMLValue(fileContent(startPos:endPos),'header_type')=='UInt64')
compressed = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor'
end if
else
if (.not. inImage) then
if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
inImage = .true.
call cellsSizeOrigin(cells,geomSize,origin,fileContent(startPos:endPos))
exit outer
end if
end if
end if
startPos = endPos + 2_pI64
end do outer
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 VTI_readCellsSizeOrigin
!--------------------------------------------------------------------------------------------------
!> @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 real of kind pReal.
!--------------------------------------------------------------------------------------------------
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_real
select case(dataType)
case('Int32')
as_real = real(prec_bytesToC_INT32_T(asBytes(base64_str,headerType,compressed)),pReal)
case('Int64')
as_real = real(prec_bytesToC_INT64_T(asBytes(base64_str,headerType,compressed)),pReal)
case('Float32')
as_real = real(prec_bytesToC_FLOAT (asBytes(base64_str,headerType,compressed)),pReal)
case('Float64')
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_real
!--------------------------------------------------------------------------------------------------
!> @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 VTI

View File

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