DAMASK_EICMD/python/damask/_grid.py

1208 lines
47 KiB
Python
Raw Normal View History

2021-06-15 23:08:01 +05:30
import os
2020-08-22 23:25:18 +05:30
import copy
2021-06-15 23:08:01 +05:30
import warnings
2020-09-18 20:02:08 +05:30
import multiprocessing as mp
2020-03-29 22:42:23 +05:30
from functools import partial
2022-01-13 01:04:29 +05:30
import typing
2022-01-30 20:16:26 +05:30
from typing import Union, Optional, TextIO, List, Sequence
2021-12-06 18:52:52 +05:30
from pathlib import Path
2019-05-26 15:33:21 +05:30
import numpy as np
import pandas as pd
import h5py
2020-12-11 01:23:11 +05:30
from scipy import ndimage, spatial
2019-05-26 15:33:21 +05:30
2020-03-11 12:02:03 +05:30
from . import VTK
2019-05-28 01:30:26 +05:30
from . import util
2020-03-29 22:42:23 +05:30
from . import grid_filters
from . import Rotation
2021-12-06 18:52:52 +05:30
from . import Table
2022-01-13 01:04:29 +05:30
from ._typehints import FloatSequence, IntSequence
2019-05-26 15:33:21 +05:30
2020-12-04 11:42:18 +05:30
class Grid:
"""
Geometry definition for grid solvers.
Create and manipulate geometry definitions for storage as VTK
2021-06-15 23:08:01 +05:30
image data files ('.vti' extension). A grid contains the
material ID (referring to the entry in 'material.yaml') and
the physical size.
"""
2019-11-23 01:22:36 +05:30
2021-12-06 18:52:52 +05:30
def __init__(self,
material: np.ndarray,
2022-01-13 01:04:29 +05:30
size: FloatSequence,
origin: FloatSequence = np.zeros(3),
comments: Union[str, Sequence[str]] = None):
2019-11-23 01:22:36 +05:30
"""
New geometry definition for grid solvers.
2019-11-23 01:22:36 +05:30
Parameters
----------
2022-01-13 03:43:38 +05:30
material : numpy.ndarray, shape (:,:,:)
Material indices. The shape of the material array defines
the number of cells.
2022-01-13 01:04:29 +05:30
size : sequence of float, len (3)
2021-06-17 21:56:37 +05:30
Physical size of grid in meter.
2022-01-13 01:04:29 +05:30
origin : sequence of float, len (3), optional
Coordinates of grid origin in meter. Defaults to [0.0,0.0,0.0].
comments : (list of) str, optional
2021-06-17 21:56:37 +05:30
Comments, e.g. history of operations.
2019-11-23 01:22:36 +05:30
"""
self.material = material
2022-01-13 01:04:29 +05:30
self.size = size # type: ignore
self.origin = origin # type: ignore
self.comments = [] if comments is None else comments # type: ignore
2020-03-21 15:37:21 +05:30
2021-12-06 18:52:52 +05:30
def __repr__(self) -> str:
2020-12-04 11:42:18 +05:30
"""Basic information on grid definition."""
mat_min = np.nanmin(self.material)
mat_max = np.nanmax(self.material)
mat_N = self.N_materials
2019-11-23 01:22:36 +05:30
return util.srepr([
f'cells : {util.srepr(self.cells, " x ")}',
2021-11-16 20:58:23 +05:30
f'size : {util.srepr(self.size, " x ")} / m³',
f'origin: {util.srepr(self.origin," ")} / m',
2020-12-11 01:23:11 +05:30
f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else
f' (min: {mat_min}, max: {mat_max})')
2019-11-23 01:22:36 +05:30
])
2022-01-22 12:35:49 +05:30
def __copy__(self) -> 'Grid':
2021-01-03 16:33:40 +05:30
"""Create deep copy."""
2020-08-22 23:25:18 +05:30
return copy.deepcopy(self)
2021-01-03 16:33:40 +05:30
copy = __copy__
2020-08-22 23:25:18 +05:30
2022-01-27 04:07:07 +05:30
def __eq__(self,
other: object) -> bool:
"""
Test equality of other.
Parameters
----------
2020-12-04 11:42:18 +05:30
other : damask.Grid
Grid to compare self against.
"""
2021-12-06 18:52:52 +05:30
if not isinstance(other, Grid):
2022-01-13 01:04:29 +05:30
return NotImplemented
return bool(np.allclose(other.size,self.size)
and np.allclose(other.origin,self.origin)
and np.all(other.cells == self.cells)
and np.all(other.material == self.material))
2019-11-23 01:22:36 +05:30
@property
2021-12-06 18:52:52 +05:30
def material(self) -> np.ndarray:
"""Material indices."""
return self._material
@material.setter
def material(self,
material: np.ndarray):
if len(material.shape) != 3:
2020-12-11 01:23:11 +05:30
raise ValueError(f'invalid material shape {material.shape}')
if material.dtype not in np.sctypes['float'] and material.dtype not in np.sctypes['int']:
2020-12-11 01:23:11 +05:30
raise TypeError(f'invalid material data type {material.dtype}')
self._material = np.copy(material)
if self.material.dtype in np.sctypes['float'] and \
np.all(self.material == self.material.astype(int).astype(float)):
self._material = self.material.astype(int)
@property
2021-12-06 18:52:52 +05:30
def size(self) -> np.ndarray:
2020-12-04 11:42:18 +05:30
"""Physical size of grid in meter."""
return self._size
@size.setter
def size(self,
size: FloatSequence):
if len(size) != 3 or any(np.array(size) < 0):
2020-12-11 01:23:11 +05:30
raise ValueError(f'invalid size {size}')
self._size = np.array(size)
@property
2022-01-13 01:04:29 +05:30
def origin(self) -> np.ndarray:
2020-12-04 11:42:18 +05:30
"""Coordinates of grid origin in meter."""
return self._origin
@origin.setter
def origin(self,
origin: FloatSequence):
if len(origin) != 3:
2020-12-11 01:23:11 +05:30
raise ValueError(f'invalid origin {origin}')
self._origin = np.array(origin)
@property
2021-12-06 18:52:52 +05:30
def comments(self) -> List[str]:
2020-12-04 11:42:18 +05:30
"""Comments, e.g. history of operations."""
return self._comments
@comments.setter
def comments(self,
comments: Union[str, Sequence[str]]):
self._comments = [str(c) for c in comments] if isinstance(comments,list) else [str(comments)]
2019-12-08 13:47:57 +05:30
@property
2021-12-06 18:52:52 +05:30
def cells(self) -> np.ndarray:
2020-12-04 02:28:24 +05:30
"""Number of cells in x,y,z direction."""
return np.asarray(self.material.shape)
2019-11-23 01:22:36 +05:30
2020-03-18 18:19:53 +05:30
@property
2021-12-06 18:52:52 +05:30
def N_materials(self) -> int:
2020-12-04 11:42:18 +05:30
"""Number of (unique) material indices within grid."""
return np.unique(self.material).size
2019-11-23 01:22:36 +05:30
@staticmethod
2022-01-22 12:35:49 +05:30
def load(fname: Union[str, Path]) -> 'Grid':
"""
Load from VTK image data file.
Parameters
----------
2022-01-13 01:04:29 +05:30
fname : str or pathlib.Path
Grid file to read. Valid extension is .vti, which will be appended
2020-12-04 11:42:18 +05:30
if not given.
Returns
-------
loaded : damask.Grid
2021-04-24 18:17:52 +05:30
Grid-based geometry from file.
"""
2021-06-16 18:05:54 +05:30
v = VTK.load(fname if str(fname).endswith(('.vti','.vtr')) else str(fname)+'.vti') # compatibility hack
comments = v.get_comments()
2020-12-04 02:28:24 +05:30
cells = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
2020-12-04 11:42:18 +05:30
return Grid(material = v.get('material').reshape(cells,order='F'),
size = bbox[1] - bbox[0],
origin = bbox[0],
comments=comments)
2022-01-13 01:04:29 +05:30
@typing. no_type_check
@staticmethod
2022-01-22 12:35:49 +05:30
def load_ASCII(fname)-> 'Grid':
2019-11-23 01:22:36 +05:30
"""
2020-12-04 02:28:24 +05:30
Load from geom file.
2019-11-23 01:22:36 +05:30
2020-11-14 22:24:47 +05:30
Storing geometry files in ASCII format is deprecated.
This function will be removed in a future version of DAMASK.
2019-11-23 01:22:36 +05:30
Parameters
----------
fname : str, pathlib.Path, or file handle
2020-08-08 23:12:34 +05:30
Geometry file to read.
2019-11-23 01:22:36 +05:30
Returns
-------
loaded : damask.Grid
2021-04-24 18:17:52 +05:30
Grid-based geometry from file.
2019-11-23 01:22:36 +05:30
"""
2022-01-05 18:36:06 +05:30
warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.0.0', DeprecationWarning,2)
2021-12-06 18:52:52 +05:30
if isinstance(fname, (str, Path)):
f = open(fname)
2021-12-06 18:52:52 +05:30
elif isinstance(fname, TextIO):
f = fname
2021-12-06 18:52:52 +05:30
else:
raise TypeError
f.seek(0)
try:
2021-12-06 18:52:52 +05:30
header_length_,keyword = f.readline().split()[:2]
header_length = int(header_length_)
except ValueError:
header_length,keyword = (-1, 'invalid')
2019-11-23 01:22:36 +05:30
if not keyword.startswith('head') or header_length < 3:
2020-12-11 01:23:11 +05:30
raise TypeError('header length information missing or invalid')
2019-11-23 01:22:36 +05:30
2020-03-21 15:37:21 +05:30
comments = []
2020-11-20 00:56:15 +05:30
content = f.readlines()
2019-11-23 01:22:36 +05:30
for i,line in enumerate(content[:header_length]):
2022-01-13 01:04:29 +05:30
items = line.split('#')[0].lower().strip().split()
if (key := items[0] if items else '') == 'grid':
2022-01-13 01:04:29 +05:30
cells = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']])
2019-11-23 01:22:36 +05:30
elif key == 'size':
size = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
elif key == 'origin':
origin = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
else:
comments.append(line.strip())
2022-01-13 01:04:29 +05:30
material = np.empty(int(cells.prod())) # initialize as flat array
2019-11-23 01:22:36 +05:30
i = 0
for line in content[header_length:]:
if len(items := line.split('#')[0].split()) == 3:
2021-12-06 18:52:52 +05:30
if items[1].lower() == 'of':
material_entry = np.ones(int(items[0]))*float(items[2])
2019-11-23 01:22:36 +05:30
elif items[1].lower() == 'to':
2021-12-06 18:52:52 +05:30
material_entry = np.linspace(int(items[0]),int(items[2]),
2019-11-23 01:22:36 +05:30
abs(int(items[2])-int(items[0]))+1,dtype=float)
2021-12-06 18:52:52 +05:30
else: material_entry = list(map(float, items))
else: material_entry = list(map(float, items))
material[i:i+len(material_entry)] = material_entry
2019-11-23 01:22:36 +05:30
i += len(items)
2020-03-21 15:37:21 +05:30
2020-12-04 02:28:24 +05:30
if i != cells.prod():
2020-12-11 01:23:11 +05:30
raise TypeError(f'invalid file: expected {cells.prod()} entries, found {i}')
2020-03-21 15:37:21 +05:30
if not np.any(np.mod(material,1) != 0.0): # no float present
material = material.astype('int') - (1 if material.min() > 0 else 0)
2020-03-21 15:37:21 +05:30
2020-12-04 11:42:18 +05:30
return Grid(material.reshape(cells,order='F'),size,origin,comments)
2019-11-23 01:22:36 +05:30
2021-06-16 00:53:10 +05:30
@staticmethod
2022-01-22 12:35:49 +05:30
def load_Neper(fname: Union[str, Path]) -> 'Grid':
2021-06-16 00:53:10 +05:30
"""
Load from Neper VTK file.
Parameters
----------
2022-01-13 01:04:29 +05:30
fname : str or pathlib.Path
2021-06-16 00:53:10 +05:30
Geometry file to read.
Returns
-------
loaded : damask.Grid
Grid-based geometry from file.
"""
2022-01-22 12:35:49 +05:30
v = VTK.load(fname,'ImageData')
2021-06-16 00:53:10 +05:30
cells = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
return Grid(v.get('MaterialId').reshape(cells,order='F').astype('int32',casting='unsafe') - 1,
bbox[1] - bbox[0], bbox[0],
2021-06-16 00:53:10 +05:30
util.execution_stamp('Grid','load_Neper'))
2020-10-09 17:49:19 +05:30
@staticmethod
2022-01-13 01:04:29 +05:30
def load_DREAM3D(fname: Union[str, Path],
2021-12-06 18:52:52 +05:30
feature_IDs: str = None, cell_data: str = None,
phases: str = 'Phases', Euler_angles: str = 'EulerAngles',
2022-01-22 12:35:49 +05:30
base_group: str = None) -> 'Grid':
2020-10-09 17:49:19 +05:30
"""
2021-03-20 18:07:06 +05:30
Load DREAM.3D (HDF5) file.
Data in DREAM.3D files can be stored per cell ('CellData') and/or
per grain ('Grain Data'). Per default, cell-wise data is assumed.
2021-03-20 18:07:06 +05:30
damask.ConfigMaterial.load_DREAM3D gives the corresponding material definition.
2020-10-09 17:49:19 +05:30
Parameters
----------
2022-01-13 01:04:29 +05:30
fname : str or or pathlib.Path
Filename of the DREAM.3D (HDF5) file.
2022-01-13 01:04:29 +05:30
feature_IDs : str, optional
2021-03-20 18:07:06 +05:30
Name of the dataset containing the mapping between cells and
grain-wise data. Defaults to 'None', in which case cell-wise
data is used.
2022-01-13 01:04:29 +05:30
cell_data : str, optional
2021-03-23 18:58:56 +05:30
Name of the group (folder) containing cell-wise data. Defaults to
None in wich case it is automatically detected.
2022-01-13 01:04:29 +05:30
phases : str, optional
2021-03-20 18:07:06 +05:30
Name of the dataset containing the phase ID. It is not used for
grain-wise data, i.e. when feature_IDs is not None.
Defaults to 'Phases'.
2022-01-13 01:04:29 +05:30
Euler_angles : str, optional
2021-03-20 18:07:06 +05:30
Name of the dataset containing the crystallographic orientation as
Euler angles in radians It is not used for grain-wise data, i.e.
when feature_IDs is not None. Defaults to 'EulerAngles'.
2022-01-13 01:04:29 +05:30
base_group : str, optional
2021-03-20 18:07:06 +05:30
Path to the group (folder) that contains geometry (_SIMPL_GEOMETRY),
and grain- or cell-wise data. Defaults to None, in which case
2021-03-20 04:19:41 +05:30
it is set as the path that contains _SIMPL_GEOMETRY/SPACING.
Returns
-------
loaded : damask.Grid
2021-04-24 18:17:52 +05:30
Grid-based geometry from file.
2021-03-20 18:07:06 +05:30
2020-10-09 17:49:19 +05:30
"""
b = util.DREAM3D_base_group(fname) if base_group is None else base_group
c = util.DREAM3D_cell_data_group(fname) if cell_data is None else cell_data
2020-10-09 17:49:19 +05:30
f = h5py.File(fname, 'r')
2021-03-23 18:58:56 +05:30
cells = f['/'.join([b,'_SIMPL_GEOMETRY','DIMENSIONS'])][()]
size = f['/'.join([b,'_SIMPL_GEOMETRY','SPACING'])] * cells
origin = f['/'.join([b,'_SIMPL_GEOMETRY','ORIGIN'])][()]
2021-03-20 04:19:41 +05:30
if feature_IDs is None:
phase = f['/'.join([b,c,phases])][()].reshape(-1,1)
O = Rotation.from_Euler_angles(f['/'.join([b,c,Euler_angles])]).as_quaternion().reshape(-1,4) # noqa
unique,unique_inverse = np.unique(np.hstack([O,phase]),return_inverse=True,axis=0)
ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \
np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse]
else:
ma = f['/'.join([b,c,feature_IDs])][()].flatten()
2020-10-09 17:49:19 +05:30
2020-12-04 11:42:18 +05:30
return Grid(ma.reshape(cells,order='F'),size,origin,util.execution_stamp('Grid','load_DREAM3D'))
2020-10-09 17:49:19 +05:30
@staticmethod
2022-01-13 01:04:29 +05:30
def from_table(table: Table,
coordinates: str,
2022-01-22 12:35:49 +05:30
labels: Union[str, Sequence[str]]) -> 'Grid':
2020-10-09 17:49:19 +05:30
"""
Create grid from ASCII table.
2020-10-09 17:49:19 +05:30
Parameters
----------
table : damask.Table
Table that contains material information.
2020-10-09 17:49:19 +05:30
coordinates : str
Label of the vector column containing the spatial coordinates.
Need to be ordered (1./x fast, 3./z slow).
2022-01-13 01:04:29 +05:30
labels : (list of) str
2020-10-09 17:49:19 +05:30
Label(s) of the columns containing the material definition.
Each unique combination of values results in one material ID.
2020-10-09 17:49:19 +05:30
Returns
-------
new : damask.Grid
2021-04-24 18:17:52 +05:30
Grid-based geometry from values in table.
2020-10-09 17:49:19 +05:30
"""
cells,size,origin = grid_filters.cellsSizeOrigin_coordinates0_point(table.get(coordinates))
2020-10-09 17:49:19 +05:30
labels_ = [labels] if isinstance(labels,str) else labels
unique,unique_inverse = np.unique(np.hstack([table.get(l) for l in labels_]),return_inverse=True,axis=0)
2020-10-31 21:53:58 +05:30
2020-12-04 02:28:24 +05:30
ma = np.arange(cells.prod()) if len(unique) == cells.prod() else \
2020-10-31 21:53:58 +05:30
np.arange(unique.size)[np.argsort(pd.unique(unique_inverse))][unique_inverse]
2020-10-09 17:49:19 +05:30
2020-12-04 11:42:18 +05:30
return Grid(ma.reshape(cells,order='F'),size,origin,util.execution_stamp('Grid','from_table'))
2020-10-09 17:49:19 +05:30
2020-03-29 22:42:23 +05:30
@staticmethod
def _find_closest_seed(seeds: np.ndarray,
weights: np.ndarray,
point: np.ndarray) -> np.integer:
2020-03-29 22:42:23 +05:30
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
2020-06-24 21:35:12 +05:30
2020-03-29 22:42:23 +05:30
@staticmethod
2022-01-13 01:04:29 +05:30
def from_Laguerre_tessellation(cells: IntSequence,
size: FloatSequence,
seeds: np.ndarray,
weights: FloatSequence,
material: IntSequence = None,
periodic: bool = True):
2020-03-29 22:42:23 +05:30
"""
Create grid from Laguerre tessellation.
2020-03-29 22:42:23 +05:30
Parameters
----------
2022-01-13 01:04:29 +05:30
cells : sequence of int, len (3)
2020-12-04 02:28:24 +05:30
Number of cells in x,y,z direction.
2022-01-13 01:04:29 +05:30
size : sequence of float, len (3)
2020-12-04 11:42:18 +05:30
Physical size of the grid in meter.
2022-01-13 03:43:38 +05:30
seeds : numpy.ndarray, shape (:,3)
2020-08-08 23:12:34 +05:30
Position of the seed points in meter. All points need to lay within the box.
2022-01-13 01:04:29 +05:30
weights : sequence of float, len (seeds.shape[0])
2020-08-08 23:12:34 +05:30
Weights of the seeds. Setting all weights to 1.0 gives a standard Voronoi tessellation.
2022-01-13 01:04:29 +05:30
material : sequence of int, len (seeds.shape[0]), optional
Material ID of the seeds.
Defaults to None, in which case materials are consecutively numbered.
2022-01-13 03:43:38 +05:30
periodic : bool, optional
2020-12-04 11:42:18 +05:30
Assume grid to be periodic. Defaults to True.
2020-03-29 22:42:23 +05:30
Returns
-------
new : damask.Grid
2021-04-24 18:17:52 +05:30
Grid-based geometry from tessellation.
2020-03-29 22:42:23 +05:30
"""
weights_p: FloatSequence
2020-03-29 22:42:23 +05:30
if periodic:
2020-04-20 23:46:25 +05:30
weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
2020-03-29 22:42:23 +05:30
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
else:
2022-01-20 17:22:56 +05:30
weights_p = np.array(weights,float)
2020-03-29 22:42:23 +05:30
seeds_p = seeds
coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
2022-01-13 01:04:29 +05:30
pool = mp.Pool(int(os.environ.get('OMP_NUM_THREADS',4)))
2021-06-29 14:14:13 +05:30
result = pool.map_async(partial(Grid._find_closest_seed,seeds_p,weights_p), coords)
2020-03-29 22:42:23 +05:30
pool.close()
pool.join()
material_ = np.array(result.get()).reshape(cells)
2020-03-29 22:42:23 +05:30
if periodic: material_ %= len(weights)
2020-03-29 22:42:23 +05:30
2022-01-13 01:04:29 +05:30
return Grid(material = material_ if material is None else np.array(material)[material_],
size = size,
2020-12-04 11:42:18 +05:30
comments = util.execution_stamp('Grid','from_Laguerre_tessellation'),
)
2020-03-29 22:42:23 +05:30
@staticmethod
2022-01-13 01:04:29 +05:30
def from_Voronoi_tessellation(cells: IntSequence,
size: FloatSequence,
2021-12-06 18:52:52 +05:30
seeds: np.ndarray,
2022-01-13 01:04:29 +05:30
material: IntSequence = None,
2022-01-22 12:35:49 +05:30
periodic: bool = True) -> 'Grid':
2020-03-29 22:42:23 +05:30
"""
Create grid from Voronoi tessellation.
2020-03-29 22:42:23 +05:30
Parameters
----------
2022-01-13 01:04:29 +05:30
cells : sequence of int, len (3)
2020-12-04 02:28:24 +05:30
Number of cells in x,y,z direction.
2022-01-13 01:04:29 +05:30
size : sequence of float, len (3)
2020-12-04 11:42:18 +05:30
Physical size of the grid in meter.
2022-01-13 03:43:38 +05:30
seeds : numpy.ndarray, shape (:,3)
2020-08-08 23:12:34 +05:30
Position of the seed points in meter. All points need to lay within the box.
2022-01-13 01:04:29 +05:30
material : sequence of int, len (seeds.shape[0]), optional
Material ID of the seeds.
Defaults to None, in which case materials are consecutively numbered.
2022-01-13 03:43:38 +05:30
periodic : bool, optional
2020-12-04 11:42:18 +05:30
Assume grid to be periodic. Defaults to True.
2020-03-29 22:42:23 +05:30
Returns
-------
new : damask.Grid
2021-04-24 18:17:52 +05:30
Grid-based geometry from tessellation.
2020-03-29 22:42:23 +05:30
"""
coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
2021-04-29 12:26:40 +05:30
tree = spatial.cKDTree(seeds,boxsize=size) if periodic else \
spatial.cKDTree(seeds)
try:
material_ = tree.query(coords, workers = int(os.environ.get('OMP_NUM_THREADS',4)))[1]
except TypeError:
material_ = tree.query(coords, n_jobs = int(os.environ.get('OMP_NUM_THREADS',4)))[1] # scipy <1.6
2020-03-29 22:42:23 +05:30
2022-01-13 01:04:29 +05:30
return Grid(material = (material_ if material is None else np.array(material)[material_]).reshape(cells),
size = size,
2020-12-04 11:42:18 +05:30
comments = util.execution_stamp('Grid','from_Voronoi_tessellation'),
)
2020-03-29 22:42:23 +05:30
_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) ),
}
2020-09-30 10:40:15 +05:30
@staticmethod
2022-01-13 01:04:29 +05:30
def from_minimal_surface(cells: IntSequence,
size: FloatSequence,
2021-12-06 18:52:52 +05:30
surface: str,
threshold: float = 0.0,
periods: int = 1,
2022-01-22 12:35:49 +05:30
materials: IntSequence = (0,1)) -> 'Grid':
2020-09-30 10:40:15 +05:30
"""
Create grid from definition of triply periodic minimal surface.
2020-09-30 10:40:15 +05:30
Parameters
----------
2022-01-13 01:04:29 +05:30
cells : sequence of int, len (3)
2020-12-04 02:28:24 +05:30
Number of cells in x,y,z direction.
2022-01-13 01:04:29 +05:30
size : sequence of float, len (3)
2020-12-04 11:42:18 +05:30
Physical size of the grid in meter.
surface : str
Type of the minimal surface. See notes for details.
2020-09-30 10:40:15 +05:30
threshold : float, optional.
Threshold of the minimal surface. Defaults to 0.0.
periods : integer, optional.
Number of periods per unit cell. Defaults to 1.
2022-01-13 01:04:29 +05:30
materials : sequence of int, len (2)
Material IDs. Defaults to (0,1).
2020-10-28 14:01:55 +05:30
Returns
-------
new : damask.Grid
2021-04-24 18:17:52 +05:30
Grid-based geometry from definition of minimal surface.
Notes
-----
The following triply-periodic minimal surfaces are implemented:
- Schwarz P
- Double Primitive
- Schwarz D
- Complementary D
- Double Diamond
- Dprime
- Gyroid
- Gprime
- Karcher K
- Lidinoid
- Neovius
- Fisher-Koch S
2020-09-30 10:40:15 +05:30
References
----------
2021-03-18 20:06:40 +05:30
S.B.G. Blanquer et al., Biofabrication 9(2):025001, 2017
https://doi.org/10.1088/1758-5090/aa6553
2021-03-18 20:06:40 +05:30
M. Wohlgemuth et al., Macromolecules 34(17):6083-6089, 2001
https://doi.org/10.1021/ma0019499
2021-03-18 20:06:40 +05:30
M.-T. Hsieh and L. Valdevit, Software Impacts 6:100026, 2020
https://doi.org/10.1016/j.simpa.2020.100026
Examples
--------
Minimal surface of 'Gyroid' type.
>>> import numpy as np
>>> import damask
2022-01-13 01:04:29 +05:30
>>> damask.Grid.from_minimal_surface([64]*3,np.ones(3)*1.e-4,'Gyroid')
cells : 64 x 64 x 64
size : 0.0001 x 0.0001 x 0.0001 /
origin: 0.0 0.0 0.0 / m
# materials: 2
Minimal surface of 'Neovius' type. non-default material IDs.
>>> import numpy as np
>>> import damask
2022-01-13 01:04:29 +05:30
>>> damask.Grid.from_minimal_surface([80]*3,np.ones(3)*5.e-4,
... 'Neovius',materials=(1,5))
2022-01-13 01:04:29 +05:30
cells : 80 x 80 x 80
size : 0.0005 x 0.0005 x 0.0005 /
origin: 0.0 0.0 0.0 / m
# materials: 2 (min: 1, max: 5)
2020-09-30 10:40:15 +05:30
"""
2020-12-04 02:28:24 +05:30
x,y,z = np.meshgrid(periods*2.0*np.pi*(np.arange(cells[0])+0.5)/cells[0],
periods*2.0*np.pi*(np.arange(cells[1])+0.5)/cells[1],
periods*2.0*np.pi*(np.arange(cells[2])+0.5)/cells[2],
2020-09-30 10:40:15 +05:30
indexing='ij',sparse=True)
2020-12-04 11:42:18 +05:30
return Grid(material = np.where(threshold < Grid._minimal_surface[surface](x,y,z),materials[1],materials[0]),
2020-09-30 10:40:15 +05:30
size = size,
2020-12-04 11:42:18 +05:30
comments = util.execution_stamp('Grid','from_minimal_surface'),
2020-09-30 10:40:15 +05:30
)
def save(self,
fname: Union[str, Path],
compress: bool = True):
"""
Save as VTK image data file.
Parameters
----------
2020-11-20 00:56:15 +05:30
fname : str or pathlib.Path
Filename to write. Valid extension is .vti, it will be appended if not given.
compress : bool, optional
Compress with zlib algorithm. Defaults to True.
"""
v = VTK.from_image_data(self.cells,self.size,self.origin)
v.add(self.material.flatten(order='F'),'material')
v.add_comments(self.comments)
2022-01-13 01:04:29 +05:30
v.save(fname,parallel=False,compress=compress)
def save_ASCII(self,
fname: Union[str, TextIO]):
2019-11-23 01:22:36 +05:30
"""
2020-12-04 02:28:24 +05:30
Save as geom file.
2019-11-23 01:22:36 +05:30
2020-11-14 22:24:47 +05:30
Storing geometry files in ASCII format is deprecated.
This function will be removed in a future version of DAMASK.
2019-11-23 01:22:36 +05:30
Parameters
----------
fname : str or file handle
Geometry file to write with extension '.geom'.
2020-09-18 20:02:08 +05:30
compress : bool, optional
Compress geometry with 'x of y' and 'a to b'.
2019-11-23 01:22:36 +05:30
"""
2022-01-05 18:36:06 +05:30
warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.0.0', DeprecationWarning,2)
header = [f'{len(self.comments)+4} header'] + self.comments \
2020-12-04 02:28:24 +05:30
+ ['grid a {} b {} c {}'.format(*self.cells),
'size x {} y {} z {}'.format(*self.size),
'origin x {} y {} z {}'.format(*self.origin),
'homogenization 1',
]
format_string = '%g' if self.material.dtype in np.sctypes['float'] else \
'%{}i'.format(1+int(np.floor(np.log10(np.nanmax(self.material)))))
np.savetxt(fname,
2020-12-04 02:28:24 +05:30
self.material.reshape([self.cells[0],np.prod(self.cells[1:])],order='F').T,
header='\n'.join(header), fmt=format_string, comments='')
2022-01-13 01:04:29 +05:30
def show(self) -> None:
2020-09-18 20:02:08 +05:30
"""Show on screen."""
VTK.from_image_data(self.cells,self.size,self.origin).show()
2021-12-06 18:52:52 +05:30
def add_primitive(self,
2022-01-13 01:04:29 +05:30
dimension: Union[FloatSequence, IntSequence],
center: Union[FloatSequence, IntSequence],
exponent: Union[FloatSequence, float],
2021-12-06 18:52:52 +05:30
fill: int = None,
R: Rotation = Rotation(),
inverse: bool = False,
2022-01-22 12:35:49 +05:30
periodic: bool = True) -> 'Grid':
2020-08-08 23:44:30 +05:30
"""
Insert a primitive geometric object at a given position.
2020-08-08 23:44:30 +05:30
Parameters
----------
2022-01-13 01:04:29 +05:30
dimension : sequence of int or float, len (3)
Dimension (diameter/side length) of the primitive.
If given as integers, cell centers are addressed.
If given as floats, physical coordinates are addressed.
center : sequence of int or float, len (3)
Center of the primitive.
If given as integers, cell centers are addressed.
If given as floats, physical coordinates are addressed.
exponent : float or sequence of float, len (3)
2020-09-26 21:32:25 +05:30
Exponents for the three axes.
2020-09-23 19:51:20 +05:30
0 gives octahedron (ǀxǀ^(2^0) + ǀyǀ^(2^0) + ǀzǀ^(2^0) < 1)
2020-09-26 21:32:25 +05:30
1 gives sphere (ǀxǀ^(2^1) + ǀyǀ^(2^1) + ǀzǀ^(2^1) < 1)
2020-08-09 00:26:17 +05:30
fill : int, optional
Fill value for primitive. Defaults to material.max()+1.
2020-08-08 23:44:30 +05:30
R : damask.Rotation, optional
Rotation of primitive. Defaults to no rotation.
2022-01-13 03:43:38 +05:30
inverse : bool, optional
Retain original materials within primitive and fill outside.
Defaults to False.
2022-01-13 03:43:38 +05:30
periodic : bool, optional
2020-12-04 11:42:18 +05:30
Assume grid to be periodic. Defaults to True.
2020-08-08 23:44:30 +05:30
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
2021-07-02 09:46:12 +05:30
Examples
--------
Add a sphere at the center.
>>> import numpy as np
>>> import damask
>>> g = damask.Grid(np.zeros([64]*3,int), np.ones(3)*1e-4)
>>> g.add_primitive(np.ones(3)*5e-5,np.ones(3)*5e-5,1)
2022-01-13 01:04:29 +05:30
cells : 64 x 64 x 64
size : 0.0001 x 0.0001 x 0.0001 /
origin: 0.0 0.0 0.0 / m
2021-07-02 09:46:12 +05:30
# materials: 2
Add a cube at the origin.
>>> import numpy as np
>>> import damask
>>> g = damask.Grid(np.zeros([64]*3,int), np.ones(3)*1e-4)
>>> g.add_primitive(np.ones(3,int)*32,np.zeros(3),np.inf)
2022-01-13 01:04:29 +05:30
cells : 64 x 64 x 64
size : 0.0001 x 0.0001 x 0.0001 /
origin: 0.0 0.0 0.0 / m
2021-07-02 09:46:12 +05:30
# materials: 2
2020-08-08 23:44:30 +05:30
"""
# radius and center
2020-12-04 02:28:24 +05:30
r = np.array(dimension)/2.0*self.size/self.cells if np.array(dimension).dtype in np.sctypes['int'] else \
np.array(dimension)/2.0
2020-12-04 02:28:24 +05:30
c = (np.array(center) + .5)*self.size/self.cells if np.array(center).dtype in np.sctypes['int'] else \
(np.array(center) - self.origin)
2020-08-08 23:44:30 +05:30
coords = grid_filters.coordinates0_point(self.cells,self.size,
2020-12-04 02:28:24 +05:30
-(0.5*(self.size + (self.size/self.cells
if np.array(center).dtype in np.sctypes['int'] else
0)) if periodic else c))
2020-12-04 02:28:24 +05:30
coords_rot = R.broadcast_to(tuple(self.cells))@coords
2020-08-08 23:44:30 +05:30
2020-08-27 13:02:49 +05:30
with np.errstate(all='ignore'):
mask = np.sum(np.power(coords_rot/r,2.0**np.array(exponent)),axis=-1) > 1.0
2020-08-08 23:44:30 +05:30
if periodic: # translate back to center
2020-12-04 02:28:24 +05:30
mask = np.roll(mask,((c/self.size-0.5)*self.cells).round().astype(int),(0,1,2))
2020-08-08 23:44:30 +05:30
2020-12-04 11:42:18 +05:30
return Grid(material = np.where(np.logical_not(mask) if inverse else mask,
self.material,
np.nanmax(self.material)+1 if fill is None else fill),
size = self.size,
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','add_primitive')],
)
2020-08-08 23:44:30 +05:30
def mirror(self,
directions: Sequence[str],
reflect: bool = False) -> 'Grid':
"""
2020-12-04 11:42:18 +05:30
Mirror grid along given directions.
2019-11-24 13:22:46 +05:30
Parameters
----------
2022-01-30 03:46:57 +05:30
directions : (sequence of) {'x', 'y', 'z'}
2020-12-04 11:42:18 +05:30
Direction(s) along which the grid is mirrored.
reflect : bool, optional
2020-08-25 12:04:04 +05:30
Reflect (include) outermost layers. Defaults to False.
2019-11-24 13:22:46 +05:30
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
2021-07-02 09:46:12 +05:30
Examples
--------
Mirror along x- and y-direction.
>>> import numpy as np
>>> import damask
>>> g = damask.Grid(np.zeros([32]*3,int), np.ones(3)*1e-4)
>>> g.mirror('xy',True)
2022-01-13 01:04:29 +05:30
cells : 64 x 64 x 32
size : 0.0002 x 0.0002 x 0.0001 /
origin: 0.0 0.0 0.0 / m
2021-07-02 09:46:12 +05:30
# materials: 1
"""
if not set(directions).issubset(valid := ['x', 'y', 'z']):
2020-12-11 01:23:11 +05:30
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
2021-12-06 18:52:52 +05:30
limits: Sequence[Optional[int]] = [None,None] if reflect else [-2,0]
mat = self.material.copy()
if 'x' in directions:
mat = np.concatenate([mat,mat[limits[0]:limits[1]:-1,:,:]],0)
if 'y' in directions:
mat = np.concatenate([mat,mat[:,limits[0]:limits[1]:-1,:]],1)
if 'z' in directions:
mat = np.concatenate([mat,mat[:,:,limits[0]:limits[1]:-1]],2)
2020-03-21 15:37:21 +05:30
2020-12-04 11:42:18 +05:30
return Grid(material = mat,
2020-12-04 02:28:24 +05:30
size = self.size/self.cells*np.asarray(mat.shape),
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','mirror')],
)
def flip(self,
2022-01-30 03:46:57 +05:30
directions: Sequence[str]) -> 'Grid':
"""
2020-12-04 11:42:18 +05:30
Flip grid along given directions.
Parameters
----------
2022-01-30 03:46:57 +05:30
directions : (sequence of) {'x', 'y', 'z'}
2020-12-04 11:42:18 +05:30
Direction(s) along which the grid is flipped.
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
"""
if not set(directions).issubset(valid := ['x', 'y', 'z']):
2020-12-11 01:23:11 +05:30
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
mat = np.flip(self.material, [valid.index(d) for d in directions if d in valid])
2020-12-04 11:42:18 +05:30
return Grid(material = mat,
size = self.size,
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','flip')],
)
def scale(self,
cells: IntSequence,
periodic: bool = True) -> 'Grid':
2019-11-24 19:43:26 +05:30
"""
2020-12-04 11:42:18 +05:30
Scale grid to new cells.
2019-11-24 19:43:26 +05:30
Parameters
----------
2022-01-13 01:04:29 +05:30
cells : sequence of int, len (3)
2020-12-04 02:28:24 +05:30
Number of cells in x,y,z direction.
2022-01-13 03:43:38 +05:30
periodic : bool, optional
2020-12-04 11:42:18 +05:30
Assume grid to be periodic. Defaults to True.
2019-11-24 19:43:26 +05:30
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
2021-07-02 09:46:12 +05:30
Examples
--------
Double resolution.
>>> import numpy as np
>>> import damask
2021-07-02 11:18:01 +05:30
>>> g = damask.Grid(np.zeros([32]*3,int),np.ones(3)*1e-4)
>>> g.scale(g.cells*2)
2022-01-13 01:04:29 +05:30
cells : 64 x 64 x 64
size : 0.0001 x 0.0001 x 0.0001 /
origin: 0.0 0.0 0.0 / m
2021-07-02 09:46:12 +05:30
# materials: 1
2019-11-24 19:43:26 +05:30
"""
2020-12-04 11:42:18 +05:30
return Grid(material = ndimage.interpolation.zoom(
self.material,
2020-12-04 02:28:24 +05:30
cells/self.cells,
output=self.material.dtype,
order=0,
mode=('wrap' if periodic else 'nearest'),
prefilter=False
),
size = self.size,
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','scale')],
)
2019-11-24 19:43:26 +05:30
2022-01-13 01:04:29 +05:30
def clean(self,
stencil: int = 3,
selection: IntSequence = None,
2022-01-22 12:35:49 +05:30
periodic: bool = True) -> 'Grid':
"""
2020-12-04 11:42:18 +05:30
Smooth grid by selecting most frequent material index within given stencil at each location.
2019-11-24 13:22:46 +05:30
Parameters
----------
stencil : int, optional
2020-08-08 23:12:34 +05:30
Size of smoothing stencil.
2022-01-13 01:04:29 +05:30
selection : sequence of int, optional
Field values that can be altered. Defaults to all.
2022-01-13 03:43:38 +05:30
periodic : bool, optional
2020-12-04 11:42:18 +05:30
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
"""
2022-01-13 01:04:29 +05:30
def mostFrequent(arr: np.ndarray, selection = None):
me = arr[arr.size//2]
if selection is None or me in selection:
unique, inverse = np.unique(arr, return_inverse=True)
return unique[np.argmax(np.bincount(inverse))]
else:
return me
2020-12-04 11:42:18 +05:30
return Grid(material = ndimage.filters.generic_filter(
self.material,
mostFrequent,
size=(stencil if selection is None else stencil//2*2+1,)*3,
mode=('wrap' if periodic else 'nearest'),
extra_keywords=dict(selection=selection),
).astype(self.material.dtype),
size = self.size,
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','clean')],
)
2019-11-24 22:51:05 +05:30
2022-01-22 12:35:49 +05:30
def renumber(self) -> 'Grid':
"""
Renumber sorted material indices as 0,...,N-1.
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
"""
2020-10-28 14:13:20 +05:30
_,renumbered = np.unique(self.material,return_inverse=True)
2020-12-04 11:42:18 +05:30
return Grid(material = renumbered.reshape(self.cells),
size = self.size,
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','renumber')],
)
2020-05-24 12:36:42 +05:30
def rotate(self,
R: Rotation,
fill: int = None) -> 'Grid':
2020-05-30 21:01:50 +05:30
"""
2020-12-04 11:42:18 +05:30
Rotate grid (pad if required).
2020-05-30 21:01:50 +05:30
Parameters
----------
R : damask.Rotation
2020-12-04 11:42:18 +05:30
Rotation to apply to the grid.
2022-01-13 01:04:29 +05:30
fill : int, optional
Material index to fill the corners. Defaults to material.max() + 1.
2020-05-30 21:01:50 +05:30
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
2020-05-30 21:01:50 +05:30
"""
2021-01-13 16:03:28 +05:30
material = self.material
2020-05-24 12:36:42 +05:30
# These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
2020-05-24 22:00:45 +05:30
# see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
2021-01-13 16:03:28 +05:30
for angle,axes in zip(R.as_Euler_angles(degrees=True)[::-1], [(0,1),(1,2),(0,1)]):
2022-01-20 17:22:56 +05:30
material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,
output=self.material.dtype,
cval=np.nanmax(self.material) + 1 if fill is None else fill)
2021-01-13 16:03:28 +05:30
# avoid scipy interpolation errors for rotations close to multiples of 90°
material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \
np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes)
2020-05-24 12:36:42 +05:30
2021-01-13 16:03:28 +05:30
origin = self.origin-(np.asarray(material.shape)-self.cells)*.5 * self.size/self.cells
2020-05-24 12:36:42 +05:30
2021-01-13 16:03:28 +05:30
return Grid(material = material,
size = self.size/self.cells*np.asarray(material.shape),
origin = origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','rotate')],
)
2020-05-24 12:36:42 +05:30
2022-01-13 01:04:29 +05:30
def canvas(self,
cells: IntSequence = None,
offset: IntSequence = None,
2022-01-22 12:35:49 +05:30
fill: int = None) -> 'Grid':
2020-05-30 21:01:50 +05:30
"""
2020-12-04 11:42:18 +05:30
Crop or enlarge/pad grid.
2020-05-30 21:01:50 +05:30
Parameters
----------
2022-01-13 01:04:29 +05:30
cells : sequence of int, len (3), optional
2020-12-04 02:28:24 +05:30
Number of cells x,y,z direction.
2022-01-13 01:04:29 +05:30
offset : sequence of int, len (3), optional
2020-12-04 11:42:18 +05:30
Offset (measured in cells) from old to new grid [0,0,0].
2022-01-13 01:04:29 +05:30
fill : int, optional
Material index to fill the background. Defaults to material.max() + 1.
2020-05-30 21:01:50 +05:30
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
2021-07-02 09:46:12 +05:30
Examples
--------
Remove 1/2 of the microstructure in z-direction.
>>> import numpy as np
>>> import damask
>>> g = damask.Grid(np.zeros([32]*3,int),np.ones(3)*1e-4)
2022-01-13 01:04:29 +05:30
>>> g.canvas([32,32,16])
cells : 33 x 32 x 16
size : 0.0001 x 0.0001 x 5e-05 /
origin: 0.0 0.0 0.0 / m
2021-07-02 09:46:12 +05:30
# materials: 1
2020-05-30 21:01:50 +05:30
"""
2022-01-13 01:04:29 +05:30
offset_ = np.array(offset,int) if offset is not None else np.zeros(3,int)
cells_ = np.array(cells,int) if cells is not None else self.cells
2020-05-24 12:36:42 +05:30
2022-01-20 17:22:56 +05:30
canvas = np.full(cells_,np.nanmax(self.material) + 1 if fill is None else fill,self.material.dtype)
2020-05-24 12:36:42 +05:30
2022-01-13 01:04:29 +05:30
LL = np.clip( offset_, 0,np.minimum(self.cells, cells_+offset_))
UR = np.clip( offset_+cells_, 0,np.minimum(self.cells, cells_+offset_))
ll = np.clip(-offset_, 0,np.minimum( cells_,self.cells-offset_))
ur = np.clip(-offset_+self.cells,0,np.minimum( cells_,self.cells-offset_))
2020-05-24 12:36:42 +05:30
canvas[ll[0]:ur[0],ll[1]:ur[1],ll[2]:ur[2]] = self.material[LL[0]:UR[0],LL[1]:UR[1],LL[2]:UR[2]]
2020-05-24 12:36:42 +05:30
2020-12-04 11:42:18 +05:30
return Grid(material = canvas,
2020-12-04 02:28:24 +05:30
size = self.size/self.cells*np.asarray(canvas.shape),
2022-01-13 01:04:29 +05:30
origin = self.origin+offset_*self.size/self.cells,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','canvas')],
)
2020-05-24 12:36:42 +05:30
def substitute(self,
from_material: IntSequence,
to_material: IntSequence) -> 'Grid':
2020-05-30 21:01:50 +05:30
"""
Substitute material indices.
2020-05-30 21:01:50 +05:30
Parameters
----------
2022-01-13 01:04:29 +05:30
from_material : sequence of int
Material indices to be substituted.
2022-01-13 01:04:29 +05:30
to_material : sequence of int
New material indices.
2020-05-30 21:01:50 +05:30
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
2020-05-30 21:01:50 +05:30
"""
2022-01-15 01:14:34 +05:30
material = self.material.copy()
for f,t in zip(from_material,to_material): # ToDo Python 3.10 has strict mode for zip
material[self.material==f] = t
2020-05-24 12:36:42 +05:30
2022-01-15 01:14:34 +05:30
return Grid(material = material,
size = self.size,
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','substitute')],
)
2020-08-08 23:12:34 +05:30
2022-01-22 12:35:49 +05:30
def sort(self) -> 'Grid':
"""
Sort material indices such that min(material) is located at (0,0,0).
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
"""
a = self.material.flatten(order='F')
from_ma = pd.unique(a)
sort_idx = np.argsort(from_ma)
ma = np.unique(a)[sort_idx][np.searchsorted(from_ma,a,sorter = sort_idx)]
2020-12-04 11:42:18 +05:30
return Grid(material = ma.reshape(self.cells,order='F'),
size = self.size,
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','sort')],
)
2021-12-06 18:52:52 +05:30
def vicinity_offset(self,
vicinity: int = 1,
offset: int = None,
2022-01-13 01:04:29 +05:30
trigger: IntSequence = [],
2022-01-22 12:35:49 +05:30
periodic: bool = True) -> 'Grid':
2020-08-08 23:12:34 +05:30
"""
Offset material index of points in the vicinity of xxx.
2020-08-08 23:12:34 +05:30
Different from themselves (or listed as triggers) within a given (cubic) vicinity,
i.e. within the region close to a grain/phase boundary.
2020-12-04 11:42:18 +05:30
ToDo: use include/exclude as in seeds.from_grid
2020-08-08 23:12:34 +05:30
Parameters
----------
vicinity : int, optional
Voxel distance checked for presence of other materials.
2020-08-08 23:12:34 +05:30
Defaults to 1.
offset : int, optional
Offset (positive or negative) to tag material indices,
defaults to material.max()+1.
2022-01-13 01:04:29 +05:30
trigger : sequence of int, optional
List of material indices that trigger a change.
Defaults to [], meaning that any different neighbor triggers a change.
2022-01-13 03:43:38 +05:30
periodic : bool, optional
2020-12-04 11:42:18 +05:30
Assume grid to be periodic. Defaults to True.
2020-08-08 23:12:34 +05:30
Returns
-------
updated : damask.Grid
2021-04-24 18:17:52 +05:30
Updated grid-based geometry.
2020-08-08 23:12:34 +05:30
"""
2022-01-13 01:04:29 +05:30
def tainted_neighborhood(stencil: np.ndarray, trigger):
2020-08-08 23:12:34 +05:30
me = stencil[stencil.shape[0]//2]
return np.any(stencil != me if len(trigger) == 0 else
2020-11-03 04:50:52 +05:30
np.in1d(stencil,np.array(list(set(trigger) - {me}))))
2020-08-08 23:12:34 +05:30
offset_ = np.nanmax(self.material)+1 if offset is None else offset
mask = ndimage.filters.generic_filter(self.material,
2020-08-08 23:12:34 +05:30
tainted_neighborhood,
size=1+2*vicinity,
mode='wrap' if periodic else 'nearest',
2020-08-08 23:12:34 +05:30
extra_keywords={'trigger':trigger})
2020-12-04 11:42:18 +05:30
return Grid(material = np.where(mask, self.material + offset_,self.material),
size = self.size,
origin = self.origin,
2020-12-04 11:42:18 +05:30
comments = self.comments+[util.execution_stamp('Grid','vicinity_offset')],
)
2020-11-28 00:46:06 +05:30
def get_grain_boundaries(self,
periodic: bool = True,
directions: Sequence[str] = 'xyz') -> VTK:
2020-11-18 16:55:08 +05:30
"""
2020-11-28 00:46:06 +05:30
Create VTK unstructured grid containing grain boundaries.
2020-11-18 16:44:12 +05:30
Parameters
----------
2022-01-13 03:43:38 +05:30
periodic : bool, optional
2020-12-04 11:42:18 +05:30
Assume grid to be periodic. Defaults to True.
2022-01-30 03:46:57 +05:30
directions : (sequence of) {'x', 'y', 'z'}, optional
2020-12-04 11:42:18 +05:30
Direction(s) along which the boundaries are determined.
2022-01-30 03:46:57 +05:30
Defaults to 'xyz'.
2020-11-18 16:44:12 +05:30
2021-04-25 11:17:00 +05:30
Returns
-------
grain_boundaries : damask.VTK
VTK-based geometry of grain boundary network.
2020-11-18 16:44:12 +05:30
"""
if not set(directions).issubset(valid := ['x', 'y', 'z']):
2020-12-11 01:23:11 +05:30
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
2020-11-19 00:40:04 +05:30
2020-12-04 02:28:24 +05:30
o = [[0, self.cells[0]+1, np.prod(self.cells[:2]+1)+self.cells[0]+1, np.prod(self.cells[:2]+1)],
[0, np.prod(self.cells[:2]+1), np.prod(self.cells[:2]+1)+1, 1],
[0, 1, self.cells[0]+1+1, self.cells[0]+1]] # offset for connectivity
2020-11-28 00:46:06 +05:30
connectivity = []
for i,d in enumerate(['x','y','z']):
if d not in directions: continue
2020-11-28 00:46:06 +05:30
mask = self.material != np.roll(self.material,1,i)
for j in [0,1,2]:
mask = np.concatenate((mask,np.take(mask,[0],j)*(i==j)),j)
if i == 0 and not periodic: mask[0,:,:] = mask[-1,:,:] = False
if i == 1 and not periodic: mask[:,0,:] = mask[:,-1,:] = False
if i == 2 and not periodic: mask[:,:,0] = mask[:,:,-1] = False
2020-11-28 00:46:06 +05:30
base_nodes = np.argwhere(mask.flatten(order='F')).reshape(-1,1)
connectivity.append(np.block([base_nodes + o[i][k] for k in range(4)]))
coords = grid_filters.coordinates0_node(self.cells,self.size,self.origin).reshape(-1,3,order='F')
return VTK.from_unstructured_grid(coords,np.vstack(connectivity),'QUAD')