diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index bb4a2e7a0..40662cb40 100644 --- a/.gitlab-ci.yml +++ b/.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: 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/_geom.py b/python/damask/_geom.py index 6af265021..50d4e2b33 100644 --- a/python/damask/_geom.py +++ b/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[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0] + 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] 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'), ) diff --git a/python/damask/seeds.py b/python/damask/seeds.py new file mode 100644 index 000000000..ea40794da --- /dev/null +++ b/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) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index ac4a8eac3..9aff60186 100644 --- a/python/tests/test_Geom.py +++ b/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) - Laguerre = Geom.from_Laguerre_tessellation(grid,size,seeds,np.ones(N_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),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) diff --git a/python/tests/test_seeds.py b/python/tests/test_seeds.py new file mode 100644 index 000000000..49d9a6bfd --- /dev/null +++ b/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=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()