From a59e64a8e459b779a439edb0fc4022084d1de1d5 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 30 Sep 2020 17:27:49 -0400 Subject: [PATCH] renamed TPMS and added more from additional references --- python/damask/_geom.py | 154 +++++++++++++++++++++++++------------- python/tests/test_Geom.py | 36 +++++++-- 2 files changed, 134 insertions(+), 56 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index aaa534df0..e760e40c3 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -120,6 +120,29 @@ class Geom: return np.unique(self.material).size + @staticmethod + def load(fname): + """ + Read a VTK rectilinear grid. + + Parameters + ---------- + fname : str or or pathlib.Path + Geometry file to read. + Valid extension is .vtr, it will be appended if not given. + + """ + v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr') + comments = v.get_comments() + grid = np.array(v.vtk_data.GetDimensions())-1 + bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T + + return Geom(material = v.get('material').reshape(grid,order='F'), + size = bbox[1] - bbox[0], + origin = bbox[0], + comments=comments) + + @staticmethod def load_ASCII(fname): """ @@ -184,29 +207,6 @@ class Geom: return Geom(material.reshape(grid,order='F'),size,origin,comments) - @staticmethod - def load(fname): - """ - Read a VTK rectilinear grid. - - Parameters - ---------- - fname : str or or pathlib.Path - Geometry file to read. - Valid extension is .vtr, it will be appended if not given. - - """ - v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr') - comments = v.get_comments() - grid = np.array(v.vtk_data.GetDimensions())-1 - bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T - - return Geom(material = v.get('material').reshape(grid,order='F'), - size = bbox[1] - bbox[0], - origin = bbox[0], - comments=comments) - - @staticmethod def _find_closest_seed(seeds, weights, point): return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights) @@ -292,10 +292,52 @@ class Geom: ) + _minimal_surface = \ + {'Schwarz P': lambda x,y,z: np.cos(x) + np.cos(y) + np.cos(z), + 'Double Primitive': lambda x,y,z: ( 0.5 * (np.cos(x)*np.cos(y) + np.cos(y)*np.cos(z) + np.cos(z)*np.cos(x)) + + 0.2 * (np.cos(2*x) + np.cos(2*y) + np.cos(2*z)) ), + 'Schwarz D': lambda x,y,z: ( np.sin(x)*np.sin(y)*np.sin(z) + + np.sin(x)*np.cos(y)*np.cos(z) + + np.cos(x)*np.cos(y)*np.sin(z) + + np.cos(x)*np.sin(y)*np.cos(z) ), + 'Complementary D': lambda x,y,z: ( np.cos(3*x+y)*np.cos(z) - np.sin(3*x-y)*np.sin(z) + np.cos(x+3*y)*np.cos(z) + + np.sin(x-3*y)*np.sin(z) + np.cos(x-y)*np.cos(3*z) - np.sin(x+y)*np.sin(3*z) ), + 'Double Diamond': lambda x,y,z: 0.5 * (np.sin(x)*np.sin(y) + + np.sin(y)*np.sin(z) + + np.sin(z)*np.sin(x) + + np.cos(x) * np.cos(y) * np.cos(z) ), + 'Dprime': lambda x,y,z: 0.5 * ( np.cos(x)*np.cos(y)*np.cos(z) + + np.cos(x)*np.sin(y)*np.sin(z) + + np.sin(x)*np.cos(y)*np.sin(z) + + np.sin(x)*np.sin(y)*np.cos(z) + - np.sin(2*x)*np.sin(2*y) + - np.sin(2*y)*np.sin(2*z) + - np.sin(2*z)*np.sin(2*x) ) - 0.2, + 'Gyroid': lambda x,y,z: np.cos(x)*np.sin(y) + np.cos(y)*np.sin(z) + np.cos(z)*np.sin(x), + 'Gprime': lambda x,y,z : ( np.sin(2*x)*np.cos(y)*np.sin(z) + + np.sin(2*y)*np.cos(z)*np.sin(x) + + np.sin(2*z)*np.cos(x)*np.sin(y) ) + 0.32, + 'Karcher K': lambda x,y,z: ( 0.3 * ( np.cos(x) + np.cos(y) + np.cos(z) + + np.cos(x)*np.cos(y) + np.cos(y)*np.cos(z) + np.cos(z)*np.cos(x) ) + - 0.4 * ( np.cos(2*x) + np.cos(2*y) + np.cos(2*z) ) ) + 0.2, + 'Lidinoid': lambda x,y,z: 0.5 * ( np.sin(2*x)*np.cos(y)*np.sin(z) + + np.sin(2*y)*np.cos(z)*np.sin(x) + + np.sin(2*z)*np.cos(x)*np.sin(y) + - np.cos(2*x)*np.cos(2*y) + - np.cos(2*y)*np.cos(2*z) + - np.cos(2*z)*np.cos(2*x) ) + 0.15, + 'Neovius': lambda x,y,z: ( 3 * (np.cos(x)+np.cos(y)+np.cos(z)) + + 4 * np.cos(x)*np.cos(y)*np.cos(z) ), + 'Fisher-Koch S': lambda x,y,z: ( np.cos(2*x)*np.sin( y)*np.cos( z) + + np.cos( x)*np.cos(2*y)*np.sin( z) + + np.sin( x)*np.cos( y)*np.cos(2*z) ), + } + + @staticmethod def from_minimal_surface(grid,size,surface,threshold=0.0,periods=1,materials=(1,2)): """ - Generate geometry from definition of minimal surface. + Generate geometry from definition of triply periodic minimal surface. Parameters ---------- @@ -303,7 +345,7 @@ class Geom: 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'} + surface : {'primitive', 'gyroid', 'lidinoid', 'neovius', 'diamond', 'doublediamond'} Type of the minimal surface. threshold : float, optional. Threshold of the minimal surface. Defaults to 0.0. @@ -312,21 +354,53 @@ class Geom: materials : (int, int), optional Material IDs. Defaults to (1,2). + References + ---------- + Surface curvature in triply-periodic minimal surface architectures as + a distinct design parameter in preparing advanced tissue engineering scaffolds + Sébastien B G Blanquer, Maike Werner, Markus Hannula, Shahriar Sharifi, + Guillaume P R Lajoinie, David Eglin, Jari Hyttinen, André A Poot, and Dirk W Grijpma + 10.1088/1758-5090/aa6553 + + Triply Periodic Bicontinuous Cubic Microdomain Morphologies by Symmetries + Meinhard Wohlgemuth, Nataliya Yufa, James Hoffman, and Edwin L. Thomas + 10.1021/ma0019499 + + Minisurf – A minimal surface generator for finite element modeling and additive manufacturing + Meng-Ting Hsieh, Lorenzo Valdevit + 10.1016/j.simpa.2020.100026 + """ - 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]), + return Geom(material = np.where(threshold < Geom._minimal_surface[surface](x,y,z),materials[1],materials[0]), size = size, comments = util.execution_stamp('Geom','from_minimal_surface'), ) + def save(self,fname,compress=True): + """ + Generates vtk rectilinear grid. + + Parameters + ---------- + fname : str, optional + Filename to write. If no file is given, a string is returned. + Valid extension is .vtr, it will be appended if not given. + compress : bool, optional + Compress with zlib algorithm. Defaults to True. + + """ + v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) + v.add(self.material.flatten(order='F'),'material') + v.add_comments(self.comments) + + v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress) + + def save_ASCII(self,fname,compress=None): """ Writes a geom file. @@ -399,26 +473,6 @@ class Geom: f.write(f'{reps} of {former}\n') - def save(self,fname,compress=True): - """ - Generates vtk rectilinear grid. - - Parameters - ---------- - fname : str, optional - Filename to write. If no file is given, a string is returned. - Valid extension is .vtr, it will be appended if not given. - compress : bool, optional - Compress with zlib algorithm. Defaults to True. - - """ - v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) - v.add(self.material.flatten(order='F'),'material') - v.add_comments(self.comments) - - v.save(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr',parallel=False,compress=compress) - - def show(self): """Show on screen.""" v = VTK.from_rectilinearGrid(self.grid,self.size,self.origin) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 27829ac34..884ae72e7 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -362,7 +362,19 @@ class TestGeom: assert np.all(geom.material == material) - @pytest.mark.parametrize('surface',['primitive','gyroid','diamond']) + @pytest.mark.parametrize('surface',['Schwarz P', + 'Double Primitive', + 'Schwarz D', + 'Complementary D', + 'Double Diamond', + 'Dprime', + 'Gyroid', + 'Gprime', + 'Karcher K', + 'Lidinoid', + 'Neovius', + 'Fisher-Koch S', + ]) def test_minimal_surface_basic_properties(self,surface): grid = np.random.randint(60,100,3) size = np.ones(3)+np.random.rand(3) @@ -373,8 +385,20 @@ class TestGeom: 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) + @pytest.mark.parametrize('surface,threshold',[('Schwarz P',0), + ('Double Primitive',-1./6.), + ('Schwarz D',0), + ('Complementary D',0), + ('Double Diamond',-0.133), + ('Dprime',-0.0395), + ('Gyroid',0), + ('Gprime',0.22913), + ('Karcher K',0.17045), + ('Lidinoid',0.14455), + ('Neovius',0), + ('Fisher-Koch S',0), + ]) + def test_minimal_surface_volume(self,surface,threshold): + 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)