From b2abaef0b35fad97b2c1857c6b7bf07d7393c9a7 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 9 Mar 2020 18:28:25 -0400 Subject: [PATCH 01/19] added placing of data within geometry --- python/damask/result.py | 114 +++++++++++++++++++++++++++++++++------- 1 file changed, 96 insertions(+), 18 deletions(-) diff --git a/python/damask/result.py b/python/damask/result.py index f58b57f23..c40206f89 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -11,6 +11,7 @@ import numpy as np from . import util from . import version +from . import table from . import mechanics from . import Rotation from . import Orientation @@ -151,12 +152,11 @@ class Result(): selected = [] for i,time in enumerate(self.times): if start <= time <= end: - selected.append(self.increments[i]) + selected.append(self.times[i]) return selected - - def iter_selection(self,what): + def iterate(self,what): """ Iterate over selection items by setting each one selected. @@ -226,6 +226,83 @@ class Result(): self._manage_selection('del',what,datasets) + # def createGeometry4VTK(): + # """ reads geometry and fummels... to return VTK.XXXobject (pointcloud fallback)""" + + # myVTK = new + # myVTK.geom = result.getGeom() + # myVTK.table/data = result.spatiocondense(['f','p']) # of type damask.Table + # myVTK.writeme + + # def datamerger(regular expression to filter groups into one copy) + + + def place(self,datasets,component=0,tagged=False,split=True): + """ + Distribute datasets onto geometry and return Table or (split) dictionary of Tables. + + Must not mix nodal end cell data. + + Only data within + - inc?????/constituent/*_*/* + - inc?????/materialpoint/*_*/* + - inc?????/geometry/* + are considered. + + Parameters + ---------- + datasets : iterable or str + component : int + homogenization component to consider for constituent data + tagged : Boolean + tag Table.column name with '#component' + defaults to False + split : Boolean + split Table by increment and return dictionary of Tables + defaults to True + + """ + sets = datasets if hasattr(datasets,'__iter__') and not isinstance(datasets,str) \ + else [datasets] + tag = f'#{component}' if tagged else '' + tbl = {} if split else None + inGeom = {} + inData = {} + with h5py.File(self.fname,'r') as f: + for dataset in sets: + for group in self.groups_with_datasets(dataset): + path = os.path.join(group,dataset) + inc,prop,name,cat,item = (path.split('/') + ['']*5)[:5] + key = '/'.join([prop,name+tag]) + if key not in inGeom: + if prop == 'geometry': + inGeom[key] = inData[key] = np.arange(self.Nmaterialpoints) + elif prop == 'constituent': + inGeom[key] = np.where(f['mapping/cellResults/constituent'][:,component]['Name'] == str.encode(name))[0] + inData[key] = f['mapping/cellResults/constituent'][inGeom[key],component]['Position'] + else: + inGeom[key] = np.where(f['mapping/cellResults/materialpoint']['Name'] == str.encode(name))[0] + inData[key] = f['mapping/cellResults/materialpoint'][inGeom[key].tolist()]['Position'] + shape = np.shape(f[path]) + data = np.full((self.Nmaterialpoints,) + (shape[1:] if len(shape)>1 else (1,)), + np.nan, + dtype=np.dtype(f[path])) + data[inGeom[key]] = (f[path] if len(shape)>1 else np.expand_dims(f[path],1))[inData[key]] + path = (os.path.join(*([prop,name]+([cat] if cat else [])+([item] if item else []))) if split else path)+tag + if split: + try: + tbl[inc].add(path,data) + except KeyError: + tbl[inc] = table.Table(data.reshape(self.Nmaterialpoints,-1),{path:data.shape[1:]}) + else: + try: + tbl.add(path,data) + except AttributeError: + tbl = table.Table(data.reshape(self.Nmaterialpoints,-1),{path:data.shape[1:]}) + + return tbl + + def groups_with_datasets(self,datasets): """ Return groups that contain all requested datasets. @@ -262,10 +339,10 @@ class Result(): groups = [] with h5py.File(self.fname,'r') as f: - for i in self.iter_selection('increments'): + for i in self.iterate('increments'): for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): - for oo in self.iter_selection(o): - for pp in self.iter_selection(p): + for oo in self.iterate(o): + for pp in self.iterate(p): group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue if sets is True: groups.append(group) @@ -279,12 +356,12 @@ class Result(): """Return information on all active datasets in the file.""" message = '' with h5py.File(self.fname,'r') as f: - for i in self.iter_selection('increments'): + for i in self.iterate('increments'): message+='\n{} ({}s)\n'.format(i,self.times[self.increments.index(i)]) for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): - for oo in self.iter_selection(o): + for oo in self.iterate(o): message+=' {}\n'.format(oo) - for pp in self.iter_selection(p): + for pp in self.iterate(p): message+=' {}\n'.format(pp) group = '/'.join([i,o[:-1],oo,pp]) # o[:-1]: plural/singular issue for d in f[group].keys(): @@ -301,7 +378,7 @@ class Result(): """Return the location of all active datasets with given label.""" path = [] with h5py.File(self.fname,'r') as f: - for i in self.iter_selection('increments'): + for i in self.iterate('increments'): k = '/'.join([i,'geometry',label]) try: f[k] @@ -309,8 +386,8 @@ class Result(): except KeyError as e: pass for o,p in zip(['constituents','materialpoints'],['con_physics','mat_physics']): - for oo in self.iter_selection(o): - for pp in self.iter_selection(p): + for oo in self.iterate(o): + for pp in self.iterate(p): k = '/'.join([i,o[:-1],oo,pp,label]) try: f[k] @@ -1013,15 +1090,15 @@ class Result(): N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1 - for i,inc in enumerate(self.iter_selection('increments')): + for i,inc in enumerate(self.iterate('increments')): vtk_data = [] materialpoints_backup = self.selection['materialpoints'].copy() self.pick('materialpoints',False) for label in (labels if isinstance(labels,list) else [labels]): - for p in self.iter_selection('con_physics'): + for p in self.iterate('con_physics'): if p != 'generic': - for c in self.iter_selection('constituents'): + for c in self.iterate('constituents'): x = self.get_dataset_location(label) if len(x) == 0: continue @@ -1048,9 +1125,9 @@ class Result(): constituents_backup = self.selection['constituents'].copy() self.pick('constituents',False) for label in (labels if isinstance(labels,list) else [labels]): - for p in self.iter_selection('mat_physics'): + for p in self.iterate('mat_physics'): if p != 'generic': - for m in self.iter_selection('materialpoints'): + for m in self.iterate('materialpoints'): x = self.get_dataset_location(label) if len(x) == 0: continue @@ -1095,7 +1172,8 @@ class Result(): ################################################################################################### # BEGIN DEPRECATED - iter_visible = iter_selection + iter_visible = iterate + iter_selection = iterate def _time_to_inc(self,start,end): From bffce1ab9cf0882274970e3b040a4160754a420a Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 9 Mar 2020 18:50:27 -0400 Subject: [PATCH 02/19] use "nan" to represent np.nan in ASCIItable output --- python/damask/table.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/table.py b/python/damask/table.py index a4badb3dc..b704bdc24 100644 --- a/python/damask/table.py +++ b/python/damask/table.py @@ -350,4 +350,4 @@ class Table(): f = fname for line in header + [' '.join(labels)]: f.write(line+'\n') - self.data.to_csv(f,sep=' ',index=False,header=False) + self.data.to_csv(f,sep=' ',na_rep='nan',index=False,header=False) From 8a4bc3dda4e14227d4ec856d517b3816796db725 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 10 Mar 2020 23:52:02 +0100 Subject: [PATCH 03/19] separating vtk from results we should discuss the naming! --- python/damask/ktv.py | 57 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 57 insertions(+) create mode 100644 python/damask/ktv.py diff --git a/python/damask/ktv.py b/python/damask/ktv.py new file mode 100644 index 000000000..e14d80310 --- /dev/null +++ b/python/damask/ktv.py @@ -0,0 +1,57 @@ +import numpy as np +import vtk +#from vtk.util import numpy_support + +class VTK: # capitals needed/preferred? + """ + Manage vtk files. + + tbd + """ + + def __init__(self,geom): + """tbd.""" + self.geom = geom + + @staticmethod + def from_rectilinearGrid(grid,size,origin=np.zeros(3)): + coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] + for dim in [0,1,2]: + for c in np.linspace(0,size[dim],1+grid[dim]): + coordArray[dim].InsertNextValue(c) + + geom = vtk.vtkRectilinearGrid() + geom.SetDimensions(*(grid+1)) + geom.SetXCoordinates(coordArray[0]) + geom.SetYCoordinates(coordArray[1]) + geom.SetZCoordinates(coordArray[2]) + + return VTK(geom) + + + @staticmethod + def from_unstructuredGrid(nodes,connectivity,elem): + geom = vtk.vtkUnstructuredGrid() + geom.SetPoints(nodes) + geom.Allocate(connectivity.shape[0]) + + if elem == 'TRIANGLE': + vtk_type = vtk.VTK_TRIANGLE + n_nodes = 3 + elif elem == 'QUAD': + vtk_type = vtk.VTK_QUAD + n_nodes = 4 + elif elem == 'TETRA': + vtk_type = vtk.VTK_TETRA + n_nodes = 4 + elif elem == 'HEXAHEDRON': + vtk_type = vtk.VTK_HEXAHEDRON + n_nodes = 8 + + for i in connectivity: + geom.InsertNextCell(vtk_type,n_nodes,i-1) + + return VTK(geom) + + def write(self,fname): + print('tbd',fname) From 9878ddc5503e037bcfaacf242bf5dc0fdf058635 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 06:28:13 +0100 Subject: [PATCH 04/19] easier way to show data. needs information on geometry (structured/unstructured) and probably we should list not more than 5 incs --- python/damask/result.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/python/damask/result.py b/python/damask/result.py index 2c0972ffe..47803ae58 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -84,6 +84,11 @@ class Result: self.fname = fname + def __repr__(self): + """Show selected data.""" + return util.srepr(self.list_data()) + + def _manage_selection(self,action,what,datasets): """ Manages the visibility of the groups. From a024ec378a4b1c35af662e298fce14127dc4e78e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 06:49:14 +0100 Subject: [PATCH 05/19] new class should be accesible as damask.VTK --- python/damask/__init__.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 4aa28853e..d0ebd45a7 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -9,20 +9,23 @@ with open(os.path.join(os.path.dirname(__file__),'VERSION')) as f: # classes from .environment import Environment # noqa from .table import Table # noqa -from .asciitable import ASCIItable # noqa - -from .config import Material # noqa +from .ktv import VTK # noqa from .colormaps import Colormap, Color # noqa from .rotation import Rotation # noqa from .lattice import Symmetry, Lattice# noqa from .orientation import Orientation # noqa from .result import Result # noqa -from .result import Result as DADF5 # noqa - from .geom import Geom # noqa from .solver import Solver # noqa -from .test import Test # noqa + +# compatibility hack +from .result import Result as DADF5 # noqa + +# deprecated +from .asciitable import ASCIItable # noqa from .util import extendableOption # noqa +from .config import Material # noqa +from .test import Test # noqa # functions in modules from . import mechanics # noqa From b3e8a4405e0952fb338de06b0e01d8bdd035a28f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 06:50:13 +0100 Subject: [PATCH 06/19] better use centralized functionality --- python/damask/result.py | 49 +++++++++-------------------------------- 1 file changed, 10 insertions(+), 39 deletions(-) diff --git a/python/damask/result.py b/python/damask/result.py index 47803ae58..3c1cb15a1 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -9,6 +9,7 @@ from vtk.util import numpy_support import h5py import numpy as np +from . import VTK from . import util from . import version from . import table @@ -245,15 +246,15 @@ class Result: def place(self,datasets,component=0,tagged=False,split=True): """ Distribute datasets onto geometry and return Table or (split) dictionary of Tables. - + Must not mix nodal end cell data. - + Only data within - inc?????/constituent/*_*/* - inc?????/materialpoint/*_*/* - inc?????/geometry/* are considered. - + Parameters ---------- datasets : iterable or str @@ -1040,44 +1041,14 @@ class Result: if mode.lower()=='cell': if self.structured: - coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] - for dim in [0,1,2]: - for c in np.linspace(0,self.size[dim],1+self.grid[dim]): - coordArray[dim].InsertNextValue(c) - - vtk_geom = vtk.vtkRectilinearGrid() - vtk_geom.SetDimensions(*(self.grid+1)) - vtk_geom.SetXCoordinates(coordArray[0]) - vtk_geom.SetYCoordinates(coordArray[1]) - vtk_geom.SetZCoordinates(coordArray[2]) + v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) + vtk_geom = v.geom else: - nodes = vtk.vtkPoints() with h5py.File(self.fname,'r') as f: - nodes.SetData(numpy_support.numpy_to_vtk(f['/geometry/x_n'][()],deep=True)) - - vtk_geom = vtk.vtkUnstructuredGrid() - vtk_geom.SetPoints(nodes) - vtk_geom.Allocate(f['/geometry/T_c'].shape[0]) - - if self.version_major == 0 and self.version_minor <= 5: - vtk_type = vtk.VTK_HEXAHEDRON - n_nodes = 8 - else: - if f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TRIANGLE': - vtk_type = vtk.VTK_TRIANGLE - n_nodes = 3 - elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'QUAD': - vtk_type = vtk.VTK_QUAD - n_nodes = 4 - elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'TETRA': # not tested - vtk_type = vtk.VTK_TETRA - n_nodes = 4 - elif f['/geometry/T_c'].attrs['VTK_TYPE'] == b'HEXAHEDRON': - vtk_type = vtk.VTK_HEXAHEDRON - n_nodes = 8 - - for i in f['/geometry/T_c']: - vtk_geom.InsertNextCell(vtk_type,n_nodes,i-1) + v = VTK.from_unstructuredGrid(f['/geometry/x_n'][()], + f['/geometry/T_c'], + f['/geometry/T_c'].attrs['VTK_TYPE'].decode()) + vtk_geom = v.geom elif mode.lower()=='point': Points = vtk.vtkPoints() From 32734e7dce6b04110c10e9f7aa90e6c560ff7c45 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 07:32:03 +0100 Subject: [PATCH 07/19] using central functionality --- python/damask/geom.py | 56 +++++++++++++++---------------------------- python/damask/ktv.py | 32 ++++++++++++++++++++++--- 2 files changed, 48 insertions(+), 40 deletions(-) diff --git a/python/damask/geom.py b/python/damask/geom.py index 9a1677191..dc6196157 100644 --- a/python/damask/geom.py +++ b/python/damask/geom.py @@ -6,6 +6,7 @@ from scipy import ndimage import vtk from vtk.util import numpy_support +from . import VTK from . import util from . import version @@ -36,7 +37,7 @@ class Geom(): self.set_origin(origin) self.set_homogenization(homogenization) self.set_comments(comments) - + def __repr__(self): """Basic information on geometry definition.""" @@ -70,7 +71,7 @@ class Geom(): origin_old = self.get_origin() unique_old = len(np.unique(self.microstructure)) max_old = np.nanmax(self.microstructure) - + if size is not None and rescale: raise ValueError('Either set size explicitly or rescale automatically') @@ -108,7 +109,7 @@ class Geom(): if max_old != np.nanmax(self.microstructure): message[-1] = util.delete(message[-1]) message.append(util.emph('max microstructure: {}'.format(np.nanmax(self.microstructure)))) - + return util.return_message(message) def set_comments(self,comments): @@ -123,7 +124,7 @@ class Geom(): """ self.comments = [] self.add_comments(comments) - + def add_comments(self,comments): """ Appends comments to existing comments. @@ -241,7 +242,7 @@ class Geom(): header.append('origin x {} y {} z {}'.format(*self.get_origin())) header.append('homogenization {}'.format(self.get_homogenization())) return header - + @staticmethod def from_file(fname): """ @@ -266,7 +267,7 @@ class Geom(): if not keyword.startswith('head') or header_length < 3: raise TypeError('Header length information missing or invalid') - comments = [] + comments = [] for i,line in enumerate(content[:header_length]): items = line.lower().strip().split() key = items[0] if items else '' @@ -295,14 +296,14 @@ class Geom(): else: items = list(map(float,items)) microstructure[i:i+len(items)] = items i += len(items) - + if i != grid.prod(): raise TypeError('Invalid file: expected {} entries, found {}'.format(grid.prod(),i)) - + microstructure = microstructure.reshape(grid,order='F') if not np.any(np.mod(microstructure.flatten(),1) != 0.0): # no float present microstructure = microstructure.astype('int') - + return Geom(microstructure.reshape(grid),size,origin,homogenization,comments) @@ -320,7 +321,7 @@ class Geom(): """ header = self.get_header() grid = self.get_grid() - + if pack is None: plain = grid.prod()/np.unique(self.microstructure).size < 250 else: @@ -371,7 +372,7 @@ class Geom(): elif compressType == 'of': f.write('{} of {}\n'.format(reps,former)) - + def to_vtk(self,fname=None): """ Generates vtk file. @@ -382,28 +383,9 @@ class Geom(): vtk file to write. If no file is given, a string is returned. """ - grid = self.get_grid() + np.ones(3,dtype=int) - size = self.get_size() - origin = self.get_origin() + v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) + rGrid = v.geom - coords = [ - np.linspace(0,size[0],grid[0]) + origin[0], - np.linspace(0,size[1],grid[1]) + origin[1], - np.linspace(0,size[2],grid[2]) + origin[2] - ] - - rGrid = vtk.vtkRectilinearGrid() - coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] - - rGrid.SetDimensions(*grid) - for d,coord in enumerate(coords): - for c in coord: - coordArray[d].InsertNextValue(c) - - rGrid.SetXCoordinates(coordArray[0]) - rGrid.SetYCoordinates(coordArray[1]) - rGrid.SetZCoordinates(coordArray[2]) - ms = numpy_support.numpy_to_vtk(num_array=self.microstructure.flatten(order='F'), array_type=vtk.VTK_INT if self.microstructure.dtype == int else vtk.VTK_FLOAT) ms.SetName('microstructure') @@ -418,7 +400,7 @@ class Geom(): writer = vtk.vtkXMLRectilinearGridWriter() writer.SetCompressorTypeToZLib() writer.SetDataModeToBinary() - + ext = os.path.splitext(fname)[1] if ext == '': name = fname + '.' + writer.GetDefaultFileExtension() @@ -427,13 +409,13 @@ class Geom(): else: raise ValueError("unknown extension {}".format(ext)) writer.SetFileName(name) - + writer.SetInputData(rGrid) writer.Write() - if fname is None: return writer.GetOutputString() + if not fname: return writer.GetOutputString() + - def show(self): """Show raw content (as in file).""" f=StringIO() @@ -469,7 +451,7 @@ class Geom(): ms = np.concatenate([ms,ms[:,limits[0]:limits[1]:-1,:]],1) if 'x' in directions: ms = np.concatenate([ms,ms[limits[0]:limits[1]:-1,:,:]],0) - + return self.update(ms,rescale=True) #self.add_comments('tbd') diff --git a/python/damask/ktv.py b/python/damask/ktv.py index e14d80310..08e8f21a9 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -1,5 +1,7 @@ +import os import numpy as np import vtk +import vtkmodules #from vtk.util import numpy_support class VTK: # capitals needed/preferred? @@ -15,6 +17,7 @@ class VTK: # capitals needed/preferred? @staticmethod def from_rectilinearGrid(grid,size,origin=np.zeros(3)): + """Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data.""" coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] for dim in [0,1,2]: for c in np.linspace(0,size[dim],1+grid[dim]): @@ -25,7 +28,7 @@ class VTK: # capitals needed/preferred? geom.SetXCoordinates(coordArray[0]) geom.SetYCoordinates(coordArray[1]) geom.SetZCoordinates(coordArray[2]) - + return VTK(geom) @@ -53,5 +56,28 @@ class VTK: # capitals needed/preferred? return VTK(geom) - def write(self,fname): - print('tbd',fname) + + def write(self,fname): #ToDo: Discuss how to handle consistently filename extensions + if (isinstance(self.geom,vtkmodules.vtkCommonDataModel.vtkRectilinearGrid)): + writer = vtk.vtkXMLRectilinearGridWriter() + elif(isinstance(self.geom,vtkmodules.vtkCommonDataModel.vtkUnstructuredGrid)): + writer = vtk.vtkUnstructuredGrid() + elif(isinstance(self.geom,vtkmodules.vtkCommonDataModel.vtkPolyData)): + writer = vtk.vtkXMLPolyDataWriter() + + writer.SetFileName('{}.{}'.format(os.path.splitext(fname)[0], + writer.GetDefaultFileExtension())) + writer.SetCompressorTypeToZLib() + writer.SetDataModeToBinary() + writer.SetInputData(self.geom) + + writer.Write() + + def __repr__(self): + """ASCII representation of the VTK data.""" + writer = vtk.vtkDataSetWriter() + #writer.SetHeader('damask.Geom '+version) + writer.WriteToOutputStringOn() + writer.SetInputData(self.geom) + writer.Write() + return writer.GetOutputString() From c92a6ad45991f0e0404911ad27811b0b90a977fd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 07:50:11 +0100 Subject: [PATCH 08/19] simpler and system independent --- python/damask/ktv.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/python/damask/ktv.py b/python/damask/ktv.py index 08e8f21a9..2eb1008a4 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -1,7 +1,6 @@ import os import numpy as np import vtk -import vtkmodules #from vtk.util import numpy_support class VTK: # capitals needed/preferred? @@ -58,11 +57,11 @@ class VTK: # capitals needed/preferred? def write(self,fname): #ToDo: Discuss how to handle consistently filename extensions - if (isinstance(self.geom,vtkmodules.vtkCommonDataModel.vtkRectilinearGrid)): + if (isinstance(self.geom,vtk.vtkRectilinearGrid)): writer = vtk.vtkXMLRectilinearGridWriter() - elif(isinstance(self.geom,vtkmodules.vtkCommonDataModel.vtkUnstructuredGrid)): + elif(isinstance(self.geom,vtk.vtkUnstructuredGrid)): writer = vtk.vtkUnstructuredGrid() - elif(isinstance(self.geom,vtkmodules.vtkCommonDataModel.vtkPolyData)): + elif(isinstance(self.geom,vtk.vtkPolyData)): writer = vtk.vtkXMLPolyDataWriter() writer.SetFileName('{}.{}'.format(os.path.splitext(fname)[0], From f324e67f7be64a43fc1ce029087e3af2351a6a39 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 18:07:31 +0100 Subject: [PATCH 09/19] fix: nodes need to be converted to vtk type some stub definitions --- python/damask/ktv.py | 21 ++++++++++++++++++--- 1 file changed, 18 insertions(+), 3 deletions(-) diff --git a/python/damask/ktv.py b/python/damask/ktv.py index 2eb1008a4..7afbffcc6 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -1,7 +1,12 @@ import os + +import pandas as pd import numpy as np import vtk -#from vtk.util import numpy_support +from vtk.util import numpy_support + +from . import table +from . import version class VTK: # capitals needed/preferred? """ @@ -34,7 +39,7 @@ class VTK: # capitals needed/preferred? @staticmethod def from_unstructuredGrid(nodes,connectivity,elem): geom = vtk.vtkUnstructuredGrid() - geom.SetPoints(nodes) + geom.SetPoints(numpy_support.numpy_to_vtk(nodes)) #,deep=True) geom.Allocate(connectivity.shape[0]) if elem == 'TRIANGLE': @@ -72,10 +77,20 @@ class VTK: # capitals needed/preferred? writer.Write() + + def add(data,label=None): + if isinstance(data,np.ndarray): + pass + elif isinstance(data,pd.DataFrame): + pass + elif isinstance(data,table): + pass + + def __repr__(self): """ASCII representation of the VTK data.""" writer = vtk.vtkDataSetWriter() - #writer.SetHeader('damask.Geom '+version) + writer.SetHeader('DAMASK.VTK v{}'.format(version)) writer.WriteToOutputStringOn() writer.SetInputData(self.geom) writer.Write() From 744e3bb50bc309e53c59c62ce9ebc813c98d4270 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 20:56:58 +0100 Subject: [PATCH 10/19] vectorized cell assignment + bugfix for writing out --- python/damask/ktv.py | 51 ++++++++++++++++++++++------------------- python/damask/result.py | 4 ++-- 2 files changed, 30 insertions(+), 25 deletions(-) diff --git a/python/damask/ktv.py b/python/damask/ktv.py index 7afbffcc6..5e1919300 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -3,7 +3,7 @@ import os import pandas as pd import numpy as np import vtk -from vtk.util import numpy_support +from vtk.util import numpy_support as nps from . import table from . import version @@ -24,11 +24,11 @@ class VTK: # capitals needed/preferred? """Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data.""" coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] for dim in [0,1,2]: - for c in np.linspace(0,size[dim],1+grid[dim]): + for c in np.linspace(origin[dim],origin[dim]+size[dim],grid[dim]): coordArray[dim].InsertNextValue(c) geom = vtk.vtkRectilinearGrid() - geom.SetDimensions(*(grid+1)) + geom.SetDimensions(*grid) geom.SetXCoordinates(coordArray[0]) geom.SetYCoordinates(coordArray[1]) geom.SetZCoordinates(coordArray[2]) @@ -37,35 +37,40 @@ class VTK: # capitals needed/preferred? @staticmethod - def from_unstructuredGrid(nodes,connectivity,elem): + def from_unstructuredGrid(nodes,connectivity,cell_type): + """ + Create an unstructured grid (mesh). + + connectivity: 1 based at the moment + cell_type: TRIANGLE, 'QUAD', 'TETRA','HEXAHEDRON' + + """ + vtk_nodes = vtk.vtkPoints() + vtk_nodes.SetData(nps.numpy_to_vtk(nodes)) + + cells = vtk.vtkCellArray() + cells.SetNumberOfCells(connectivity.shape[0]) + T = np.concatenate((np.ones((connectivity.shape[0],1),dtype=np.int64)*connectivity.shape[1], + connectivity),axis=1).ravel() + cells.SetCells(connectivity.shape[0],nps.numpy_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE)) + geom = vtk.vtkUnstructuredGrid() - geom.SetPoints(numpy_support.numpy_to_vtk(nodes)) #,deep=True) - geom.Allocate(connectivity.shape[0]) - - if elem == 'TRIANGLE': - vtk_type = vtk.VTK_TRIANGLE - n_nodes = 3 - elif elem == 'QUAD': - vtk_type = vtk.VTK_QUAD - n_nodes = 4 - elif elem == 'TETRA': - vtk_type = vtk.VTK_TETRA - n_nodes = 4 - elif elem == 'HEXAHEDRON': - vtk_type = vtk.VTK_HEXAHEDRON - n_nodes = 8 - - for i in connectivity: - geom.InsertNextCell(vtk_type,n_nodes,i-1) + geom.SetPoints(vtk_nodes) + geom.SetCells(eval('vtk.VTK_{}'.format(cell_type.upper())),cells) return VTK(geom) + @staticmethod + def from_points(nodes,connectivity,cell_type): + pass + + def write(self,fname): #ToDo: Discuss how to handle consistently filename extensions if (isinstance(self.geom,vtk.vtkRectilinearGrid)): writer = vtk.vtkXMLRectilinearGridWriter() elif(isinstance(self.geom,vtk.vtkUnstructuredGrid)): - writer = vtk.vtkUnstructuredGrid() + writer = vtk.vtkXMLUnstructuredGridWriter() elif(isinstance(self.geom,vtk.vtkPolyData)): writer = vtk.vtkXMLPolyDataWriter() diff --git a/python/damask/result.py b/python/damask/result.py index 3c1cb15a1..3ac442f5c 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -1041,12 +1041,12 @@ class Result: if mode.lower()=='cell': if self.structured: - v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) + v = VTK.from_rectilinearGrid(self.grid+1,self.size,self.origin) vtk_geom = v.geom else: with h5py.File(self.fname,'r') as f: v = VTK.from_unstructuredGrid(f['/geometry/x_n'][()], - f['/geometry/T_c'], + f['/geometry/T_c'][()]-1, f['/geometry/T_c'].attrs['VTK_TYPE'].decode()) vtk_geom = v.geom From 575da581a9403a0a50e5b5d8d382374ecce02f71 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 22:35:58 +0100 Subject: [PATCH 11/19] support for point cloud --- python/damask/ktv.py | 22 +++++++++++++++++----- python/damask/result.py | 25 +++++++------------------ 2 files changed, 24 insertions(+), 23 deletions(-) diff --git a/python/damask/ktv.py b/python/damask/ktv.py index 5e1919300..67425fd73 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -41,13 +41,12 @@ class VTK: # capitals needed/preferred? """ Create an unstructured grid (mesh). - connectivity: 1 based at the moment + connectivity: 0 based at the moment, shape Ncell x N nodes cell_type: TRIANGLE, 'QUAD', 'TETRA','HEXAHEDRON' """ vtk_nodes = vtk.vtkPoints() vtk_nodes.SetData(nps.numpy_to_vtk(nodes)) - cells = vtk.vtkCellArray() cells.SetNumberOfCells(connectivity.shape[0]) T = np.concatenate((np.ones((connectivity.shape[0],1),dtype=np.int64)*connectivity.shape[1], @@ -62,8 +61,21 @@ class VTK: # capitals needed/preferred? @staticmethod - def from_points(nodes,connectivity,cell_type): - pass + def from_polyData(points): + vtk_points= vtk.vtkPoints() + vtk_points.SetData(nps.numpy_to_vtk(points)) + + vertices = vtk.vtkCellArray() + vertices.SetNumberOfCells(points.shape[0]) + T = np.concatenate((np.ones((points.shape[0],1),dtype=np.int64), + np.arange(points.shape[0],dtype=np.int64).reshape(-1,1)),axis=1).ravel() + vertices.SetCells(points.shape[0],nps.numpy_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE)) + + geom = vtk.vtkPolyData() + geom.SetPoints(vtk_points) + geom.SetVerts(vertices) + + return VTK(geom) def write(self,fname): #ToDo: Discuss how to handle consistently filename extensions @@ -95,7 +107,7 @@ class VTK: # capitals needed/preferred? def __repr__(self): """ASCII representation of the VTK data.""" writer = vtk.vtkDataSetWriter() - writer.SetHeader('DAMASK.VTK v{}'.format(version)) + writer.SetHeader('# DAMASK.VTK v{}'.format(version)) writer.WriteToOutputStringOn() writer.SetInputData(self.geom) writer.Write() diff --git a/python/damask/result.py b/python/damask/result.py index 3ac442f5c..9f4c21cf6 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -975,11 +975,11 @@ class Result: datasets_in = {} lock.acquire() with h5py.File(self.fname,'r') as f: - for arg,label in datasets.items(): - loc = f[group+'/'+label] - datasets_in[arg]={'data' :loc[()], - 'label':label, - 'meta': {k:v.decode() for k,v in loc.attrs.items()}} + for arg,label in datasets.items(): + loc = f[group+'/'+label] + datasets_in[arg]={'data' :loc[()], + 'label':label, + 'meta': {k:v.decode() for k,v in loc.attrs.items()}} lock.release() r = func(**datasets_in,**args) return [group,r] @@ -1042,26 +1042,15 @@ class Result: if self.structured: v = VTK.from_rectilinearGrid(self.grid+1,self.size,self.origin) - vtk_geom = v.geom else: with h5py.File(self.fname,'r') as f: v = VTK.from_unstructuredGrid(f['/geometry/x_n'][()], f['/geometry/T_c'][()]-1, f['/geometry/T_c'].attrs['VTK_TYPE'].decode()) - vtk_geom = v.geom - elif mode.lower()=='point': - Points = vtk.vtkPoints() - Vertices = vtk.vtkCellArray() - for c in self.cell_coordinates(): - pointID = Points.InsertNextPoint(c) - Vertices.InsertNextCell(1) - Vertices.InsertCellPoint(pointID) + v = VTK.from_polyData(self.cell_coordinates()) - vtk_geom = vtk.vtkPolyData() - vtk_geom.SetPoints(Points) - vtk_geom.SetVerts(Vertices) - vtk_geom.Modified() + vtk_geom = v.geom N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1 From bb2019810a9cfba14be5e5206cb5b945be39220c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 11 Mar 2020 23:54:36 +0100 Subject: [PATCH 12/19] centralizing functionality --- python/damask/geom.py | 37 +++------------------- python/damask/ktv.py | 48 ++++++++++++++++++++++++++-- python/damask/result.py | 69 ++++++++--------------------------------- 3 files changed, 63 insertions(+), 91 deletions(-) diff --git a/python/damask/geom.py b/python/damask/geom.py index dc6196157..a2b5249f2 100644 --- a/python/damask/geom.py +++ b/python/damask/geom.py @@ -1,14 +1,11 @@ -import os +import sys from io import StringIO import numpy as np from scipy import ndimage -import vtk -from vtk.util import numpy_support from . import VTK from . import util -from . import version class Geom(): @@ -384,36 +381,12 @@ class Geom(): """ v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) - rGrid = v.geom + v.add(self.microstructure.flatten(order='F'),'microstructure') - ms = numpy_support.numpy_to_vtk(num_array=self.microstructure.flatten(order='F'), - array_type=vtk.VTK_INT if self.microstructure.dtype == int else vtk.VTK_FLOAT) - ms.SetName('microstructure') - rGrid.GetCellData().AddArray(ms) - - - if fname is None: - writer = vtk.vtkDataSetWriter() - writer.SetHeader('damask.Geom '+version) - writer.WriteToOutputStringOn() + if fname: + v.write(fname) else: - writer = vtk.vtkXMLRectilinearGridWriter() - writer.SetCompressorTypeToZLib() - writer.SetDataModeToBinary() - - ext = os.path.splitext(fname)[1] - if ext == '': - name = fname + '.' + writer.GetDefaultFileExtension() - elif ext[1:] == writer.GetDefaultFileExtension(): - name = fname - else: - raise ValueError("unknown extension {}".format(ext)) - writer.SetFileName(name) - - writer.SetInputData(rGrid) - writer.Write() - - if not fname: return writer.GetOutputString() + sys.stdout.write(v.__repr__()) def show(self): diff --git a/python/damask/ktv.py b/python/damask/ktv.py index 67425fd73..de8392339 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -77,8 +77,40 @@ class VTK: # capitals needed/preferred? return VTK(geom) + @staticmethod + def from_file(fname,ftype=None): + ext = os.path.splitext(fname)[1] + if ext == '.vtk': + reader = vtk.vtkGenericDataObjectReader() + reader.SetFileName(fname) + reader.Update() + if ftype.lower() == 'rectilineargrid': + geom = reader.GetRectilinearGridOutput() + elif ftype.lower() == 'unstructuredgrid': + geom = reader.GetUnstructuredGridOutput() + elif ftype.lower() == 'polydata': + geom = reader.GetPolyDataOutput() + else: + raise Exception + else: + if ext == '.vtr': + reader = vtk.vtkXMLRectilinearGridReader() + elif ext == '.vtu': + reader = vtk.vtkXMLUnstructuredGridReader() + elif ext == '.vtp': + reader = vtk.vtkXMLPolyDataReader() + else: + raise Exception - def write(self,fname): #ToDo: Discuss how to handle consistently filename extensions + reader.SetFileName(fname) + reader.Update() + geom = reader.GetOutput() + + return VTK(geom) + + + def write(self,fname): + """ToDo: Check if given fileextension makes sense.""" if (isinstance(self.geom,vtk.vtkRectilinearGrid)): writer = vtk.vtkXMLRectilinearGridWriter() elif(isinstance(self.geom,vtk.vtkUnstructuredGrid)): @@ -95,9 +127,19 @@ class VTK: # capitals needed/preferred? writer.Write() - def add(data,label=None): + def add(self,data,label=None): + + Npoints = self.geom.GetNumberOfPoints() + Ncells = self.geom.GetNumberOfCells() + if isinstance(data,np.ndarray): - pass + shape = [data.shape[0],np.product(data.shape[1:],dtype=np.int)] + d = nps.numpy_to_vtk(num_array=data.reshape(shape),deep=True) + d.SetName(label) + if shape[0] == Ncells: + self.geom.GetCellData().AddArray(d) + elif shape[0] == Npoints: + self.geom.GetPointData().AddArray(d) elif isinstance(data,pd.DataFrame): pass elif isinstance(data,table): diff --git a/python/damask/result.py b/python/damask/result.py index 9f4c21cf6..20ee803fd 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -4,8 +4,6 @@ import glob import os from functools import partial -import vtk -from vtk.util import numpy_support import h5py import numpy as np @@ -77,10 +75,9 @@ class Result: self.mat_physics = list(set(self.mat_physics)) # make unique self.selection= {'increments': self.increments, - 'constituents': self.constituents, - 'materialpoints': self.materialpoints, - 'con_physics': self.con_physics, - 'mat_physics': self.mat_physics} + 'constituents': self.constituents,'materialpoints': self.materialpoints, + 'con_physics': self.con_physics, 'mat_physics': self.mat_physics + } self.fname = fname @@ -231,15 +228,6 @@ class Result: """ self._manage_selection('del',what,datasets) - - # def createGeometry4VTK(): - # """ reads geometry and fummels... to return VTK.XXXobject (pointcloud fallback)""" - - # myVTK = new - # myVTK.geom = result.getGeom() - # myVTK.table/data = result.spatiocondense(['f','p']) # of type damask.Table - # myVTK.writeme - # def datamerger(regular expression to filter groups into one copy) @@ -1050,12 +1038,9 @@ class Result: elif mode.lower()=='point': v = VTK.from_polyData(self.cell_coordinates()) - vtk_geom = v.geom - - N_digits = int(np.floor(np.log10(int(self.increments[-1][3:]))))+1 + N_digits = int(np.floor(np.log10(min(int(self.increments[-1][3:]),1))))+1 for i,inc in enumerate(self.iterate('increments')): - vtk_data = [] materialpoints_backup = self.selection['materialpoints'].copy() self.pick('materialpoints',False) @@ -1067,23 +1052,15 @@ class Result: if len(x) == 0: continue array = self.read_dataset(x,0) - shape = [array.shape[0],np.product(array.shape[1:])] - vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True)) - vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1! - vtk_geom.GetCellData().AddArray(vtk_data[-1]) - + v.add(array,'1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1! else: x = self.get_dataset_location(label) if len(x) == 0: continue array = self.read_dataset(x,0) - shape = [array.shape[0],np.product(array.shape[1:])] - vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True)) ph_name = re.compile(r'(?<=(constituent\/))(.*?)(?=(generic))') # identify phase name dset_name = '1_' + re.sub(ph_name,r'',x[0].split('/',1)[1]) # removing phase name - vtk_data[-1].SetName(dset_name) - vtk_geom.GetCellData().AddArray(vtk_data[-1]) - + v.add(array,dset_name) self.pick('materialpoints',materialpoints_backup) constituents_backup = self.selection['constituents'].copy() @@ -1096,43 +1073,23 @@ class Result: if len(x) == 0: continue array = self.read_dataset(x,0) - shape = [array.shape[0],np.product(array.shape[1:])] - vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True)) - vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) #ToDo: why 1_? - vtk_geom.GetCellData().AddArray(vtk_data[-1]) + v.add(array,'1_'+x[0].split('/',1)[1]) #ToDo: why 1_? else: x = self.get_dataset_location(label) if len(x) == 0: continue array = self.read_dataset(x,0) - shape = [array.shape[0],np.product(array.shape[1:])] - vtk_data.append(numpy_support.numpy_to_vtk(num_array=array.reshape(shape),deep=True)) - vtk_data[-1].SetName('1_'+x[0].split('/',1)[1]) - vtk_geom.GetCellData().AddArray(vtk_data[-1]) + v.add(array,'1_'+x[0].split('/',1)[1]) self.pick('constituents',constituents_backup) if mode.lower()=='cell': - writer = vtk.vtkXMLRectilinearGridWriter() if self.structured else \ - vtk.vtkXMLUnstructuredGridWriter() - x = self.get_dataset_location('u_n') - vtk_data.append(numpy_support.numpy_to_vtk(num_array=self.read_dataset(x,0),deep=True)) - vtk_data[-1].SetName('u') - vtk_geom.GetPointData().AddArray(vtk_data[-1]) - elif mode.lower()=='point': - writer = vtk.vtkXMLPolyDataWriter() + u = self.read_dataset(self.get_dataset_location('u_n')) + v.add(u,'u') + file_out = '{}_inc{}'.format(os.path.splitext(os.path.basename(self.fname))[0], + inc[3:].zfill(N_digits)) - file_out = '{}_inc{}.{}'.format(os.path.splitext(os.path.basename(self.fname))[0], - inc[3:].zfill(N_digits), - writer.GetDefaultFileExtension()) - - writer.SetCompressorTypeToZLib() - writer.SetDataModeToBinary() - writer.SetFileName(file_out) - writer.SetInputData(vtk_geom) - - writer.Write() - + v.write(file_out) ################################################################################################### # BEGIN DEPRECATED From a6a73cdc0f07c9eb162aa880ff2e9bb3e6cf9fb5 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 12 Mar 2020 06:51:52 +0100 Subject: [PATCH 13/19] polishing grid is interpreted again in the DAMASK meaning, i.e it specifies the number of cells, not the number of nodes --- python/damask/ktv.py | 106 ++++++++++++++++++++++++++++++++-------- python/damask/result.py | 6 +-- 2 files changed, 88 insertions(+), 24 deletions(-) diff --git a/python/damask/ktv.py b/python/damask/ktv.py index de8392339..5cef8af21 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -8,27 +8,50 @@ from vtk.util import numpy_support as nps from . import table from . import version -class VTK: # capitals needed/preferred? +class VTK: """ - Manage vtk files. + Spatial visualization (and potentially manipulation). - tbd + High-level interface to VTK. """ def __init__(self,geom): - """tbd.""" + """ + Set geometry and topology. + + Parameters + ---------- + geom : subclass of vtk.vtkDataSet + Description of geometry and topology. Valid types are vtk.vtkRectilinearGrid, + vtk.vtkUnstructuredGrid, or vtk.vtkPolyData. + + """ self.geom = geom @staticmethod def from_rectilinearGrid(grid,size,origin=np.zeros(3)): - """Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data.""" + """ + Create VTK of type vtk.vtkRectilinearGrid. + + This is the common type for results from the grid solver. + + Parameters + ---------- + grid : numpy.ndarray of shape (3) of np.dtype = int + Number of cells. + size : numpy.ndarray of shape (3) + Physical length. + origin : numpy.ndarray of shape (3), optional + Spatial origin. + + """ coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] for dim in [0,1,2]: - for c in np.linspace(origin[dim],origin[dim]+size[dim],grid[dim]): + for c in np.linspace(origin[dim],origin[dim]+size[dim],grid[dim]+1): coordArray[dim].InsertNextValue(c) geom = vtk.vtkRectilinearGrid() - geom.SetDimensions(*grid) + geom.SetDimensions(*(grid+1)) geom.SetXCoordinates(coordArray[0]) geom.SetYCoordinates(coordArray[1]) geom.SetZCoordinates(coordArray[2]) @@ -39,10 +62,18 @@ class VTK: # capitals needed/preferred? @staticmethod def from_unstructuredGrid(nodes,connectivity,cell_type): """ - Create an unstructured grid (mesh). + Create VTK of type vtk.vtkUnstructuredGrid. - connectivity: 0 based at the moment, shape Ncell x N nodes - cell_type: TRIANGLE, 'QUAD', 'TETRA','HEXAHEDRON' + This is the common type for results from FEM solvers. + + Parameters + ---------- + nodes : numpy.ndarray of shape (:,3) + Spatial position of the nodes. + connectivity : numpy.ndarray of np.dtype = int + Cell connectivity (0-based), first dimension determines #Cells, second dimension determines #Nodes/Cell. + cell_type : str + Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, and HEXAHEDRON. """ vtk_nodes = vtk.vtkPoints() @@ -62,6 +93,17 @@ class VTK: # capitals needed/preferred? @staticmethod def from_polyData(points): + """ + Create VTK of type vtk.polyData. + + This is the common type for point-wise data. + + Parameters + ---------- + points : numpy.ndarray of shape (:,3) + Spatial position of the points. + + """ vtk_points= vtk.vtkPoints() vtk_points.SetData(nps.numpy_to_vtk(points)) @@ -79,16 +121,28 @@ class VTK: # capitals needed/preferred? @staticmethod def from_file(fname,ftype=None): + """ + Create VTK from file. + + Parameters + ---------- + fname : str + Filename for reading. Valid extensions are .vtk, .vtr, .vtu, and .vtp. + ftype : str, optional + Name of the vtk.vtkDataSet subclass when opening an .vtk file. Valid types are vtkRectilinearGrid, + vtkUnstructuredGrid, and vtkPolyData. + + """ ext = os.path.splitext(fname)[1] if ext == '.vtk': reader = vtk.vtkGenericDataObjectReader() reader.SetFileName(fname) reader.Update() - if ftype.lower() == 'rectilineargrid': + if 'rectilineargrid' in ftype.lower(): geom = reader.GetRectilinearGridOutput() - elif ftype.lower() == 'unstructuredgrid': + elif 'unstructuredgrid' in ftype.lower(): geom = reader.GetUnstructuredGridOutput() - elif ftype.lower() == 'polydata': + elif 'polydata' in ftype.lower(): geom = reader.GetPolyDataOutput() else: raise Exception @@ -109,8 +163,17 @@ class VTK: # capitals needed/preferred? return VTK(geom) + # ToDo: If extension is given, check for consistency. def write(self,fname): - """ToDo: Check if given fileextension makes sense.""" + """ + Write to file. + + Parameters + ---------- + fname : str + Filename for writing. + + """ if (isinstance(self.geom,vtk.vtkRectilinearGrid)): writer = vtk.vtkXMLRectilinearGridWriter() elif(isinstance(self.geom,vtk.vtkUnstructuredGrid)): @@ -127,18 +190,19 @@ class VTK: # capitals needed/preferred? writer.Write() + # Check https://blog.kitware.com/ghost-and-blanking-visibility-changes/ for missing data + # Needs support for pd.DataFrame and/or table def add(self,data,label=None): - - Npoints = self.geom.GetNumberOfPoints() - Ncells = self.geom.GetNumberOfCells() + """Add data to either cells or points.""" + N_points = self.geom.GetNumberOfPoints() + N_cells = self.geom.GetNumberOfCells() if isinstance(data,np.ndarray): - shape = [data.shape[0],np.product(data.shape[1:],dtype=np.int)] - d = nps.numpy_to_vtk(num_array=data.reshape(shape),deep=True) + d = nps.numpy_to_vtk(num_array=data.reshape(data.shape[0],-1),deep=True) d.SetName(label) - if shape[0] == Ncells: + if data.shape[0] == N_cells: self.geom.GetCellData().AddArray(d) - elif shape[0] == Npoints: + elif data.shape[0] == N_points: self.geom.GetPointData().AddArray(d) elif isinstance(data,pd.DataFrame): pass diff --git a/python/damask/result.py b/python/damask/result.py index 20ee803fd..37bd90b42 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -26,12 +26,12 @@ class Result: def __init__(self,fname): """ - Opens an existing DADF5 file. + Open an existing DADF5 file. Parameters ---------- fname : str - name of the DADF5 file to be openend. + name of the DADF5 file to be openend. """ with h5py.File(fname,'r') as f: @@ -1029,7 +1029,7 @@ class Result: if mode.lower()=='cell': if self.structured: - v = VTK.from_rectilinearGrid(self.grid+1,self.size,self.origin) + v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) else: with h5py.File(self.fname,'r') as f: v = VTK.from_unstructuredGrid(f['/geometry/x_n'][()], From dafc48dcc74485adf41431a20bf10c61e2794283 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 12 Mar 2020 08:15:12 +0100 Subject: [PATCH 14/19] polishing --- python/damask/ktv.py | 37 +++++++++++++++++++------------------ 1 file changed, 19 insertions(+), 18 deletions(-) diff --git a/python/damask/ktv.py b/python/damask/ktv.py index 5cef8af21..1a05a6611 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -3,7 +3,7 @@ import os import pandas as pd import numpy as np import vtk -from vtk.util import numpy_support as nps +from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk from . import table from . import version @@ -47,8 +47,8 @@ class VTK: """ coordArray = [vtk.vtkDoubleArray(),vtk.vtkDoubleArray(),vtk.vtkDoubleArray()] for dim in [0,1,2]: - for c in np.linspace(origin[dim],origin[dim]+size[dim],grid[dim]+1): - coordArray[dim].InsertNextValue(c) + coords = np.linspace(origin[dim],origin[dim]+size[dim],grid[dim]+1) + coordArray[dim].SetArray(np_to_vtk(coords),grid[dim]+1,1) geom = vtk.vtkRectilinearGrid() geom.SetDimensions(*(grid+1)) @@ -77,16 +77,16 @@ class VTK: """ vtk_nodes = vtk.vtkPoints() - vtk_nodes.SetData(nps.numpy_to_vtk(nodes)) + vtk_nodes.SetData(np_to_vtk(nodes)) cells = vtk.vtkCellArray() cells.SetNumberOfCells(connectivity.shape[0]) T = np.concatenate((np.ones((connectivity.shape[0],1),dtype=np.int64)*connectivity.shape[1], - connectivity),axis=1).ravel() - cells.SetCells(connectivity.shape[0],nps.numpy_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE)) + connectivity),axis=1).ravel() + cells.SetCells(connectivity.shape[0],np_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE)) geom = vtk.vtkUnstructuredGrid() geom.SetPoints(vtk_nodes) - geom.SetCells(eval('vtk.VTK_{}'.format(cell_type.upper())),cells) + geom.SetCells(eval('vtk.VTK_{}'.format(cell_type.split('_',1)[-1].upper())),cells) return VTK(geom) @@ -105,13 +105,13 @@ class VTK: """ vtk_points= vtk.vtkPoints() - vtk_points.SetData(nps.numpy_to_vtk(points)) + vtk_points.SetData(np_to_vtk(points)) vertices = vtk.vtkCellArray() vertices.SetNumberOfCells(points.shape[0]) T = np.concatenate((np.ones((points.shape[0],1),dtype=np.int64), - np.arange(points.shape[0],dtype=np.int64).reshape(-1,1)),axis=1).ravel() - vertices.SetCells(points.shape[0],nps.numpy_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE)) + np.arange(points.shape[0],dtype=np.int64).reshape(-1,1)),axis=1).ravel() + vertices.SetCells(points.shape[0],np_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE)) geom = vtk.vtkPolyData() geom.SetPoints(vtk_points) @@ -119,8 +119,9 @@ class VTK: return VTK(geom) + @staticmethod - def from_file(fname,ftype=None): + def from_file(fname,dataset_type=None): """ Create VTK from file. @@ -128,7 +129,7 @@ class VTK: ---------- fname : str Filename for reading. Valid extensions are .vtk, .vtr, .vtu, and .vtp. - ftype : str, optional + dataset_type : str, optional Name of the vtk.vtkDataSet subclass when opening an .vtk file. Valid types are vtkRectilinearGrid, vtkUnstructuredGrid, and vtkPolyData. @@ -138,14 +139,14 @@ class VTK: reader = vtk.vtkGenericDataObjectReader() reader.SetFileName(fname) reader.Update() - if 'rectilineargrid' in ftype.lower(): + if 'rectilineargrid' in dataset_type.lower(): geom = reader.GetRectilinearGridOutput() - elif 'unstructuredgrid' in ftype.lower(): + elif 'unstructuredgrid' in dataset_type.lower(): geom = reader.GetUnstructuredGridOutput() - elif 'polydata' in ftype.lower(): + elif 'polydata' in dataset_type.lower(): geom = reader.GetPolyDataOutput() else: - raise Exception + raise TypeError('Unknown dataset type for vtk file {}'.format(dataset_type)) else: if ext == '.vtr': reader = vtk.vtkXMLRectilinearGridReader() @@ -154,7 +155,7 @@ class VTK: elif ext == '.vtp': reader = vtk.vtkXMLPolyDataReader() else: - raise Exception + raise TypeError('Unknown file extension {}'.format(ext)) reader.SetFileName(fname) reader.Update() @@ -198,7 +199,7 @@ class VTK: N_cells = self.geom.GetNumberOfCells() if isinstance(data,np.ndarray): - d = nps.numpy_to_vtk(num_array=data.reshape(data.shape[0],-1),deep=True) + d = np_to_vtk(num_array=data.reshape(data.shape[0],-1),deep=True) d.SetName(label) if data.shape[0] == N_cells: self.geom.GetCellData().AddArray(d) From 4c915eddbcc2a2fd89b3bbe8227e50c6c9c08b6a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 12 Mar 2020 08:30:33 +0100 Subject: [PATCH 15/19] inform the user --- python/damask/result.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/result.py b/python/damask/result.py index 37bd90b42..b2938a68e 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -1040,7 +1040,7 @@ class Result: N_digits = int(np.floor(np.log10(min(int(self.increments[-1][3:]),1))))+1 - for i,inc in enumerate(self.iterate('increments')): + for i,inc in enumerate(util.show_progress(self.iterate('increments'),len(self.selection['increments']))): materialpoints_backup = self.selection['materialpoints'].copy() self.pick('materialpoints',False) From 827f354435cb81b88f541fc8ea049683bb7259db Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 12 Mar 2020 12:28:28 +0100 Subject: [PATCH 16/19] only store points, no vertices visualize via 'Points Gaussian', not 'Points' in paraview. + adding displacements for points --- python/damask/ktv.py | 7 ------- python/damask/result.py | 39 +++++++++++++++++++-------------------- 2 files changed, 19 insertions(+), 27 deletions(-) diff --git a/python/damask/ktv.py b/python/damask/ktv.py index 1a05a6611..c4c7fa567 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -107,15 +107,8 @@ class VTK: vtk_points= vtk.vtkPoints() vtk_points.SetData(np_to_vtk(points)) - vertices = vtk.vtkCellArray() - vertices.SetNumberOfCells(points.shape[0]) - T = np.concatenate((np.ones((points.shape[0],1),dtype=np.int64), - np.arange(points.shape[0],dtype=np.int64).reshape(-1,1)),axis=1).ravel() - vertices.SetCells(points.shape[0],np_to_vtk(T, deep=True, array_type=vtk.VTK_ID_TYPE)) - geom = vtk.vtkPolyData() geom.SetPoints(vtk_points) - geom.SetVerts(vertices) return VTK(geom) diff --git a/python/damask/result.py b/python/damask/result.py index b2938a68e..7d6a06441 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -446,7 +446,7 @@ class Result: def cell_coordinates(self): """Return initial coordinates of the cell centers.""" if self.structured: - return grid_filters.cell_coord0(self.grid,self.size,self.origin) + return grid_filters.cell_coord0(self.grid,self.size,self.origin).reshape(-1,3) else: with h5py.File(self.fname,'r') as f: return f['geometry/x_c'][()] @@ -1045,14 +1045,14 @@ class Result: materialpoints_backup = self.selection['materialpoints'].copy() self.pick('materialpoints',False) for label in (labels if isinstance(labels,list) else [labels]): - for p in self.iterate('con_physics'): - if p != 'generic': - for c in self.iterate('constituents'): - x = self.get_dataset_location(label) - if len(x) == 0: - continue - array = self.read_dataset(x,0) - v.add(array,'1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1! + for p in self.iterate('con_physics'): + if p != 'generic': + for c in self.iterate('constituents'): + x = self.get_dataset_location(label) + if len(x) == 0: + continue + array = self.read_dataset(x,0) + v.add(array,'1_'+x[0].split('/',1)[1]) #ToDo: hard coded 1! else: x = self.get_dataset_location(label) if len(x) == 0: @@ -1066,14 +1066,14 @@ class Result: constituents_backup = self.selection['constituents'].copy() self.pick('constituents',False) for label in (labels if isinstance(labels,list) else [labels]): - for p in self.iterate('mat_physics'): - if p != 'generic': - for m in self.iterate('materialpoints'): - x = self.get_dataset_location(label) - if len(x) == 0: - continue - array = self.read_dataset(x,0) - v.add(array,'1_'+x[0].split('/',1)[1]) #ToDo: why 1_? + for p in self.iterate('mat_physics'): + if p != 'generic': + for m in self.iterate('materialpoints'): + x = self.get_dataset_location(label) + if len(x) == 0: + continue + array = self.read_dataset(x,0) + v.add(array,'1_'+x[0].split('/',1)[1]) #ToDo: why 1_? else: x = self.get_dataset_location(label) if len(x) == 0: @@ -1082,9 +1082,8 @@ class Result: v.add(array,'1_'+x[0].split('/',1)[1]) self.pick('constituents',constituents_backup) - if mode.lower()=='cell': - u = self.read_dataset(self.get_dataset_location('u_n')) - v.add(u,'u') + u = self.read_dataset(self.get_dataset_location('u_n' if mode.lower() == 'cell' else 'u_p')) + v.add(u,'u') file_out = '{}_inc{}'.format(os.path.splitext(os.path.basename(self.fname))[0], inc[3:].zfill(N_digits)) From 6fbace8220e02d36d971fdcca9a5de8725395c54 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 12 Mar 2020 13:45:47 +0100 Subject: [PATCH 17/19] stub: show rendered geometry --- python/damask/ktv.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/python/damask/ktv.py b/python/damask/ktv.py index c4c7fa567..12252e093 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -212,3 +212,31 @@ class VTK: writer.SetInputData(self.geom) writer.Write() return writer.GetOutputString() + + + def show(self): + """ + Render. + + See http://compilatrix.com/article/vtk-1 for further ideas. + """ + mapper = vtk.vtkDataSetMapper() + mapper.SetInputData(self.geom) + actor = vtk.vtkActor() + actor.SetMapper(mapper) + + ren = vtk.vtkRenderer() + + renWin = vtk.vtkRenderWindow() + renWin.AddRenderer(ren) + + ren.AddActor(actor) + ren.SetBackground(0.2,0.2,0.2) + renWin.SetSize(1024, 1024) + + iren = vtk.vtkRenderWindowInteractor() + iren.SetRenderWindow(renWin) + + iren.Initialize() + renWin.Render() + iren.Start() From 81e98055dd3f6df8abf50bc4251a0854eb54a7d9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 12 Mar 2020 19:52:33 +0100 Subject: [PATCH 18/19] polishing --- python/damask/environment.py | 22 +++++++++++++++------- python/damask/ktv.py | 9 ++++++--- python/damask/result.py | 13 +++++++------ 3 files changed, 28 insertions(+), 16 deletions(-) diff --git a/python/damask/environment.py b/python/damask/environment.py index 4c9ab762c..b5feb4fcd 100644 --- a/python/damask/environment.py +++ b/python/damask/environment.py @@ -1,24 +1,32 @@ +import tkinter import os class Environment(): - __slots__ = [ \ - 'options', - ] def __init__(self): """Read and provide values of DAMASK configuration.""" - self.options = {} - self.__get_options() + self.options = self._get_options() + tk = tkinter.Tk() + self.screen_width = tk.winfo_screenwidth() + self.screen_height = tk.winfo_screenheight() + def relPath(self,relative = '.'): + """Return absolute path from path relative to DAMASK root.""" return os.path.join(self.rootDir(),relative) + def rootDir(self): + """Return DAMASK root path.""" return os.path.normpath(os.path.join(os.path.realpath(__file__),'../../../')) - def __get_options(self): + + def _get_options(self): + options = {} for item in ['DAMASK_NUM_THREADS', 'MSC_ROOT', 'MARC_VERSION', ]: - self.options[item] = os.environ[item] if item in os.environ else None + options[item] = os.environ[item] if item in os.environ else None + + return options diff --git a/python/damask/ktv.py b/python/damask/ktv.py index 12252e093..19b87efdf 100644 --- a/python/damask/ktv.py +++ b/python/damask/ktv.py @@ -5,9 +5,11 @@ import numpy as np import vtk from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk -from . import table +from . import Table +from . import Environment from . import version + class VTK: """ Spatial visualization (and potentially manipulation). @@ -200,7 +202,7 @@ class VTK: self.geom.GetPointData().AddArray(d) elif isinstance(data,pd.DataFrame): pass - elif isinstance(data,table): + elif isinstance(data,Table): pass @@ -232,7 +234,8 @@ class VTK: ren.AddActor(actor) ren.SetBackground(0.2,0.2,0.2) - renWin.SetSize(1024, 1024) + + renWin.SetSize(Environment().screen_width,Environment().screen_height) iren = vtk.vtkRenderWindowInteractor() iren.SetRenderWindow(renWin) diff --git a/python/damask/result.py b/python/damask/result.py index 7d6a06441..748ea63fd 100644 --- a/python/damask/result.py +++ b/python/damask/result.py @@ -8,14 +8,15 @@ import h5py import numpy as np from . import VTK -from . import util -from . import version -from . import table -from . import mechanics +from . import Table from . import Rotation from . import Orientation from . import Environment from . import grid_filters +from . import mechanics +from . import util +from . import version + class Result: """ @@ -287,12 +288,12 @@ class Result: try: tbl[inc].add(path,data) except KeyError: - tbl[inc] = table.Table(data.reshape(self.Nmaterialpoints,-1),{path:data.shape[1:]}) + tbl[inc] = Table(data.reshape(self.Nmaterialpoints,-1),{path:data.shape[1:]}) else: try: tbl.add(path,data) except AttributeError: - tbl = table.Table(data.reshape(self.Nmaterialpoints,-1),{path:data.shape[1:]}) + tbl = Table(data.reshape(self.Nmaterialpoints,-1),{path:data.shape[1:]}) return tbl From 3d8e8cfe7a8e4aa85a8557331833f6feb95a1092 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 14 Mar 2020 00:07:49 +0100 Subject: [PATCH 19/19] tkinter might not work --- python/damask/environment.py | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/python/damask/environment.py b/python/damask/environment.py index b5feb4fcd..a4e4ed7b3 100644 --- a/python/damask/environment.py +++ b/python/damask/environment.py @@ -1,4 +1,3 @@ -import tkinter import os class Environment(): @@ -6,10 +5,14 @@ class Environment(): def __init__(self): """Read and provide values of DAMASK configuration.""" self.options = self._get_options() - tk = tkinter.Tk() - self.screen_width = tk.winfo_screenwidth() - self.screen_height = tk.winfo_screenheight() - + try: + import tkinter + tk = tkinter.Tk() + self.screen_width = tk.winfo_screenwidth() + self.screen_height = tk.winfo_screenheight() + except Exception: + self.screen_width = 1024 + self.screen_height = 768 def relPath(self,relative = '.'): """Return absolute path from path relative to DAMASK root."""