Merge branch 'export_DAMASK_to_DREAM3D' into 'development'

Export Dream3D output

See merge request damask/DAMASK!835

processing/legacy/DADF5toDREAM3D.py

@@ -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<int32_t>'
-        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<float>'
-        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<uint32_t>'
-            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)

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)
+        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,xdmf_dir.resolve())))/hdf5_name
+        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)
+        out_dir.mkdir(parents=True,exist_ok=True)
-        with util.open_text((xdmf_dir/hdf5_name).with_suffix('.xdmf'),'w') as f:
+        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)
+        out_dir = Path.cwd() if target_dir is None else Path(target_dir)
-        vtk_dir.mkdir(parents=True,exist_ok=True)
+        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<int32_t>')
+                    add_attribute(cell['EulerAngles'], 'ComponentDimensions', np.array([3],np.uint64))
+                    add_attribute(cell['EulerAngles'], 'ObjectType', 'DataArray<float>')
+
+                    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<uint32_t>')
+                    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,

python/tests/resources/Result/2phase_irregularGrid.dream3d

@@ -0,0 +1 @@
+../GeomGrid/2phase_irregularGrid.dream3d

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']: