Merge branch 'seeds-module' into 'development'

Seeds module

See merge request damask/DAMASK!236

.gitlab-ci.yml

@@ -344,6 +344,7 @@ Phenopowerlaw_singleSlip:
 Pytest_grid:
   stage: grid
   script:
+    - module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
     - cd pytest
     - pytest
   except:

python/damask/__init__.py

@@ -20,6 +20,12 @@ from ._result      import Result           # noqa
 from ._geom        import Geom             # noqa
 from ._material    import Material         # noqa
 from .             import solver           # noqa
+from .             import util             # noqa
+from .             import seeds            # noqa
+from .             import grid_filters     # noqa
+from .             import mechanics        # noqa
+
+
 
 # deprecated
 Environment = _

python/damask/_geom.py

@@ -212,7 +212,7 @@ class Geom:
         return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
 
     @staticmethod
-    def from_Laguerre_tessellation(grid,size,seeds,weights,periodic=True):
+    def from_Laguerre_tessellation(grid,size,seeds,weights,material=None,periodic=True):
         """
         Generate geometry from Laguerre tessellation.
 
@@ -226,6 +226,9 @@ class Geom:
             Position of the seed points in meter. All points need to lay within the box.
         weights : numpy.ndarray of shape (seeds.shape[0])
             Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
+        material : numpy.ndarray of shape (seeds.shape[0]), optional
+            Material ID of the seeds. Defaults to None, in which case materials are
+            consecutively numbered.
         periodic : Boolean, optional
             Perform a periodic tessellation. Defaults to True.
 
@@ -245,22 +248,22 @@ class Geom:
         result = pool.map_async(partial(Geom._find_closest_seed,seeds_p,weights_p), [coord for coord in coords])
         pool.close()
         pool.join()
-        material = np.array(result.get())
+        material_ = np.array(result.get())
 
         if periodic:
-            material = material.reshape(grid*3)
+            material_ = material_.reshape(grid*3)
-            material = material[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
+            material_ = material_[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
         else:
-            material = material.reshape(grid)
+            material_ = material_.reshape(grid)
 
-        return Geom(material = material+1,
+        return Geom(material = material_+1 if material is None else material[material_],
                     size     = size,
                     comments = util.execution_stamp('Geom','from_Laguerre_tessellation'),
                    )
 
 
     @staticmethod
-    def from_Voronoi_tessellation(grid,size,seeds,periodic=True):
+    def from_Voronoi_tessellation(grid,size,seeds,material=None,periodic=True):
         """
         Generate geometry from Voronoi tessellation.
 
@@ -272,15 +275,18 @@ class Geom:
             Physical size of the geometry in meter.
         seeds : numpy.ndarray of shape (:,3)
             Position of the seed points in meter. All points need to lay within the box.
+        material : numpy.ndarray of shape (seeds.shape[0]), optional
+            Material ID of the seeds. Defaults to None, in which case materials are
+            consecutively numbered.
         periodic : Boolean, optional
             Perform a periodic tessellation. Defaults to True.
 
         """
         coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
         KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
-        devNull,material = KDTree.query(coords)
+        devNull,material_ = KDTree.query(coords)
 
-        return Geom(material = material.reshape(grid)+1,
+        return Geom(material = (material_+1 if material is None else material[material_]).reshape(grid),
                     size     = size,
                     comments = util.execution_stamp('Geom','from_Voronoi_tessellation'),
                    )

python/damask/seeds.py

@@ -0,0 +1,113 @@
+from scipy import spatial as _spatial
+import numpy as _np
+
+from . import util
+from . import grid_filters
+
+
+def from_random(size,N_seeds,grid=None,seed=None):
+    """
+    Random seeding in space.
+
+    Parameters
+    ----------
+    size : numpy.ndarray of shape (3)
+        Physical size of the seeding domain.
+    N_seeds : int
+        Number of seeds.
+    grid : numpy.ndarray of shape (3), optional.
+        If given, ensures that all seeds initiate one grain if using a
+        standard Voronoi tessellation.
+    seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
+        A seed to initialize the BitGenerator. Defaults to None.
+        If None, then fresh, unpredictable entropy will be pulled from the OS.
+
+    """
+    rng = _np.random.default_rng(seed)
+    if grid is None:
+        coords = rng.random((N_seeds,3)) * size
+    else:
+        grid_coords = grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
+        coords = grid_coords[rng.choice(_np.prod(grid),N_seeds, replace=False)] \
+               + _np.broadcast_to(size/grid,(N_seeds,3))*(rng.random((N_seeds,3))*.5-.25)           # wobble without leaving grid
+
+    return coords
+
+
+def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,seed=None):
+    """
+    Seeding in space according to a Poisson disc distribution.
+
+    Parameters
+    ----------
+    size : numpy.ndarray of shape (3)
+        Physical size of the seeding domain.
+    N_seeds : int
+        Number of seeds.
+    N_candidates : int
+        Number of candidates to consider for finding best candidate.
+    distance : float
+        Minimum acceptable distance to other seeds.
+    periodic : boolean, optional
+        Calculate minimum distance for periodically repeated grid.
+    seed : {None, int, array_like[ints], SeedSequence, BitGenerator, Generator}, optional
+        A seed to initialize the BitGenerator. Defaults to None.
+        If None, then fresh, unpredictable entropy will be pulled from the OS.
+
+    """
+    rng = _np.random.default_rng(seed)
+    coords = _np.empty((N_seeds,3))
+    coords[0] = rng.random(3) * size
+
+    i = 1
+    progress = util._ProgressBar(N_seeds+1,'',50)
+    while i < N_seeds:
+        candidates = rng.random((N_candidates,3))*_np.broadcast_to(size,(N_candidates,3))
+        tree = _spatial.cKDTree(coords[:i],boxsize=size) if periodic else \
+               _spatial.cKDTree(coords[:i])
+        distances, dev_null = tree.query(candidates)
+        best = distances.argmax()
+        if distances[best] > distance:                                                              # require minimum separation
+            coords[i] = candidates[best]                                                            # maximum separation to existing point cloud
+            i += 1
+            progress.update(i)
+
+    return coords
+
+
+def from_geom(geom,selection=None,invert=False,average=False,periodic=True):
+    """
+    Create seed from existing geometry description.
+
+    Parameters
+    ----------
+    geom : damask.Geom
+        Geometry, from which the material IDs are used as seeds.
+    selection : iterable of integers, optional
+        Material IDs to consider.
+    invert : boolean, false
+        Do not consider the material IDs given in selection. Defaults to False.
+    average : boolean, optional
+        Seed corresponds to center of gravity of material ID cloud.
+    periodic : boolean, optional
+        Center of gravity with periodic boundaries.
+
+    """
+    material = geom.material.reshape((-1,1),order='F')
+    mask = _np.full(geom.grid.prod(),True,dtype=bool) if selection is None else \
+           _np.isin(material,selection,invert=invert)
+    coords = grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
+
+    if not average:
+        return (coords[mask],material[mask])
+    else:
+        materials = _np.unique(material[mask])
+        coords_ = _np.zeros((materials.size,3),dtype=float)
+        for i,mat in enumerate(materials):
+            pc = (2*_np.pi*coords[material[:,0]==mat,:]-geom.origin)/geom.size
+            coords_[i] = geom.origin + geom.size / 2 / _np.pi * (_np.pi +
+                         _np.arctan2(-_np.average(_np.sin(pc),axis=0),
+                                     -_np.average(_np.cos(pc),axis=0))) \
+                         if periodic else \
+                         _np.average(coords[material[:,0]==mat,:],axis=0)
+        return (coords_,materials)

python/tests/test_Geom.py

@@ -83,6 +83,13 @@ class TestGeom:
         with pytest.raises(ValueError):
             Geom.load(tmpdir/'no_materialpoint.vtr')
 
+    def test_invalid_material(self):
+        with pytest.raises(TypeError):
+            Geom(np.zeros((3,3,3),dtype='complex'),np.ones(3))
+
+    def test_cast_to_int(self):
+        g = Geom(np.zeros((3,3,3)),np.ones(3))
+        assert g.material.dtype in np.sctypes['int']
 
     @pytest.mark.parametrize('compress',[True,False])
     def test_compress(self,default,tmpdir,compress):
@@ -324,8 +331,8 @@ class TestGeom:
         size   = np.random.random(3) + 1.0
         N_seeds= np.random.randint(10,30)
         seeds  = np.random.rand(N_seeds,3) * np.broadcast_to(size,(N_seeds,3))
-        Voronoi  = Geom.from_Voronoi_tessellation( grid,size,seeds,                 periodic)
+        Voronoi  = Geom.from_Voronoi_tessellation( grid,size,seeds,                 np.arange(N_seeds)+5,periodic)
-        Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),periodic)
+        Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_seeds),np.arange(N_seeds)+5,periodic)
         assert geom_equal(Laguerre,Voronoi)
 
 
@@ -337,7 +344,7 @@ class TestGeom:
         weights= np.full((N_seeds),-np.inf)
         ms     = np.random.randint(1, N_seeds+1)
         weights[ms-1] = np.random.random()
-        Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,np.random.random()>0.5)
+        Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,weights,periodic=np.random.random()>0.5)
         assert np.all(Laguerre.material == ms)
 
 
@@ -349,7 +356,7 @@ class TestGeom:
         material = np.ones(grid)
         material[:,grid[1]//2:,:] = 2
         if   approach == 'Laguerre':
-            geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),np.random.random()>0.5)
+            geom = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(2),periodic=np.random.random()>0.5)
         elif approach == 'Voronoi':
-            geom = Geom.from_Voronoi_tessellation(grid,size,seeds,            np.random.random()>0.5)
+            geom = Geom.from_Voronoi_tessellation(grid,size,seeds,            periodic=np.random.random()>0.5)
         assert np.all(geom.material == material)

python/tests/test_seeds.py

@@ -0,0 +1,36 @@
+import pytest
+import numpy as np
+from scipy.spatial import cKDTree
+
+from damask import seeds
+from damask import Geom
+
+class TestSeeds:
+
+    @pytest.mark.parametrize('grid',[None,np.ones(3,dtype='i')*10])
+    def test_from_random(self,grid):
+        N_seeds = np.random.randint(30,300)
+        size = np.ones(3) + np.random.random(3)
+        coords = seeds.from_random(size,N_seeds,grid)
+        assert (0<=coords).all() and (coords<size).all()
+
+    @pytest.mark.parametrize('periodic',[True,False])
+    def test_from_Poisson_disc(self,periodic):
+        N_seeds = np.random.randint(30,300)
+        N_candidates = N_seeds//15
+        distance = np.random.random()
+        size = np.ones(3)*distance*N_seeds
+        coords = seeds.from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=periodic)
+        min_dists, _ = cKDTree(coords,boxsize=size).query(coords, 2) if periodic else \
+                       cKDTree(coords).query(coords, 2)
+        assert (0<= coords).all() and (coords<size).all() and np.min(min_dists[:,1])>=distance
+
+    def test_from_geom(self):
+        grid = np.random.randint(10,20,3)
+        N_seeds = np.random.randint(30,300)
+        size = np.ones(3) + np.random.random(3)
+        coords = seeds.from_random(size,N_seeds,grid)
+        geom_1 = Geom.from_Voronoi_tessellation(grid,size,coords)
+        coords,material = seeds.from_geom(geom_1)
+        geom_2 = Geom.from_Voronoi_tessellation(grid,size,coords,material)
+        assert (geom_2.material==geom_1.material).all()