Merge branch 'geom-improvements' into 'development'

Geom improvements

See merge request damask/DAMASK!264

PRIVATE

@@ -1 +1 @@
-Subproject commit 1e8c66897820468ab46958d995005e2b69204d0e
+Subproject commit 3112a4dbfa1e926c07b7f9443161239b8a7e85ca

examples/.gitignore

@@ -3,6 +3,3 @@
 *.xdmf
 *.sta
 *.vt*
-*.geom
-*.seeds
-postProc

python/damask/_configmaterial.py

@@ -2,7 +2,6 @@ import copy
 
 import numpy as np
 
-from . import grid_filters
 from . import Config
 from . import Lattice
 from . import Rotation
@@ -26,7 +25,7 @@ class ConfigMaterial(Config):
 
 
     @staticmethod
-    def from_table(table,coordinates=None,constituents={},**kwargs):
+    def from_table(table,constituents={},**kwargs):
         """
         Load from an ASCII table.
 
@@ -34,10 +33,6 @@ class ConfigMaterial(Config):
         ----------
         table : damask.Table
             Table that contains material information.
-        coordinates : str, optional
-            Label of spatial coordiates. Used for sorting and performing a
-            sanity check. Default to None, in which case no sorting or checking is
-            peformed.
         constituents : dict, optional
             Entries for 'constituents'. The key is the name and the value specifies
             the label of the data column in the table
@@ -54,7 +49,7 @@ class ConfigMaterial(Config):
             pos  pos  pos   qu   qu    qu    qu   phase    homog
         0    0    0    0  0.19  0.8   0.24 -0.51  Aluminum SX
         1    1    0    0  0.8   0.19  0.24 -0.51  Steel    SX
-        >>> cm.from_table(t,'pos',{'O':'qu','phase':'phase'},homogenization='homog')
+        >>> cm.from_table(t,{'O':'qu','phase':'phase'},homogenization='homog')
         material:
           - constituents:
               - O: [0.19, 0.8, 0.24, -0.51]
@@ -68,17 +63,12 @@ class ConfigMaterial(Config):
             homogenization: SX
 
         """
-        if coordinates is not None:
+        constituents_ = {k:table.get(v) for k,v in constituents.items()}
-            t = table.sort_by([f'{i}_{coordinates}' for i in range(3,0,-1)])
+        kwargs_       = {k:table.get(v) for k,v in kwargs.items()}
-            grid_filters.coord0_check(t.get(coordinates))
-        else:
-            t = table
-
-        constituents_ = {k:t.get(v) for k,v in constituents.items()}
-        kwargs_       = {k:t.get(v) for k,v in kwargs.items()}
 
         _,idx = np.unique(np.hstack(list({**constituents_,**kwargs_}.values())),return_index=True,axis=0)
 
+        idx = np.sort(idx)
         constituents_ = {k:v[idx].squeeze() for k,v in constituents_.items()}
         kwargs_       = {k:v[idx].squeeze() for k,v in kwargs_.items()}
 
@@ -229,19 +219,18 @@ class ConfigMaterial(Config):
         Parameters
         ----------
         constituents : dict
-            Entries for 'constituents'. The key is the name and the value specifies
+            Entries for 'constituents' as key-value pair.
-            the label of the data column in the table
         **kwargs
-            Keyword arguments where the key is the name and  the value specifies
+            Key-value pairs.
-            the label of the data column in the table
 
         Examples
         --------
         >>> import damask
-        >>> m = damask.ConfigMaterial()
         >>> O = damask.Rotation.from_random(3).as_quaternion()
         >>> phase = ['Aluminum','Steel','Aluminum']
-        >>> m.material_add(constituents={'phase':phase,'O':O},homogenization='SX')
+        >>> m = damask.ConfigMaterial().material_add(constituents={'phase':phase,'O':O},
+        ...                                          homogenization='SX')
+        >>> m
         material:
           - constituents:
               - O: [0.577764, -0.146299, -0.617669, 0.513010]
@@ -264,15 +253,14 @@ class ConfigMaterial(Config):
         for k,v in kwargs.items():
             if hasattr(v,'__len__') and not isinstance(v,str):
                 if len(v) != len(c):
-                    raise ValueError('len mismatch 1')
+                    raise ValueError('Cannot add entries of different length')
                 for i,vv in enumerate(v):
                     c[i][k] = [w.item() for w in vv] if isinstance(vv,np.ndarray) else vv.item()
             else:
                 for i in range(len(c)):
                     c[i][k] = v
         dup = copy.deepcopy(self)
-        if 'material' not in dup: dup['material'] = []
+        dup['material'] = dup['material'] + c if 'material' in dup else c
-        dup['material'] +=c
 
         return dup
 
@@ -288,7 +276,7 @@ class ConfigMaterial(Config):
             for k,v in kwargs.items():
                 if hasattr(v,'__len__') and not isinstance(v,str):
                     if len(v) != N_material:
-                        raise ValueError('len mismatch 2')
+                        raise ValueError('Cannot add entries of different length')
                     for i,vv in enumerate(np.array(v)):
                         m[i][0][k] = [w.item() for w in vv] if isinstance(vv,np.ndarray) else vv.item()
                 else:

python/damask/_geom.py

@@ -4,6 +4,7 @@ from functools import partial
 from os import path
 
 import numpy as np
+import pandas as pd
 import h5py
 from scipy import ndimage,spatial
 
@@ -254,21 +255,26 @@ class Geom:
         table : damask.Table
             Table that contains material information.
         coordinates : str
-            Label of the column containing the spatial coordinates.
+            Label of the column containing the vector of spatial coordinates.
+            Need to be ordered (1./x fast, 3./z slow).
         labels : str or list of str
             Label(s) of the columns containing the material definition.
             Each unique combintation of values results in a material.
 
         """
-        t = table.sort_by([f'{i}_{coordinates}' for i in range(3,0,-1)])
+        grid,size,origin = grid_filters.cell_coord0_gridSizeOrigin(table.get(coordinates))
-
-        grid,size,origin = grid_filters.cell_coord0_gridSizeOrigin(t.get(coordinates))
 
         labels_ = [labels] if isinstance(labels,str) else labels
-        _,unique_inverse = np.unique(np.hstack([t.get(l) for l in labels_]),return_inverse=True,axis=0)
+        unique,unique_inverse = np.unique(np.hstack([table.get(l) for l in labels_]),return_inverse=True,axis=0)
-        ma = unique_inverse.reshape(grid,order='F') + 1
+        if len(unique) == grid.prod():
+            ma = np.arange(grid.prod())
+        else:
+            from_ma = pd.unique(unique_inverse)
+            sort_idx = np.argsort(from_ma)
+            idx = np.searchsorted(from_ma,unique_inverse,sorter = sort_idx)
+            ma = np.arange(from_ma.size)[sort_idx][idx]
 
-        return Geom(ma,size,origin,util.execution_stamp('Geom','from_table'))
+        return Geom(ma.reshape(grid,order='F'),size,origin,util.execution_stamp('Geom','from_table'))
 
 
     @staticmethod
@@ -687,12 +693,10 @@ class Geom:
 
 
     def renumber(self):
-        """Renumber sorted material indices to 1,...,N."""
+        """Renumber sorted material indices to 0,...,N-1."""
-        renumbered = np.empty(self.grid,dtype=self.material.dtype)
+        _,renumbered = np.unique(self.material,return_inverse=True)
-        for i, oldID in enumerate(np.unique(self.material)):
-            renumbered = np.where(self.material == oldID, i+1, renumbered)
 
-        return Geom(material = renumbered,
+        return Geom(material = renumbered.reshape(self.grid),
                     size     = self.size,
                     origin   = self.origin,
                     comments = self.comments+[util.execution_stamp('Geom','renumber')],
@@ -783,11 +787,13 @@ class Geom:
             New material indices.
 
         """
-        substituted = self.material.copy()
+        def mp(entry,mapper):
-        for from_ms,to_ms in zip(from_material,to_material):
+            return mapper[entry] if entry in mapper else entry
-            substituted[self.material==from_ms] = to_ms
        
-        return Geom(material = substituted,
+        mp = np.vectorize(mp)
+        mapper = dict(zip(from_material,to_material))
+
+        return Geom(material = mp(self.material,mapper).reshape(self.grid),
                     size     = self.size,
                     origin   = self.origin,
                     comments = self.comments+[util.execution_stamp('Geom','substitute')],

python/tests/test_ConfigMaterial.py

@@ -1,8 +1,10 @@
 import os
 
 import pytest
+import numpy as np
 
 from damask import ConfigMaterial
+from damask import Table
 
 @pytest.fixture
 def reference_dir(reference_dir_base):
@@ -74,3 +76,14 @@ class TestConfigMaterial:
         material_config = ConfigMaterial.load(reference_dir/'material.yaml')
         new = material_config.material_rename_homogenization({'Taylor':'isostrain'})
         assert not new.is_complete
+
+    def test_from_table(self):
+        N = np.random.randint(3,10)
+        a = np.vstack((np.hstack((np.arange(N),np.arange(N)[::-1])),np.ones(N*2),np.zeros(N*2),np.ones(N*2))).T
+        t = Table(a,{'varying':2,'constant':2})
+        c = ConfigMaterial.from_table(t,constituents={'a':'varying','b':'1_constant'},c='2_constant')
+        assert len(c['material']) == N
+        for i,m in enumerate(c['material']):
+            c = m['constituents'][0]
+            assert m['c'] == 1 and c['b'] == 0 and c['a'] == [i,1]
+

python/tests/test_Geom.py

@@ -3,8 +3,10 @@ import numpy as np
 
 from damask import VTK
 from damask import Geom
+from damask import Table
 from damask import Rotation
 from damask import util
+from damask import grid_filters
 
 
 def geom_equal(a,b):
@@ -176,6 +178,7 @@ class TestGeom:
         material = default.material.copy()
         for m in np.unique(material):
             material[material==m] = material.max() + np.random.randint(1,30)
+        default.material -= 1
         modified = Geom(material,
                         default.size,
                         default.origin)
@@ -194,6 +197,13 @@ class TestGeom:
                           modified.substitute(np.arange(default.material.max())+1+offset,
                                               np.arange(default.material.max())+1))
 
+    def test_substitute_invariant(self,default):
+        f = np.unique(default.material.flatten())[:np.random.randint(1,default.material.max())]
+        t = np.random.permutation(f)
+        modified = default.substitute(f,t)
+        assert np.array_equiv(t,f) or (not geom_equal(modified,default))
+        assert geom_equal(default, modified.substitute(t,f))
+
 
     @pytest.mark.parametrize('axis_angle',[np.array([1,0,0,86.7]), np.array([0,1,0,90.4]), np.array([0,0,1,90]),
                                            np.array([1,0,0,175]),np.array([0,-1,0,178]),np.array([0,0,1,180])])
@@ -363,3 +373,14 @@ class TestGeom:
         grid = np.ones(3,dtype=int)*64
         geom = Geom.from_minimal_surface(grid,np.ones(3),surface,threshold)
         assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5,rtol=1e-3)
+
+
+    def test_from_table(self):
+        grid = np.random.randint(60,100,3)
+        size = np.ones(3)+np.random.rand(3)
+        coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
+        z=np.ones(grid.prod())
+        z[grid[:2].prod()*int(grid[2]/2):]=0
+        t = Table(np.column_stack((coords,z)),{'coords':3,'z':1})
+        g = Geom.from_table(t,'coords',['1_coords','z'])
+        assert g.N_materials == g.grid[0]*2 and (g.material[:,:,-1]-g.material[:,:,0] == grid[0]).all()