out of place behavior for VTK

This commit is contained in:
Martin Diehl 2022-02-21 12:17:00 +01:00
parent c975a5b437
commit e8a4f06634
5 changed files with 32 additions and 33 deletions

View File

@ -644,8 +644,8 @@ class Grid:
Compress with zlib algorithm. Defaults to True. Compress with zlib algorithm. Defaults to True.
""" """
v = VTK.from_image_data(self.cells,self.size,self.origin) v = VTK.from_image_data(self.cells,self.size,self.origin)\
v.add(self.material.flatten(order='F'),'material') .add(self.material.flatten(order='F'),'material')
v.comments += self.comments v.comments += self.comments
v.save(fname,parallel=False,compress=compress) v.save(fname,parallel=False,compress=compress)
@ -693,8 +693,8 @@ class Grid:
Colors used to map material IDs. Colors used to map material IDs.
""" """
v = VTK.from_image_data(self.cells,self.size,self.origin) v = VTK.from_image_data(self.cells,self.size,self.origin) \
v.add(self.material.flatten(),'material') .add(self.material.flatten(),'material')
v.show('material',colormap) v.show('material',colormap)

View File

@ -1656,7 +1656,7 @@ class Result:
for inc in util.show_progress(self.visible['increments']): for inc in util.show_progress(self.visible['increments']):
u = _read(f['/'.join([inc,'geometry','u_n' if mode.lower() == 'cell' else 'u_p'])]) u = _read(f['/'.join([inc,'geometry','u_n' if mode.lower() == 'cell' else 'u_p'])])
v.add(u,'u') v = v.add(u,'u')
for ty in ['phase','homogenization']: for ty in ['phase','homogenization']:
for field in self.visible['fields']: for field in self.visible['fields']:
@ -1683,7 +1683,7 @@ class Result:
outs[out][at_cell_ho[label]] = data[in_data_ho[label]] outs[out][at_cell_ho[label]] = data[in_data_ho[label]]
for label,dataset in outs.items(): for label,dataset in outs.items():
v.add(dataset,' / '.join(['/'.join([ty,field,label]),dataset.dtype.metadata['unit']])) v = v.add(dataset,' / '.join(['/'.join([ty,field,label]),dataset.dtype.metadata['unit']]))
v.save(f'{self.fname.stem}_inc{inc[10:].zfill(N_digits)}',parallel=parallel) v.save(f'{self.fname.stem}_inc{inc[10:].zfill(N_digits)}',parallel=parallel)

View File

@ -589,7 +589,7 @@ class Table:
labels += [f'{util.srepr(self.shapes[l],"x")}:{i+1}_{l}' \ labels += [f'{util.srepr(self.shapes[l],"x")}:{i+1}_{l}' \
for i in range(np.prod(self.shapes[l]))] for i in range(np.prod(self.shapes[l]))]
f = open(fname,'w',newline='\n') if isinstance(fname, (str, Path)) else fname f = open(Path(fname).expanduser(),'w',newline='\n') if isinstance(fname, (str, Path)) else fname
f.write('\n'.join([f'# {c}' for c in self.comments] + [' '.join(labels)])+'\n') f.write('\n'.join([f'# {c}' for c in self.comments] + [' '.join(labels)])+'\n')
self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False) self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False)

View File

@ -42,11 +42,11 @@ class VTK:
if isinstance(self.vtk_data,vtk.vtkImageData): if isinstance(self.vtk_data,vtk.vtkImageData):
dup = vtk.vtkImageData() dup = vtk.vtkImageData()
elif isinstance(self.vtk_data,vtk.vtkUnstructuredGrid): elif isinstance(self.vtk_data,vtk.vtkUnstructuredGrid):
dup = vtk.vtkUnstructuredGrid dup = vtk.vtkUnstructuredGrid()
elif isinstance(self.vtk_data,vtk.vtkPolyData): elif isinstance(self.vtk_data,vtk.vtkPolyData):
dup = vtk.vtkPolyData dup = vtk.vtkPolyData()
elif isinstance(self.vtk_data,vtk.vtkRectilinearGrid): elif isinstance(self.vtk_data,vtk.vtkRectilinearGrid):
dup = vtk.vtkRectilinearGrid dup = vtk.vtkRectilinearGrid()
else: else:
raise TypeError raise TypeError
@ -222,7 +222,7 @@ class VTK:
VTK-based geometry from file. VTK-based geometry from file.
""" """
if not os.path.isfile(os.path.expanduser(fname)): # vtk has a strange error handling if not Path(fname).expanduser().is_file(): # vtk has a strange error handling
raise FileNotFoundError(f'No such file: {fname}') raise FileNotFoundError(f'No such file: {fname}')
if (ext := Path(fname).suffix) == '.vtk' or dataset_type is not None: if (ext := Path(fname).suffix) == '.vtk' or dataset_type is not None:
reader = vtk.vtkGenericDataObjectReader() reader = vtk.vtkGenericDataObjectReader()
@ -315,7 +315,6 @@ class VTK:
# Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data # Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data
# Needs support for damask.Table
def add(self, def add(self,
data: Union[np.ndarray, np.ma.MaskedArray, 'Table'], data: Union[np.ndarray, np.ma.MaskedArray, 'Table'],
label: str = None): label: str = None):
@ -332,7 +331,7 @@ class VTK:
""" """
def _add_array(self, def _add_array(vtk_data,
data: np.ndarray, data: np.ndarray,
label: str): label: str):
@ -350,27 +349,29 @@ class VTK:
d.SetName(label) d.SetName(label)
if N_data == self.N_points: if N_data == vtk_data.GetNumberOfPoints():
self.vtk_data.GetPointData().AddArray(d) vtk_data.GetPointData().AddArray(d)
elif N_data == self.N_cells: elif N_data == vtk_data.GetNumberOfCells():
self.vtk_data.GetCellData().AddArray(d) vtk_data.GetCellData().AddArray(d)
else: else:
raise ValueError(f'Data count mismatch ({N_data}{self.N_points} & {self.N_cells})') raise ValueError(f'Data count mismatch ({N_data}{self.N_points} & {self.N_cells})')
dup = self.copy()
if isinstance(data,np.ndarray): if isinstance(data,np.ndarray):
if label is not None: if label is not None:
_add_array(self, _add_array(dup.vtk_data,
np.where(data.mask,data.fill_value,data) if isinstance(data,np.ma.MaskedArray) else data, np.where(data.mask,data.fill_value,data) if isinstance(data,np.ma.MaskedArray) else data,
label) label)
else: else:
raise ValueError('No label defined for numpy.ndarray') raise ValueError('No label defined for numpy.ndarray')
elif isinstance(data,Table): elif isinstance(data,Table):
for l in data.labels: for l in data.labels:
_add_array(self,data.get(l),l) _add_array(dup.vtk_data,data.get(l),l)
else: else:
raise TypeError raise TypeError
return dup
def get(self, def get(self,
label: str) -> np.ndarray: label: str) -> np.ndarray:

View File

@ -147,8 +147,8 @@ class TestVTK:
@pytest.mark.parametrize('N_values',[5*6*7,6*7*8]) @pytest.mark.parametrize('N_values',[5*6*7,6*7*8])
def test_add_get(self,default,data_type,shape,N_values): def test_add_get(self,default,data_type,shape,N_values):
data = np.squeeze(np.random.randint(0,100,(N_values,)+shape)).astype(data_type) data = np.squeeze(np.random.randint(0,100,(N_values,)+shape)).astype(data_type)
default.add(data,'data') new = default.add(data,'data')
assert (np.squeeze(data.reshape(N_values,-1)) == default.get('data')).all() assert (np.squeeze(data.reshape(N_values,-1)) == new.get('data')).all()
@pytest.mark.parametrize('shapes',[{'scalar':(1,),'vector':(3,),'tensor':(3,3)}, @pytest.mark.parametrize('shapes',[{'scalar':(1,),'vector':(3,),'tensor':(3,3)},
@ -160,18 +160,17 @@ class TestVTK:
for k,s in shapes.items(): for k,s in shapes.items():
d[k] = dict(shape = s, d[k] = dict(shape = s,
data = np.random.random(N*np.prod(s)).reshape((N,-1))) data = np.random.random(N*np.prod(s)).reshape((N,-1)))
default.add(Table(np.column_stack([d[k]['data'] for k in shapes.keys()]),shapes)) new = default.add(Table(np.column_stack([d[k]['data'] for k in shapes.keys()]),shapes))
for k,s in shapes.items(): for k,s in shapes.items():
assert np.allclose(np.squeeze(d[k]['data']),default.get(k),rtol=1e-7) assert np.allclose(np.squeeze(d[k]['data']),new.get(k),rtol=1e-7)
def test_add_masked(self,default): def test_add_masked(self,default):
data = np.random.rand(5*6*7,3) data = np.random.rand(5*6*7,3)
masked = ma.MaskedArray(data,mask=data<.4,fill_value=42.) masked = ma.MaskedArray(data,mask=data<.4,fill_value=42.)
default.add(masked,'D') mask_auto = default.add(masked,'D')
result_masked = str(default) mask_manual = default.add(np.where(masked.mask,masked.fill_value,masked),'D')
default.add(np.where(masked.mask,masked.fill_value,masked),'D') assert str(mask_manual) == str(mask_auto)
assert result_masked == str(default)
def test_comments(self,tmp_path,default): def test_comments(self,tmp_path,default):
@ -183,8 +182,7 @@ class TestVTK:
@pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<8, reason='missing METADATA') @pytest.mark.xfail(int(vtk.vtkVersion.GetVTKVersion().split('.')[0])<8, reason='missing METADATA')
def test_compare_reference_polyData(self,update,ref_path,tmp_path): def test_compare_reference_polyData(self,update,ref_path,tmp_path):
points=np.dstack((np.linspace(0.,1.,10),np.linspace(0.,2.,10),np.linspace(-1.,1.,10))).squeeze() points=np.dstack((np.linspace(0.,1.,10),np.linspace(0.,2.,10),np.linspace(-1.,1.,10))).squeeze()
polyData = VTK.from_poly_data(points) polyData = VTK.from_poly_data(points).add(points,'coordinates')
polyData.add(points,'coordinates')
if update: if update:
polyData.save(ref_path/'polyData') polyData.save(ref_path/'polyData')
else: else:
@ -197,12 +195,12 @@ class TestVTK:
grid = [np.arange(4)**2., grid = [np.arange(4)**2.,
np.arange(5)**2., np.arange(5)**2.,
np.arange(6)**2.] # ParaView renders tetrahedral meshing unless using float coordinates! np.arange(6)**2.] # ParaView renders tetrahedral meshing unless using float coordinates!
rectilinearGrid = VTK.from_rectilinear_grid(grid)
coords = np.stack(np.meshgrid(*grid,indexing='ij'),axis=-1) coords = np.stack(np.meshgrid(*grid,indexing='ij'),axis=-1)
c = coords[:-1,:-1,:-1,:].reshape(-1,3,order='F') c = coords[:-1,:-1,:-1,:].reshape(-1,3,order='F')
n = coords[:,:,:,:].reshape(-1,3,order='F') n = coords[:,:,:,:].reshape(-1,3,order='F')
rectilinearGrid.add(np.ascontiguousarray(c),'cell') rectilinearGrid = VTK.from_rectilinear_grid(grid) \
rectilinearGrid.add(np.ascontiguousarray(n),'node') .add(np.ascontiguousarray(c),'cell') \
.add(np.ascontiguousarray(n),'node')
if update: if update:
rectilinearGrid.save(ref_path/'rectilinearGrid') rectilinearGrid.save(ref_path/'rectilinearGrid')
else: else: