Skip to content

Commit be1d6cb

Browse files
MarDiehlMarDiehl
committed
Merge branch '501-output-labels-of-vti-file' into 'development'
Output initial condition labels from VTI See merge request damask/DAMASK!1115
2 parents 21551b0 + a1960ad commit be1d6cb

File tree

2 files changed

+88
-34
lines changed

2 files changed

+88
-34
lines changed

‎src/grid/VTI.f90‎

Lines changed: 60 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -12,9 +12,9 @@ module VTI
1212
private
1313

1414
public :: &
15-
VTI_readDataset_int, &
1615
VTI_readDataset_real, &
17-
VTI_readCellsSizeOrigin
16+
VTI_readDataset_int, &
17+
VTI_readGeometry
1818

1919
contains
2020

@@ -109,7 +109,7 @@ subroutine VTI_readDataset_raw(base64Str,dataType,headerType,compressed, &
109109
getXMLValue(fileContent(startPos:endPos),'Name') == label ) then
110110

111111
if (getXMLValue(fileContent(startPos:endPos),'format') /= 'binary') &
112-
call IO_error(error_ID = 844, ext_msg='"'//label//'" not in binary format')
112+
call IO_error(844_pI16, label, 'not in binary format', emph = [1])
113113
dataType = getXMLValue(fileContent(startPos:endPos),'type')
114114

115115
startPos = endPos + 2_pI64
@@ -129,24 +129,26 @@ subroutine VTI_readDataset_raw(base64Str,dataType,headerType,compressed, &
129129

130130
end do outer
131131

132-
if (.not. allocated(base64Str)) call IO_error(error_ID = 844, ext_msg='dataset "'//label//'" not found')
132+
if (.not. allocated(base64Str)) call IO_error(844_pI16, 'dataset', label, 'not found', emph = [2])
133133

134134
end subroutine VTI_readDataset_raw
135135

136136

137137
!--------------------------------------------------------------------------------------------------
138-
!> @brief Read cells, size, and origin of an VTK image data (*.vti) file.
138+
!> @brief Read cells, size, and origin, and cell data labels of an VTK image data (*.vti) file.
139139
!> @details https://vtk.org/Wiki/VTK_XML_Formats
140140
!--------------------------------------------------------------------------------------------------
141-
subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, &
142-
fileContent)
141+
subroutine VTI_readGeometry(cells,geomSize,origin,labels, &
142+
fileContent)
143143

144144
integer, dimension(3), intent(out) :: &
145145
cells ! # of cells (across all processes!)
146146
real(pREAL), dimension(3), intent(out) :: &
147147
geomSize, & ! size (across all processes!)
148148
origin ! origin (across all processes!)
149-
character(len=*), intent(in) :: &
149+
character(len=pSTRLEN), allocatable, dimension(:), intent(out) :: &
150+
labels ! cell data labels
151+
character(len=*), intent(in) :: &
150152
fileContent
151153

152154
character(len=:), allocatable :: headerType
@@ -161,7 +163,8 @@ subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, &
161163
inFile = .false.
162164
inImage = .false.
163165
startPos = 1_pI64
164-
outer: do while (startPos < len(fileContent,kind=pI64))
166+
167+
do while (startPos < len(fileContent,kind=pI64))
165168
endPos = startPos + index(fileContent(startPos:),IO_EOL,kind=pI64) - 2_pI64
166169
if (endPos < startPos) endPos = len(fileContent,kind=pI64) ! end of file without new line
167170

@@ -170,26 +173,30 @@ subroutine VTI_readCellsSizeOrigin(cells,geomSize,origin, &
170173
inFile = .true.
171174
call checkFileFormat(fileContent(startPos:endPos))
172175
headerType = merge('UInt64','UInt32',getXMLValue(fileContent(startPos:endPos),'header_type')=='UInt64')
173-
compressed = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor'
176+
compressed = getXMLValue(fileContent(startPos:endPos),'compressor') == 'vtkZLibDataCompressor'
174177
end if
175178
else
176179
if (.not. inImage) then
177180
if (index(fileContent(startPos:endPos),'<ImageData',kind=pI64) /= 0_pI64) then
178181
inImage = .true.
179182
call cellsSizeOrigin(cells,geomSize,origin,fileContent(startPos:endPos))
180-
exit outer
183+
end if
184+
else
185+
if (index(fileContent(startPos:endPos),'<CellData',kind=pI64) /= 0_pI64) then
186+
call cell_labels(labels,fileContent(startPos:))
187+
exit
181188
end if
182189
end if
183190
end if
184191

185192
startPos = endPos + 2_pI64
186193

187-
end do outer
194+
end do
188195

189-
if (any(geomSize<=0)) call IO_error(error_ID = 844, ext_msg='one or more grid.size <= 0')
190-
if (any(cells<1)) call IO_error(error_ID = 844, ext_msg='one or more grid.cells < 1')
196+
if (any(geomSize <= 0)) call IO_error(844_pI16, 'one or more grid.size <= 0')
197+
if (any(cells < 1)) call IO_error(844_pI16, 'one or more grid.cells < 1')
191198

192-
end subroutine VTI_readCellsSizeOrigin
199+
end subroutine VTI_readGeometry
193200

194201

195202
!--------------------------------------------------------------------------------------------------
@@ -225,6 +232,39 @@ subroutine cellsSizeOrigin(c,s,o,header)
225232
end subroutine cellsSizeOrigin
226233

227234

235+
!--------------------------------------------------------------------------------------------------
236+
!> @brief Get labels of all cell-based datasets.
237+
!--------------------------------------------------------------------------------------------------
238+
subroutine cell_labels(labels,file_content)
239+
240+
character(len=pSTRLEN), allocatable, dimension(:), intent(out) :: labels !< labels of cell data
241+
character(len=*), intent(in) :: file_content
242+
243+
character(len=pSTRLEN) :: label
244+
integer(pI64) :: startPos, endPos
245+
246+
247+
startPos = 1_pI64
248+
endPos = startPos + index(file_content(startPos:),IO_EOL,kind=pI64) - 2_pI64
249+
250+
allocate(labels(0))
251+
252+
do while (index(file_content(startPos:endPos),'</CellData>',kind=pI64) == 0_pI64)
253+
if (index(file_content(startPos:endPos),'<DataArray',kind=pI64) /= 0_pI64) then
254+
label = getXMLValue(file_content(startPos:endPos),'Name')
255+
if (any(labels == label)) then
256+
call IO_error(844_pI16, 'repeated label', trim(label), emph = [2])
257+
else
258+
labels = [labels, label]
259+
end if
260+
end if
261+
startPos = endPos + 2_pI64
262+
endPos = startPos + index(file_content(startPos:),IO_EOL,kind=pI64) - 2_pI64
263+
end do
264+
265+
end subroutine cell_labels
266+
267+
228268
!--------------------------------------------------------------------------------------------------
229269
!> @brief Interpret Base64 string in vtk XML file as integer of default kind.
230270
!--------------------------------------------------------------------------------------------------
@@ -248,7 +288,7 @@ function as_Int(base64Str,headerType,compressed,dataType)
248288
case('Float64')
249289
as_Int = int(prec_bytesToC_DOUBLE (asBytes(base64Str,headerType,compressed)))
250290
case default
251-
call IO_error(844,ext_msg='unknown data type: '//trim(dataType))
291+
call IO_error(844_pI16, 'unknown data type',trim(dataType),emph=[2])
252292
end select
253293

254294
end function as_Int
@@ -277,7 +317,7 @@ function as_real(base64Str,headerType,compressed,dataType)
277317
case('Float64')
278318
as_real = real(prec_bytesToC_DOUBLE (asBytes(base64Str,headerType,compressed)),pREAL)
279319
case default
280-
call IO_error(844,ext_msg='unknown data type: '//trim(dataType))
320+
call IO_error(844_pI16,'unknown data type',trim(dataType),emph=[2])
281321
end select
282322

283323
end function as_real
@@ -436,15 +476,15 @@ subroutine checkFileFormat(line)
436476

437477
val = getXMLValue(line,'type')
438478
if (val /= 'ImageData') &
439-
call IO_error(844, ext_msg='type ("'//val//'") is not "ImageData"')
479+
call IO_error(844_pI16, 'type',val,'is not "ImageData"',emph=[2])
440480

441481
val = getXMLValue(line,'byte_order')
442482
if (val /= 'LittleEndian') &
443-
call IO_error(844, ext_msg='byte_order ("'//val//'") is not "LittleEndian"')
483+
call IO_error(844_pI16, 'byte_order',val,'is not "LittleEndian"',emph=[2])
444484

445485
val = getXMLValue(line,'compressor')
446486
if (val /= '' .and. val /= 'vtkZLibDataCompressor') &
447-
call IO_error(844, ext_msg='compressor ("'//val//'") is not "vtkZLibDataCompressor"')
487+
call IO_error(844_pI16, 'compressor',val,'is not "vtkZLibDataCompressor"',emph=[2])
448488

449489
end subroutine checkFileFormat
450490

‎src/grid/discretization_grid.f90‎

Lines changed: 28 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -64,14 +64,17 @@ subroutine discretization_grid_init()
6464
materialAt, materialAt_global
6565

6666
integer :: &
67-
j
67+
j, &
68+
n_labels !< number cell datasets in VTI file
6869
integer(MPI_INTEGER_KIND) :: err_MPI
6970
integer(C_INTPTR_T) :: &
7071
devNull, cells3_, cells3Offset_
7172
integer, dimension(worldsize) :: &
7273
displs, sendcounts
7374
character(len=:), allocatable :: &
7475
fileContent, fname
76+
character(len=pSTRLEN), dimension(:), allocatable :: &
77+
labels !< cell data labels in VTI file
7578
integer(HID_T) :: handle
7679

7780

@@ -80,12 +83,13 @@ subroutine discretization_grid_init()
8083

8184
if (worldrank == 0) then
8285
fileContent = IO_read(CLI_geomFile)
83-
call VTI_readCellsSizeOrigin(cells,geomSize,origin,fileContent)
86+
call VTI_readGeometry(cells,geomSize,origin,labels,fileContent)
87+
n_labels = size(labels)
8488
materialAt_global = VTI_readDataset_int(fileContent,'material') + 1
8589
if (any(materialAt_global < 1)) &
86-
call IO_error(180,ext_msg='material ID < 1')
90+
call IO_error(180_pI16,'material ID < 1')
8791
if (size(materialAt_global) /= product(cells)) &
88-
call IO_error(180,ext_msg='mismatch in # of material IDs and cells')
92+
call IO_error(180_pI16,'mismatch in # of material IDs and cells')
8993
fname = CLI_geomFile
9094
if (scan(fname,'/') /= 0) fname = fname(scan(fname,'/',.true.)+1:)
9195
call result_openJobFile(parallel=.false.)
@@ -95,27 +99,37 @@ subroutine discretization_grid_init()
9599
allocate(materialAt_global(0)) ! needed for IntelMPI
96100
end if
97101

98-
99-
call MPI_Bcast(cells,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
102+
call MPI_Bcast(cells,3_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
103+
call parallelization_chkerr(err_MPI)
104+
if (cells(1) < 2) call IO_error(844_pI16,'cells(1) must be larger than 1')
105+
call MPI_Bcast(geomSize,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
100106
call parallelization_chkerr(err_MPI)
101-
if (cells(1) < 2) call IO_error(844, ext_msg='cells(1) must be larger than 1')
102-
call MPI_Bcast(geomSize,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
107+
call MPI_Bcast(origin,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
103108
call parallelization_chkerr(err_MPI)
104-
call MPI_Bcast(origin,3_MPI_INTEGER_KIND,MPI_DOUBLE,0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
109+
call MPI_Bcast(n_labels,1_MPI_INTEGER_KIND,MPI_INTEGER,0_MPI_INTEGER_KIND,MPI_COMM_WORLD,err_MPI)
110+
call parallelization_chkerr(err_MPI)
111+
112+
if (worldrank /= 0) allocate(character(len=pSTRLEN) :: labels(n_labels))
113+
call MPI_Bcast(labels,int(pSTRLEN*n_labels,MPI_INTEGER_KIND),MPI_CHARACTER, &
114+
0_MPI_INTEGER_KIND,MPI_COMM_WORLD, err_MPI)
105115
call parallelization_chkerr(err_MPI)
106116

107-
print'(/,1x,a,i0,a,i0,a,i0)', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3)
108-
print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' m³'
109-
print '(1x,a,es9.2,a,es9.2,a,es9.2,a)', 'origin: ', origin(1), ' ', origin(2), ' ', origin(3), ' m'
117+
print'(/,1x,3(a,i0))', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3)
118+
print'( 1x,3(a,es9.2),a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' m³'
119+
print'( 1x,3(a,es9.2),a)', 'origin: ', origin(1), ' ', origin(2), ' ', origin(3), ' m'
120+
print'(/,1x,a)', 'cell data:'
121+
do j = 1, n_labels
122+
print '(2x,a,a)', '- ', trim(labels(j))
123+
end do
110124

111-
if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)')
125+
if (worldsize>cells(3)) call IO_error(894_pI16,'number of processes exceeds cells(3)')
112126

113127
call fftw_mpi_init()
114128
devNull = fftw_mpi_local_size_3d(int(cells(3),C_INTPTR_T),int(cells(2),C_INTPTR_T),int(cells(1)/2+1,C_INTPTR_T), &
115129
PETSC_COMM_WORLD, &
116130
cells3_, & ! domain cells size along z
117131
cells3Offset_) ! domain cells offset along z
118-
if (cells3_==0_C_INTPTR_T) call IO_error(894, ext_msg='Cannot distribute MPI processes')
132+
if (cells3_==0_C_INTPTR_T) call IO_error(894_pI16,'Cannot distribute MPI processes')
119133

120134
cells3 = int(cells3_)
121135
cells3Offset = int(cells3Offset_)

0 commit comments

Comments
 (0)