diff --git a/processing/legacy/DADF5toDREAM3D.py b/processing/legacy/DADF5toDREAM3D.py deleted file mode 100755 index 95b62c7ce..000000000 --- a/processing/legacy/DADF5toDREAM3D.py +++ /dev/null @@ -1,138 +0,0 @@ -#!/usr/bin/env python3 - -import argparse -import os - -import h5py -import numpy as np - -import damask - -class AttributeManagerNullterm(h5py.AttributeManager): - """ - Attribute management for DREAM.3D hdf5 files. - - String attribute values are stored as fixed-length string with NULLTERM - - References - ---------- - https://stackoverflow.com/questions/38267076 - https://stackoverflow.com/questions/52750232 - - """ - - def create(self, name, data, shape=None, dtype=None): - if isinstance(data,str): - tid = h5py.h5t.C_S1.copy() - tid.set_size(len(data + ' ')) - super().create(name=name,data=data+' ',dtype = h5py.Datatype(tid)) - else: - super().create(name=name,data=data,shape=shape,dtype=dtype) - - -h5py._hl.attrs.AttributeManager = AttributeManagerNullterm # 'Monkey patch' - - -# -------------------------------------------------------------------- -# Crystal structure specifications -# -------------------------------------------------------------------- -Crystal_structures = {'fcc': 1, - 'bcc': 1, - 'hcp': 0, - 'bct': 7, - 'ort': 6} #TODO: is bct Tetragonal low/Tetragonal high? -Phase_types = {'Primary': 0} #further additions to these can be done by looking at 'Create Ensemble Info' filter - - -# -------------------------------------------------------------------- -# MAIN -# -------------------------------------------------------------------- -parser = argparse.ArgumentParser(description='Creating a file for DREAM3D from DAMASK data') -parser.add_argument('filenames', nargs='+', - help='DADF5 files') -parser.add_argument('-d','--dir', dest='dir',default='postProc',metavar='string', - help='name of subdirectory relative to the location of the DADF5 file to hold output') -parser.add_argument('--inc',nargs='+', - help='Increment for which DREAM3D to be used, eg. 25',type=int) - -options = parser.parse_args() - -for filename in options.filenames: - f = damask.Result(filename) - N_digits = int(np.floor(np.log10(int(f.increments[-1][3:]))))+1 - - f.pick('increments',options.inc) - for inc in damask.util.show_progress(f.iterate('increments'),len(f.selection['increments'])): - dirname = os.path.abspath(os.path.join(os.path.dirname(filename),options.dir)) - try: - os.mkdir(dirname) - except FileExistsError: - pass - - o = h5py.File(dirname + '/' + os.path.splitext(filename)[0] \ - + '_inc_{}.dream3D'.format(inc[3:].zfill(N_digits)),'w') - o.attrs['DADF5toDREAM3D'] = '1.0' - o.attrs['FileVersion'] = '7.0' - - for g in ['DataContainerBundles','Pipeline']: # empty groups (needed) - o.create_group(g) - - data_container_label = 'DataContainers/ImageDataContainer' - cell_data_label = data_container_label + '/CellData' - - # Phase information of DREAM.3D is constituent ID in DAMASK - o[cell_data_label + '/Phases'] = f.get_constituent_ID().reshape(tuple(f.grid)+(1,)) - DAMASK_quaternion = f.read_dataset(f.get_dataset_location('orientation')) - # Convert: DAMASK uses P = -1, DREAM.3D uses P = +1. Also change position of imagninary part - DREAM_3D_quaternion = np.hstack((-DAMASK_quaternion['x'],-DAMASK_quaternion['y'],-DAMASK_quaternion['z'], - DAMASK_quaternion['w'])).astype(np.float32) - o[cell_data_label + '/Quats'] = DREAM_3D_quaternion.reshape(tuple(f.grid)+(4,)) - - # Attributes to CellData group - o[cell_data_label].attrs['AttributeMatrixType'] = np.array([3],np.uint32) - o[cell_data_label].attrs['TupleDimensions'] = f.grid.astype(np.uint64) - - # Common Attributes for groups in CellData - for group in ['/Phases','/Quats']: - o[cell_data_label + group].attrs['DataArrayVersion'] = np.array([2],np.int32) - o[cell_data_label + group].attrs['Tuple Axis Dimensions'] = 'x={},y={},z={}'.format(*f.grid) - - o[cell_data_label + '/Phases'].attrs['ComponentDimensions'] = np.array([1],np.uint64) - o[cell_data_label + '/Phases'].attrs['ObjectType'] = 'DataArray' - o[cell_data_label + '/Phases'].attrs['TupleDimensions'] = f.grid.astype(np.uint64) - - o[cell_data_label + '/Quats'].attrs['ComponentDimensions'] = np.array([4],np.uint64) - o[cell_data_label + '/Quats'].attrs['ObjectType'] = 'DataArray' - o[cell_data_label + '/Quats'].attrs['TupleDimensions'] = f.grid.astype(np.uint64) - - # Create EnsembleAttributeMatrix - ensemble_label = data_container_label + '/EnsembleAttributeMatrix' - - # Data CrystalStructures - o[ensemble_label + '/CrystalStructures'] = np.uint32(np.array([999,\ - Crystal_structures[f.get_crystal_structure()]])).reshape(2,1) - o[ensemble_label + '/PhaseTypes'] = np.uint32(np.array([999,Phase_types['Primary']])).reshape(2,1) # ToDo - - # Attributes Ensemble Matrix - o[ensemble_label].attrs['AttributeMatrixType'] = np.array([11],np.uint32) - o[ensemble_label].attrs['TupleDimensions'] = np.array([2], np.uint64) - - # Attributes for data in Ensemble matrix - for group in ['CrystalStructures','PhaseTypes']: # 'PhaseName' not required MD: But would be nice to take the phase name mapping - o[ensemble_label+'/'+group].attrs['ComponentDimensions'] = np.array([1],np.uint64) - o[ensemble_label+'/'+group].attrs['Tuple Axis Dimensions'] = 'x=2' - o[ensemble_label+'/'+group].attrs['DataArrayVersion'] = np.array([2],np.int32) - o[ensemble_label+'/'+group].attrs['ObjectType'] = 'DataArray' - o[ensemble_label+'/'+group].attrs['TupleDimensions'] = np.array([2],np.uint64) - - geom_label = data_container_label + '/_SIMPL_GEOMETRY' - - o[geom_label + '/DIMENSIONS'] = np.int64(f.grid) - o[geom_label + '/ORIGIN'] = np.float32(np.zeros(3)) - o[geom_label + '/SPACING'] = np.float32(f.size) - - o[geom_label].attrs['GeometryName'] = 'ImageGeometry' - o[geom_label].attrs['GeometryTypeName'] = 'ImageGeometry' - o[geom_label].attrs['GeometryType'] = np.array([0],np.uint32) - o[geom_label].attrs['SpatialDimensionality'] = np.array([3],np.uint32) - o[geom_label].attrs['UnitDimensionality'] = np.array([3],np.uint32) diff --git a/python/damask/_result.py b/python/damask/_result.py index 6893af1bc..6ffbd0352 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -18,6 +18,7 @@ from numpy import ma import damask from . import VTK from . import Orientation +from . import Rotation from . import grid_filters from . import mechanics from . import tensor @@ -61,6 +62,7 @@ def _empty_like(dataset: np.ma.core.MaskedArray, fill_value = fill_float if dataset.dtype in np.sctypes['float'] else fill_int, mask = True) + class Result: """ Add data to and export data from a DADF5 (DAMASK HDF5) file. @@ -1758,8 +1760,8 @@ class Result: hdf5_name = self.fname.name hdf5_dir = self.fname.parent - xdmf_dir = Path.cwd() if target_dir is None else Path(target_dir) - hdf5_link = (hdf5_dir if absolute_path else Path(os.path.relpath(hdf5_dir,xdmf_dir.resolve())))/hdf5_name + out_dir = Path.cwd() if target_dir is None else Path(target_dir) + hdf5_link = (hdf5_dir if absolute_path else Path(os.path.relpath(hdf5_dir,out_dir.resolve())))/hdf5_name with h5py.File(self.fname,'r') as f: for inc in self.visible['increments']: @@ -1819,8 +1821,8 @@ class Result: np.prod(shape))} data_items[-1].text = f'{hdf5_link}:{name}' - xdmf_dir.mkdir(parents=True,exist_ok=True) - with util.open_text((xdmf_dir/hdf5_name).with_suffix('.xdmf'),'w') as f: + out_dir.mkdir(parents=True,exist_ok=True) + with util.open_text((out_dir/hdf5_name).with_suffix('.xdmf'),'w') as f: f.write(xml.dom.minidom.parseString(ET.tostring(xdmf).decode()).toprettyxml()) @@ -1884,8 +1886,8 @@ class Result: at_cell_ph,in_data_ph,at_cell_ho,in_data_ho = self._mappings() - vtk_dir = Path.cwd() if target_dir is None else Path(target_dir) - vtk_dir.mkdir(parents=True,exist_ok=True) + out_dir = Path.cwd() if target_dir is None else Path(target_dir) + out_dir.mkdir(parents=True,exist_ok=True) with h5py.File(self.fname,'r') as f: if self.version_minor >= 13: @@ -1926,9 +1928,133 @@ class Result: v = v.set(' / '.join(['/'.join([ty,field,label]),dataset.dtype.metadata['unit']]),dataset) - v.save(vtk_dir/f'{self.fname.stem}_inc{inc.split(prefix_inc)[-1].zfill(N_digits)}', + v.save(out_dir/f'{self.fname.stem}_inc{inc.split(prefix_inc)[-1].zfill(N_digits)}', parallel=parallel) + def export_DREAM3D(self, + q: str = 'O', + target_dir: Union[None, str, Path] = None): + """ + Export the visible components to DREAM3D compatible files. + + One DREAM3D file per visible increment is created. + The geometry is based on the undeformed configuration. + + Parameters + ---------- + q : str, optional + Name of the dataset containing the crystallographic orientation as quaternions. + Defaults to 'O'. + + target_dir : str or pathlib.Path, optional + Directory to save DREAM3D files. Will be created if non-existent. + + """ + def add_attribute(obj,name,data): + """DREAM.3D requires fixed length string.""" + if isinstance(data,str): + tid = h5py.h5t.C_S1.copy() + tid.set_size(len(data)+1) + obj.attrs.create(name,data,dtype=h5py.Datatype(tid)) + else: + obj.attrs.create(name,data) + + def create_and_open(obj,name): + obj.create_group(name) + return obj[name] + + if self.N_constituents != 1 or not self.structured: + raise TypeError('DREAM3D output requires structured grid with single constituent.') + + N_digits = int(np.floor(np.log10(max(1,self.incs[-1]))))+1 + + + at_cell_ph,in_data_ph,_,_ = self._mappings() + + out_dir = Path.cwd() if target_dir is None else Path(target_dir) + out_dir.mkdir(parents=True,exist_ok=True) + + with h5py.File(self.fname,'r') as f: + for inc in util.show_progress(self.visible['increments']): + for c in range(self.N_constituents): + crystal_structure = [999] + phase_name = ['Unknown Phase Type'] + cell_orientation = np.zeros((np.prod(self.cells),3),np.float32) + phase_ID = np.zeros((np.prod(self.cells)),dtype=np.int32) + count = 1 + for label in self.visible['phases']: + try: + data = _read(f['/'.join([inc,'phase',label,'mechanical',q])]) + lattice = data.dtype.metadata['lattice'] + # Map to DREAM.3D IDs + if lattice == 'hP': + crystal_structure.append(0) + elif lattice in ['cI','cF']: + crystal_structure.append(1) + elif lattice == 'tI': + crystal_structure.append(8) + + cell_orientation[at_cell_ph[c][label],:] = \ + Rotation(data[in_data_ph[c][label],:]).as_Euler_angles().astype(np.float32) + phase_ID[at_cell_ph[c][label]] = count + phase_name.append(label) + count +=1 + except KeyError: + pass + + + with h5py.File(f'{out_dir}/{self.fname.stem}_inc{inc.split(prefix_inc)[-1].zfill(N_digits)}.dream3d','w') as f_out: + add_attribute(f_out,'FileVersion','7.0') + + for g in ['DataContainerBundles','Pipeline']: # empty groups (needed) + f_out.create_group(g) + + data_container = create_and_open(f_out,'DataContainers/SyntheticVolumeDataContainer') + + cell = create_and_open(data_container,'CellData') + add_attribute(cell,'AttributeMatrixType',np.array([3],np.uint32)) + add_attribute(cell,'TupleDimensions', np.array(self.cells,np.uint64)) + + cell['Phases'] = np.reshape(phase_ID,tuple(np.flip(self.cells))+(1,)) + cell['EulerAngles'] = cell_orientation.reshape(tuple(np.flip(self.cells))+(3,)) + for dataset in ['Phases','EulerAngles']: + add_attribute(cell[dataset],'DataArrayVersion',np.array([2],np.int32)) + add_attribute(cell[dataset],'Tuple Axis Dimensions','x={},y={},z={}'.format(*np.array(self.cells))) + add_attribute(cell[dataset],'TupleDimensions', np.array(self.cells,np.uint64)) + add_attribute(cell['Phases'], 'ComponentDimensions', np.array([1],np.uint64)) + add_attribute(cell['Phases'], 'ObjectType', 'DataArray') + add_attribute(cell['EulerAngles'], 'ComponentDimensions', np.array([3],np.uint64)) + add_attribute(cell['EulerAngles'], 'ObjectType', 'DataArray') + + cell_ensemble = create_and_open(data_container,'CellEnsembleData') + + cell_ensemble['CrystalStructures'] = np.array(crystal_structure,np.uint32).reshape(-1,1) + cell_ensemble['PhaseTypes'] = np.array([999] + [0]*(len(crystal_structure)-1),np.uint32).reshape(-1,1) + tid = h5py.h5t.C_S1.copy() + tid.set_size(h5py.h5t.VARIABLE) + tid.set_cset(h5py.h5t.CSET_ASCII) + cell_ensemble.create_dataset(name='PhaseName',data = phase_name, dtype=h5py.Datatype(tid)) + + cell_ensemble.attrs['AttributeMatrixType'] = np.array([11],np.uint32) + cell_ensemble.attrs['TupleDimensions'] = np.array([len(self.phases) + 1], np.uint64) + for group in ['CrystalStructures','PhaseTypes','PhaseName']: + add_attribute(cell_ensemble[group], 'ComponentDimensions', np.array([1],np.uint64)) + add_attribute(cell_ensemble[group], 'Tuple Axis Dimensions', f'x={len(self.phases)+1}') + add_attribute(cell_ensemble[group], 'DataArrayVersion', np.array([2],np.int32)) + add_attribute(cell_ensemble[group], 'TupleDimensions', np.array([len(self.phases) + 1],np.uint64)) + for group in ['CrystalStructures','PhaseTypes']: + add_attribute(cell_ensemble[group], 'ObjectType', 'DataArray') + add_attribute(cell_ensemble['PhaseName'], 'ObjectType', 'StringDataArray') + + geom = create_and_open(data_container,'_SIMPL_GEOMETRY') + geom['DIMENSIONS'] = np.array(self.cells,np.int64) + geom['ORIGIN'] = np.array(self.origin,np.float32) + geom['SPACING'] = np.float32(self.size/self.cells) + names = ['GeometryName', 'GeometryTypeName','GeometryType','SpatialDimensionality','UnitDimensionality'] + values = ['ImageGeometry','ImageGeometry', np.array([0],np.uint32)] + [np.array([3],np.uint32)]*2 + for name,value in zip(names,values): + add_attribute(geom,name,value) + def export_DADF5(self, fname, diff --git a/python/tests/resources/Result/2phase_irregularGrid.dream3d b/python/tests/resources/Result/2phase_irregularGrid.dream3d new file mode 120000 index 000000000..0c1611d87 --- /dev/null +++ b/python/tests/resources/Result/2phase_irregularGrid.dream3d @@ -0,0 +1 @@ +../GeomGrid/2phase_irregularGrid.dream3d \ No newline at end of file diff --git a/python/tests/resources/Result/2phase_irregularGrid_tensionX_2phase_irregularGrid.material.hdf5 b/python/tests/resources/Result/2phase_irregularGrid_tensionX_2phase_irregularGrid.material.hdf5 deleted file mode 100644 index 8370f9d6a..000000000 Binary files a/python/tests/resources/Result/2phase_irregularGrid_tensionX_2phase_irregularGrid.material.hdf5 and /dev/null differ diff --git a/python/tests/resources/Result/2phase_irregularGrid_tensionX_increment_0.dream3d b/python/tests/resources/Result/2phase_irregularGrid_tensionX_increment_0.dream3d deleted file mode 100644 index 933856445..000000000 Binary files a/python/tests/resources/Result/2phase_irregularGrid_tensionX_increment_0.dream3d and /dev/null differ diff --git a/python/tests/resources/Result/2phase_irregularGrid_tensionX_increment_40.dream3d b/python/tests/resources/Result/2phase_irregularGrid_tensionX_increment_40.dream3d deleted file mode 100644 index 42ec5da7f..000000000 Binary files a/python/tests/resources/Result/2phase_irregularGrid_tensionX_increment_40.dream3d and /dev/null differ diff --git a/python/tests/resources/Result/2phase_irregularGrid_tensionX_material.hdf5 b/python/tests/resources/Result/2phase_irregularGrid_tensionX_material.hdf5 new file mode 100644 index 000000000..83a2024d9 Binary files /dev/null and b/python/tests/resources/Result/2phase_irregularGrid_tensionX_material.hdf5 differ diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 8b92238ae..bad55c98d 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -58,6 +58,19 @@ def dict_equal(d1, d2): return False return True +@pytest.fixture +def h5py_dataset_iterator(): + """Iterate over all datasets in an HDF5 file.""" + def _h5py_dataset_iterator(g, prefix=''): + for key,item in g.items(): + path = os.path.join(prefix, key) + if isinstance(item, h5py.Dataset): # test for dataset + yield (path, item) + elif isinstance(item, h5py.Group): # test for group (go down) + yield from _h5py_dataset_iterator(item, path) + return _h5py_dataset_iterator + + class TestResult: def test_self_report(self,default): @@ -437,6 +450,35 @@ class TestResult: single_phase.export_VTK(mode='point',target_dir=export_dir,parallel=False) assert set(os.listdir(export_dir)) == set([f'{single_phase.fname.stem}_inc{i:02}.vtp' for i in range(0,40+1,4)]) + def test_export_DREAM3D(self,tmp_path,res_path,h5py_dataset_iterator): + result = Result(res_path/'2phase_irregularGrid_tensionX_material.hdf5').view(increments=0) # compare the initial data only + result.export_DREAM3D(target_dir=tmp_path) + + def ignore(path): + # features present in reference but not in exported file + for i in ['Pipeline','StatsGeneratorDataContainer','Grain Data', + 'BoundaryCells','FeatureIds','IPFColor','NumFeatures']: + if path.find(i) >= 0: return True + return False + + with h5py.File(res_path/'2phase_irregularGrid.dream3d','r') as ref, \ + h5py.File(tmp_path/'2phase_irregularGrid_tensionX_material_inc0.dream3d','r') as cur: + + for (path,dset) in h5py_dataset_iterator(ref): + if ignore(path): continue + if path.find('PhaseName') < 0: + assert np.array_equal(dset,cur[path]) + else: + c = [_.decode() for _ in cur[path]] + r = ['Unknown Phase Type'] + result.phases + assert c == r + grp = os.path.split(path)[0] + for attr in ref[grp].attrs: + assert np.array_equal(ref[grp].attrs[attr],cur[grp].attrs[attr]) + for attr in dset.attrs: + assert np.array_equal(dset.attrs[attr],cur[path].attrs[attr]) + + def test_XDMF_datatypes(self,tmp_path,single_phase,update,res_path): for what,shape in {'scalar':(),'vector':(3,),'tensor':(3,3),'matrix':(12,)}.items(): for dtype in ['f4','f8','i1','i2','i4','i8','u1','u2','u4','u8']: