diff --git a/processing/pre/geom_fromMinimalSurface.py b/processing/pre/geom_fromMinimalSurface.py index 83d1a3684..bf0e2b754 100755 --- a/processing/pre/geom_fromMinimalSurface.py +++ b/processing/pre/geom_fromMinimalSurface.py @@ -4,8 +4,6 @@ import os import sys from optparse import OptionParser -import numpy as np - import damask @@ -15,13 +13,6 @@ scriptID = ' '.join([scriptName,damask.version]) minimal_surfaces = ['primitive','gyroid','diamond'] -surface = { - 'primitive': lambda x,y,z: np.cos(x)+np.cos(y)+np.cos(z), - 'gyroid': lambda x,y,z: np.sin(x)*np.cos(y)+np.sin(y)*np.cos(z)+np.cos(x)*np.sin(z), - 'diamond': lambda x,y,z: np.cos(x-y)*np.cos(z)+np.sin(x+y)*np.sin(z), - } - - # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- @@ -71,16 +62,8 @@ parser.set_defaults(type = minimal_surfaces[0], name = None if filename == [] else filename[0] damask.util.report(scriptName,name) -x,y,z = np.meshgrid(options.periods*2.0*np.pi*(np.arange(options.grid[0])+0.5)/options.grid[0], - options.periods*2.0*np.pi*(np.arange(options.grid[1])+0.5)/options.grid[1], - options.periods*2.0*np.pi*(np.arange(options.grid[2])+0.5)/options.grid[2], - indexing='xy',sparse=True) - -microstructure = np.where(options.threshold < surface[options.type](x,y,z), - options.microstructure[1],options.microstructure[0]) - -geom=damask.Geom(microstructure,options.size, - comments=[scriptID + ' ' + ' '.join(sys.argv[1:])]) +geom=damask.Geom.from_minimal_surface(options.grid,options.size,options.type,options.threshold, + options.periods,options.microstructure) damask.util.croak(geom) geom.save_ASCII(sys.stdout if name is None else name,compress=False) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 6a7bc3c0e..aaa534df0 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -292,6 +292,41 @@ class Geom: ) + @staticmethod + def from_minimal_surface(grid,size,surface,threshold=0.0,periods=1,materials=(1,2)): + """ + Generate geometry from definition of minimal surface. + + Parameters + ---------- + grid : int numpy.ndarray of shape (3) + Number of grid points in x,y,z direction. + size : list or numpy.ndarray of shape (3) + Physical size of the geometry in meter. + surface : {'primitive', 'gyroid', 'diamond'} + Type of the minimal surface. + threshold : float, optional. + Threshold of the minimal surface. Defaults to 0.0. + periods : integer, optional. + Number of periods per unit cell. Defaults to 1. + materials : (int, int), optional + Material IDs. Defaults to (1,2). + + """ + s = {'primitive': lambda x,y,z: np.cos(x)+np.cos(y)+np.cos(z), + 'gyroid': lambda x,y,z: np.sin(x)*np.cos(y)+np.sin(y)*np.cos(z)+np.cos(x)*np.sin(z), + 'diamond': lambda x,y,z: np.cos(x-y)*np.cos(z)+np.sin(x+y)*np.sin(z), + } + x,y,z = np.meshgrid(periods*2.0*np.pi*(np.arange(grid[0])+0.5)/grid[0], + periods*2.0*np.pi*(np.arange(grid[1])+0.5)/grid[1], + periods*2.0*np.pi*(np.arange(grid[2])+0.5)/grid[2], + indexing='ij',sparse=True) + return Geom(material = np.where(threshold < s[surface](x,y,z),materials[1],materials[0]), + size = size, + comments = util.execution_stamp('Geom','from_minimal_surface'), + ) + + def save_ASCII(self,fname,compress=None): """ Writes a geom file. diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 21eb38cdd..27829ac34 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -360,3 +360,21 @@ class TestGeom: elif approach == 'Voronoi': geom = Geom.from_Voronoi_tessellation(grid,size,seeds, periodic=np.random.random()>0.5) assert np.all(geom.material == material) + + + @pytest.mark.parametrize('surface',['primitive','gyroid','diamond']) + def test_minimal_surface_basic_properties(self,surface): + grid = np.random.randint(60,100,3) + size = np.ones(3)+np.random.rand(3) + threshold = np.random.rand()-.5 + periods = np.random.randint(2)+1 + materials = np.random.randint(0,40,2) + geom = Geom.from_minimal_surface(grid,size,surface,threshold,periods,materials) + assert geom.material.max() == materials.max() and geom.material.min() == materials.min() \ + and (geom.size == size).all() and (geom.grid == grid).all() + + @pytest.mark.parametrize('surface',['primitive','gyroid','diamond']) + def test_minimal_surface_volume(self,surface): + grid = np.ones(3,dtype='i')*64 + geom = Geom.from_minimal_surface(grid,np.ones(3),surface) + assert np.isclose(np.count_nonzero(geom.material==1)/np.prod(geom.grid),.5)