From c85d27baeced2c597d510fb4a4ad60132c367c72 Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Wed, 23 Sep 2020 09:45:36 +0200
Subject: [PATCH] functionality for seed creation

---
 python/damask/__init__.py |  6 +++
 python/damask/seeds.py    | 90 +++++++++++++++++++++++++++++++++++++++
 2 files changed, 96 insertions(+)
 create mode 100644 python/damask/seeds.py

diff --git a/python/damask/__init__.py b/python/damask/__init__.py
index 3f2ff6813..e21ebd03d 100644
--- a/python/damask/__init__.py
+++ b/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 = _
diff --git a/python/damask/seeds.py b/python/damask/seeds.py
new file mode 100644
index 000000000..e2cb5fc75
--- /dev/null
+++ b/python/damask/seeds.py
@@ -0,0 +1,90 @@
+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 periodic field.
+    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,_np.arange(N_seeds,dtype='i'))
+
+
+def from_poisson_disc(size,N_seeds,N_candidates,distance,seed=None):
+    """
+    Seeding in space according to a Poisson disc distribution.
+    
+    Parameters
+    ----------
+    size : numpy.ndarray of shape (3)
+        Physical size of the periodic field.
+    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.
+    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])
+        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,_np.arange(N_seeds,dtype='i'))
+
+
+def from_geom(geom,selection=None,invert=False):
+    """
+    Create seed from existing geometry description.
+    
+    Parameters
+    ----------
+    geom : damask.Geom
+        ...
+    
+    """
+    materials = geom.materials.reshape((-1,1),order='F')
+    mask = _np.isin(materials,selection,invert=invert) if selection is not None else \
+           _np.full(geom.grid.prod(),True,dtype=bool)
+    coords = grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
+
+    return (coords[mask],materials[mask])