Merge branch 'python-documentation' into drop-old-DADF5-support

This commit is contained in:
Martin Diehl 2021-04-24 19:39:28 +02:00
commit 6e1fe712c7
27 changed files with 958 additions and 439 deletions

8
.gitattributes vendored
View File

@ -12,8 +12,12 @@
*.pbz2 binary
# ignore files from MSC.Marc in language statistics
installation/mods_MarcMentat/20*/* linguist-vendored
installation/mods_MarcMentat/** linguist-vendored
src/Marc/include/* linguist-vendored
installation/mods_MarcMentat/apply_DAMASK_modifications.py linguist-vendored=false
# ignore reference files for tests in language statistics
python/tests/reference/* linguist-vendored
python/tests/reference/** linguist-vendored
# ignore deprecated scripts
processing/legacy/** linguist-vendored

View File

@ -233,7 +233,7 @@ source_distribution:
stage: deploy
script:
- cd $(mktemp -d)
- $DAMASKROOT/PRIVATE/releasing/deployMe.sh $CI_COMMIT_SHA
- $DAMASKROOT/PRIVATE/releasing/deploy.sh $DAMASKROOT $CI_COMMIT_SHA
except:
- master
- release

@ -1 +1 @@
Subproject commit 1490f97417664d6ae11d7ceafb6b98c9cc2dade1
Subproject commit afffa8d04e110282e514a4e57d0bad9c76effe01

View File

@ -1 +1 @@
v3.0.0-alpha2-866-g1be1a72a0
v3.0.0-alpha3-2-gd0dd1fd83

View File

@ -5,6 +5,7 @@ references:
10.1016/j.ijplas.2020.102779
- K. Sedighiani et al.,
Mechanics of Materials, submitted
output: [rho_dip, rho_mob]
N_sl: [12, 12]
b_sl: [2.49e-10, 2.49e-10]
rho_mob_0: [2.81e12, 2.8e12]

View File

@ -1,11 +1,4 @@
"""
Tools for pre and post processing of DAMASK simulations.
Modules that contain only one class (of the same name),
are prefixed by a '_'. For example, '_colormap' contains
a class called 'Colormap' which is imported as 'damask.Colormap'.
"""
"""Tools for managing DAMASK simulations."""
from pathlib import Path as _Path
import re as _re
@ -15,7 +8,6 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f:
version = _re.sub(r'^v','',_f.readline().strip())
__version__ = version
# make classes directly accessible as damask.Class
from . import util # noqa
from . import seeds # noqa
from . import tensor # noqa
@ -23,6 +15,8 @@ from . import mechanics # noqa
from . import solver # noqa
from . import grid_filters # noqa
from . import lattice # noqa
#Modules that contain only one class (of the same name), are prefixed by a '_'.
#For example, '_colormap' containsa class called 'Colormap' which is imported as 'damask.Colormap'.
from ._rotation import Rotation # noqa
from ._orientation import Orientation # noqa
from ._table import Table # noqa

View File

@ -91,6 +91,16 @@ class Colormap(mpl.colors.ListedColormap):
- 'lab': CIE Lab.
- 'msh': Msh (for perceptual uniform interpolation).
Returns
-------
new : damask.Colormap
Colormap within given bounds.
Examples
--------
>>> import damask
>>> damask.Colormap.from_range((0,0,1),(0,0,0),'blue_to_black')
"""
low_high = np.vstack((low,high))
if model.lower() == 'rgb':
@ -150,6 +160,16 @@ class Colormap(mpl.colors.ListedColormap):
This parameter is not used for matplotlib colormaps
that are of type `ListedColormap`.
Returns
-------
new : damask.Colormap
Predefined colormap.
Examples
--------
>>> import damask
>>> damask.Colormap.from_predefined('strain')
"""
# matplotlib presets
try:
@ -220,6 +240,11 @@ class Colormap(mpl.colors.ListedColormap):
damask.Colormap
The reversed colormap.
Examples
--------
>>> import damask
>>> damask.Colormap.from_predefined('stress').reversed()
"""
rev = super(Colormap,self).reversed(name)
return Colormap(np.array(rev.colors),rev.name[:-4] if rev.name.endswith('_r_r') else rev.name)
@ -239,8 +264,8 @@ class Colormap(mpl.colors.ListedColormap):
Returns
-------
f
File handle
f : file object
File handle with write access.
"""
if fname is None:

View File

@ -60,54 +60,6 @@ class ConfigMaterial(Config):
return super(ConfigMaterial,cls).load(fname)
@staticmethod
def from_table(table,**kwargs):
"""
Generate from an ASCII table.
Parameters
----------
table : damask.Table
Table that contains material information.
**kwargs
Keyword arguments where the key is the name and the value specifies
the label of the data column in the table.
Examples
--------
>>> import damask
>>> import damask.ConfigMaterial as cm
>>> t = damask.Table.load('small.txt')
>>> t
pos pos pos qu qu qu qu phase homog
0 0 0 0 0.19 0.8 0.24 -0.51 Aluminum SX
1 1 0 0 0.8 0.19 0.24 -0.51 Steel SX
1 1 1 0 0.8 0.19 0.24 -0.51 Steel SX
>>> cm.from_table(t,O='qu',phase='phase',homogenization='homog')
material:
- constituents:
- O: [0.19, 0.8, 0.24, -0.51]
v: 1.0
phase: Aluminum
homogenization: SX
- constituents:
- O: [0.8, 0.19, 0.24, -0.51]
v: 1.0
phase: Steel
homogenization: SX
homogenization: {}
phase: {}
"""
kwargs_ = {k:table.get(v) for k,v in kwargs.items()}
_,idx = np.unique(np.hstack(list(kwargs_.values())),return_index=True,axis=0)
idx = np.sort(idx)
kwargs_ = {k:np.atleast_1d(v[idx].squeeze()) for k,v in kwargs_.items()}
return ConfigMaterial().material_add(**kwargs_)
@staticmethod
def load_DREAM3D(fname,
grain_data=None,cell_data=None,cell_ensemble_data='CellEnsembleData',
@ -181,9 +133,69 @@ class ConfigMaterial(Config):
return base_config.material_add(**constituent,homogenization='direct')
@staticmethod
def from_table(table,**kwargs):
"""
Generate from an ASCII table.
Parameters
----------
table : damask.Table
Table that contains material information.
**kwargs
Keyword arguments where the key is the name and the value specifies
the label of the data column in the table.
Examples
--------
>>> import damask
>>> import damask.ConfigMaterial as cm
>>> t = damask.Table.load('small.txt')
>>> t
pos pos pos qu qu qu qu phase homog
0 0 0 0 0.19 0.8 0.24 -0.51 Aluminum SX
1 1 0 0 0.8 0.19 0.24 -0.51 Steel SX
1 1 1 0 0.8 0.19 0.24 -0.51 Steel SX
>>> cm.from_table(t,O='qu',phase='phase',homogenization='homog')
material:
- constituents:
- O: [0.19, 0.8, 0.24, -0.51]
v: 1.0
phase: Aluminum
homogenization: SX
- constituents:
- O: [0.8, 0.19, 0.24, -0.51]
v: 1.0
phase: Steel
homogenization: SX
homogenization: {}
phase: {}
"""
kwargs_ = {k:table.get(v) for k,v in kwargs.items()}
_,idx = np.unique(np.hstack(list(kwargs_.values())),return_index=True,axis=0)
idx = np.sort(idx)
kwargs_ = {k:np.atleast_1d(v[idx].squeeze()) for k,v in kwargs_.items()}
return ConfigMaterial().material_add(**kwargs_)
@property
def is_complete(self):
"""Check for completeness."""
"""
Check for completeness.
Only the general file layout is considered.
This check does not consider whether parameters for
a particular phase/homogenization model are missing.
Returns
-------
complete : bool
Whether the material.yaml definition is complete.
"""
ok = True
for top_level in ['homogenization','phase','material']:
ok &= top_level in self
@ -236,7 +248,19 @@ class ConfigMaterial(Config):
@property
def is_valid(self):
"""Check for valid content."""
"""
Check for valid content.
Only the generic file content is considered.
This check does not consider whether parameters for a
particular phase/homogenization mode are out of bounds.
Returns
-------
valid : bool
Whether the material.yaml definition is valid.
"""
ok = True
if 'phase' in self:
@ -282,7 +306,7 @@ class ConfigMaterial(Config):
Returns
-------
cfg : damask.ConfigMaterial
updated : damask.ConfigMaterial
Updated material configuration.
"""
@ -311,7 +335,7 @@ class ConfigMaterial(Config):
Returns
-------
cfg : damask.ConfigMaterial
updated : damask.ConfigMaterial
Updated material configuration.
"""
@ -336,14 +360,17 @@ class ConfigMaterial(Config):
Returns
-------
cfg : damask.ConfigMaterial
updated : damask.ConfigMaterial
Updated material configuration.
Examples
--------
Create a dual-phase steel microstructure for micromechanical simulations:
>>> import numpy as np
>>> import damask
>>> m = damask.ConfigMaterial().material_add(phase = ['Aluminum','Steel'],
>>> m = damask.ConfigMaterial()
>>> m = m.material_add(phase = ['Ferrite','Martensite'],
... O = damask.Rotation.from_random(2),
... homogenization = 'SX')
>>> m
@ -351,38 +378,43 @@ class ConfigMaterial(Config):
- constituents:
- O: [0.577764, -0.146299, -0.617669, 0.513010]
v: 1.0
phase: Aluminum
phase: Ferrite
homogenization: SX
- constituents:
- O: [0.184176, 0.340305, 0.737247, 0.553840]
v: 1.0
phase: Steel
phase: Martensite
homogenization: SX
homogenization: {}
phase: {}
>>> m = damask.ConfigMaterial().material_add(phase = np.array(['Austenite','Martensite']).reshape(1,2),
Create a duplex stainless steel microstructure for forming simulations:
>>> import numpy as np
>>> import damask
>>> m = damask.ConfigMaterial()
>>> m = m.material_add(phase = np.array(['Austenite','Ferrite']).reshape(1,2),
... O = damask.Rotation.from_random((2,2)),
... v = np.array([0.2,0.8]).reshape(1,2),
... homogenization = ['A','B'])
... homogenization = 'Taylor')
>>> m
material:
- constituents:
- phase: Austenite
O: [0.659802978293224, 0.6953785848195171, 0.22426295326327111, -0.17554139512785227]
v: 0.2
- phase: Martensite
- phase: Ferrite
O: [0.49356745891301596, 0.2841806579193434, -0.7487679215072818, -0.339085707289975]
v: 0.8
homogenization: A
homogenization: Taylor
- constituents:
- phase: Austenite
O: [0.26542221365204055, 0.7268854930702071, 0.4474726435701472, -0.44828201137283735]
v: 0.2
- phase: Martensite
- phase: Ferrite
O: [0.6545817158479885, -0.08004812803625233, -0.6226561293931374, 0.4212059104577611]
v: 0.8
homogenization: B
homogenization: Taylor
homogenization: {}
phase: {}

View File

@ -161,6 +161,11 @@ class Grid:
Grid file to read. Valid extension is .vtr, which will be appended
if not given.
Returns
-------
loaded : damask.Grid
Grid-based geometry from file.
"""
v = VTK.load(fname if str(fname).endswith('.vtr') else str(fname)+'.vtr')
comments = v.get_comments()
@ -190,6 +195,11 @@ class Grid:
fname : str, pathlib.Path, or file handle
Geometry file to read.
Returns
-------
loaded : damask.Grid
Grid-based geometry from file.
"""
warnings.warn('Support for ASCII-based geom format will be removed in DAMASK 3.1.0', DeprecationWarning,2)
try:
@ -281,6 +291,10 @@ class Grid:
and grain- or cell-wise data. Defaults to None, in which case
it is set as the path that contains _SIMPL_GEOMETRY/SPACING.
Returns
-------
loaded : damask.Grid
Grid-based geometry from file.
"""
b = util.DREAM3D_base_group(fname) if base_group is None else base_group
@ -319,6 +333,11 @@ class Grid:
Label(s) of the columns containing the material definition.
Each unique combination of values results in one material ID.
Returns
-------
new : damask.Grid
Grid-based geometry from values in table.
"""
cells,size,origin = grid_filters.cellsSizeOrigin_coordinates0_point(table.get(coordinates))
@ -356,6 +375,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
new : damask.Grid
Grid-based geometry from tessellation.
"""
if periodic:
weights_p = np.tile(weights,27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
@ -405,6 +429,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
new : damask.Grid
Grid-based geometry from tessellation.
"""
coords = grid_filters.coordinates0_point(cells,size).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
@ -478,6 +507,11 @@ class Grid:
materials : (int, int), optional
Material IDs. Defaults to (1,2).
Returns
-------
new : damask.Grid
Grid-based geometry from definition of minimal surface.
Notes
-----
The following triply-periodic minimal surfaces are implemented:
@ -598,6 +632,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
# radius and center
r = np.array(dimension)/2.0*self.size/self.cells if np.array(dimension).dtype in np.sctypes['int'] else \
@ -638,6 +677,11 @@ class Grid:
reflect : bool, optional
Reflect (include) outermost layers. Defaults to False.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
@ -670,6 +714,11 @@ class Grid:
Direction(s) along which the grid is flipped.
Valid entries are 'x', 'y', 'z'.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
valid = ['x','y','z']
if not set(directions).issubset(valid):
@ -695,6 +744,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
return Grid(material = ndimage.interpolation.zoom(
self.material,
@ -723,6 +777,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
def mostFrequent(arr,selection=None):
me = arr[arr.size//2]
@ -746,7 +805,15 @@ class Grid:
def renumber(self):
"""Renumber sorted material indices as 0,...,N-1."""
"""
Renumber sorted material indices as 0,...,N-1.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
_,renumbered = np.unique(self.material,return_inverse=True)
return Grid(material = renumbered.reshape(self.cells),
@ -767,6 +834,11 @@ class Grid:
fill : int or float, optional
Material index to fill the corners. Defaults to material.max() + 1.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
if fill is None: fill = np.nanmax(self.material) + 1
dtype = float if isinstance(fill,float) or self.material.dtype in np.sctypes['float'] else int
@ -802,6 +874,11 @@ class Grid:
fill : int or float, optional
Material index to fill the background. Defaults to material.max() + 1.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
if offset is None: offset = 0
if fill is None: fill = np.nanmax(self.material) + 1
@ -834,6 +911,11 @@ class Grid:
to_material : iterable of ints
New material indices.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
def mp(entry,mapper):
return mapper[entry] if entry in mapper else entry
@ -849,7 +931,15 @@ class Grid:
def sort(self):
"""Sort material indices such that min(material) is located at (0,0,0)."""
"""
Sort material indices such that min(material) is located at (0,0,0).
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
a = self.material.flatten(order='F')
from_ma = pd.unique(a)
sort_idx = np.argsort(from_ma)
@ -884,6 +974,11 @@ class Grid:
periodic : Boolean, optional
Assume grid to be periodic. Defaults to True.
Returns
-------
updated : damask.Grid
Updated grid-based geometry.
"""
def tainted_neighborhood(stencil,trigger):

View File

@ -8,15 +8,15 @@ from . import tensor
_parameter_doc = \
"""lattice : str
Either a crystal family out of [triclinic, monoclinic, orthorhombic, tetragonal, hexagonal, cubic]
or a Bravais lattice out of [aP, mP, mS, oP, oS, oI, oF, tP, tI, hP, cP, cI, cF].
When specifying a Bravais lattice, additional lattice parameters might be required:
Either a crystal family out of {triclinic, monoclinic, orthorhombic, tetragonal, hexagonal, cubic}
or a Bravais lattice out of {aP, mP, mS, oP, oS, oI, oF, tP, tI, hP, cP, cI, cF}.
When specifying a Bravais lattice, additional lattice parameters might be required.
a : float, optional
Length of lattice parameter "a".
Length of lattice parameter 'a'.
b : float, optional
Length of lattice parameter "b".
Length of lattice parameter 'b'.
c : float, optional
Length of lattice parameter "c".
Length of lattice parameter 'c'.
alpha : float, optional
Angle between b and c lattice basis.
beta : float, optional

View File

@ -57,11 +57,30 @@ def _empty_like(dataset,N_materialpoints,fill_float,fill_int):
class Result:
"""
Manipulate and read DADF5 files.
Add data to and export data from a DADF5 file.
A DADF5 (DAMASK HDF5) file contain DAMASK results.
Its group/folder structure reflects the layout in material.yaml.
This class provides a customable view on the DADF5 file.
Upon initialization, all attributes are visible.
Derived quantities are added to the file and existing data is
exported based on the current view.
Examples
--------
Open 'my_file.hdf5', which needs to contain deformation gradient 'F'
and first Piola-Kirchhoff stress 'P', add the Mises equivalent of the
Cauchy stress, and export it to VTK (file) and numpy.ndarray (memory).
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_Cauchy()
>>> r.add_equivalent_Mises('sigma')
>>> r.save_VTK()
>>> r_last = r.view('increments',-1)
>>> sigma_vM_last = r_last.get('sigma_vM')
DADF5 (DAMASK HDF5) files contain DAMASK results.
The group/folder structure reflects the input data
in material.yaml.
"""
def __init__(self,fname):
@ -95,10 +114,10 @@ class Result:
self.N_materialpoints, self.N_constituents = np.shape(f[f'cell_to/phase'])
self.homogenizations = [m.decode() for m in np.unique(f[f'cell_to/homogenization']['label'])]
self.homogenizations = sorted(self.homogenizations,key=util.natural_sort)
self.phases = [c.decode() for c in np.unique(f[f'cell_to/phase']['label'])]
self.phases = sorted(self.phases,key=util.natural_sort)
self.homogenization = f[f'cell_to/homogenization']['label'].astype('str')
self.homogenizations = sorted(np.unique(self.homogenization),key=util.natural_sort)
self.phase = f[f'cell_to/phase']['label'].astype('str')
self.phases = sorted(np.unique(self.phase),key=util.natural_sort)
self.fields = []
for c in self.phases:
@ -154,6 +173,11 @@ class Result:
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
Returns
-------
view : damask.Result
Modified or new view on the DADF5 file.
"""
# allow True/False and string arguments
if datasets is True:
@ -166,6 +190,7 @@ class Result:
if what == 'increments':
choice = [c if isinstance(c,str) and c.startswith('increment_') else
f'increment_{c}' for c in choice]
if datasets == -1: choice = [self.increments[-1]]
elif what == 'times':
what = 'increments'
if choice == ['*']:
@ -197,15 +222,31 @@ class Result:
return dup
def allow_modification(self):
"""Allow to overwrite existing data."""
def modification_enable(self):
"""
Allow to modify existing data.
Returns
-------
modified_view : damask.Result
View where data is not write-protected.
"""
print(util.warn('Warning: Modification of existing datasets allowed!'))
dup = self.copy()
dup._allow_modification = True
return dup
def disallow_modification(self):
"""Disallow to overwrite existing data (default case)."""
def modification_disable(self):
"""
Disallow to modify existing data (default case).
Returns
-------
modified_view : damask.Result
View where data is write-protected.
"""
dup = self.copy()
dup._allow_modification = False
return dup
@ -213,7 +254,7 @@ class Result:
def increments_in_range(self,start,end):
"""
Select all increments within a given range.
Get all increments within a given range.
Parameters
----------
@ -222,6 +263,11 @@ class Result:
end : int or str
End increment.
Returns
-------
increments : list of ints
Increment number of all increments within the given bounds.
"""
selected = []
for i,inc in enumerate([int(i[10:]) for i in self.increments]):
@ -233,7 +279,7 @@ class Result:
def times_in_range(self,start,end):
"""
Select all increments within a given time range.
Get all increments within a given time range.
Parameters
----------
@ -242,6 +288,11 @@ class Result:
end : float
Time of end increment.
Returns
-------
times : list of float
Simulation time of all increments within the given bounds.
"""
selected = []
for i,time in enumerate(self.times):
@ -262,6 +313,25 @@ class Result:
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
Returns
-------
view : damask.Result
View with where selected attributes are visible.
Examples
--------
Get a view that shows only results from the initial configuration:
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r_first = r.view('increment',0)
Get a view that shows all results of in simulation time [10,40]:
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r_t10to40 = r.view('times',r.times_in_range(10.0,40.0))
"""
return self._manage_view('set',what,datasets)
@ -278,6 +348,20 @@ class Result:
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
Returns
-------
modified_view : damask.Result
View with more visible attributes.
Examples
--------
Get a view that shows only results from first and last increment:
>>> import damask
>>> r_empty = damask.Result('my_file.hdf5').view('increments',False)
>>> r_first = r_empty.view_more('increments',0)
>>> r_first_and_last = r.first.view_more('increments',-1)
"""
return self._manage_view('add',what,datasets)
@ -294,37 +378,98 @@ class Result:
Name of datasets; supports '?' and '*' wildcards.
True is equivalent to '*', False is equivalent to [].
Returns
-------
modified_view : damask.Result
View with less visible attributes.
Examples
--------
Get a view that does not show the undeformed configuration:
>>> import damask
>>> r_all = damask.Result('my_file.hdf5')
>>> r_deformed = r_all.view_less('increments',0)
"""
return self._manage_view('del',what,datasets)
def rename(self,name_old,name_new):
def rename(self,name_src,name_dst):
"""
Rename dataset.
Rename/move datasets (within the same group/folder).
This operation is discouraged because the history of the
data becomes untracable and scientific integrity cannot be
ensured.
Parameters
----------
name_old : str
Name of the dataset to be renamed.
name_new : str
New name of the dataset.
name_src : str
Name of the datasets to be renamed.
name_dst : str
New name of the datasets.
Examples
--------
Rename datasets containing the deformation gradient from 'F' to 'def_grad':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r_unprotected = r.modification_enable()
>>> r_unprotected.rename('F','def_grad')
"""
if not self._allow_modification:
raise PermissionError('Rename operation not permitted')
raise PermissionError('Renaming datasets not permitted')
with h5py.File(self.fname,'a') as f:
for inc in self.visible['increments']:
for ty in ['phase','homogenization']:
for label in self.visible[ty+'s']:
for field in self.visible['fields']:
path_old = '/'.join([inc,ty,label,field,name_old])
path_new = '/'.join([inc,ty,label,field,name_new])
if path_old in f.keys():
f[path_new] = f[path_old]
f[path_new].attrs['renamed'] = f'original name: {name_old}' if h5py3 else \
f'original name: {name_old}'.encode()
del f[path_old]
for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()):
path_src = '/'.join([inc,ty,label,field,name_src])
path_dst = '/'.join([inc,ty,label,field,name_dst])
if path_src in f.keys():
f[path_dst] = f[path_src]
f[path_dst].attrs['renamed'] = f'original name: {name_src}' if h5py3 else \
f'original name: {name_src}'.encode()
del f[path_src]
def remove(self,name):
"""
Remove/delete datasets.
This operation is discouraged because the history of the
data becomes untracable and scientific integrity cannot be
ensured.
Parameters
----------
name : str
Name of the datasets to be deleted.
Examples
--------
Delete the deformation gradient 'F':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r_unprotected = r.modification_enable()
>>> r_unprotected.remove('F')
"""
if not self._allow_modification:
raise PermissionError('Removing datasets not permitted')
with h5py.File(self.fname,'a') as f:
for inc in self.visible['increments']:
for ty in ['phase','homogenization']:
for label in self.visible[ty+'s']:
for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()):
path = '/'.join([inc,ty,label,field,name])
if path in f.keys(): del f[path]
def list_data(self):
@ -337,7 +482,7 @@ class Result:
msg = ' '.join([msg,f'{ty}\n'])
for label in self.visible[ty+'s']:
msg = ' '.join([msg,f'{label}\n'])
for field in self.visible['fields']:
for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()):
msg = ' '.join([msg,f'{field}\n'])
for d in f['/'.join([inc,ty,label,field])].keys():
dataset = f['/'.join([inc,ty,label,field,d])]
@ -357,7 +502,7 @@ class Result:
@property
def coordinates0_point(self):
"""Return initial coordinates of the cell centers."""
"""Initial/undeformed cell center coordinates."""
if self.structured:
return grid_filters.coordinates0_point(self.cells,self.size,self.origin).reshape(-1,3,order='F')
else:
@ -366,7 +511,7 @@ class Result:
@property
def coordinates0_node(self):
"""Return initial coordinates of the cell centers."""
"""Initial/undeformed nodal coordinates."""
if self.structured:
return grid_filters.coordinates0_node(self.cells,self.size,self.origin).reshape(-1,3,order='F')
else:
@ -375,6 +520,7 @@ class Result:
@property
def geometry0(self):
"""Initial/undeformed geometry."""
if self.structured:
return VTK.from_rectilinear_grid(self.cells,self.size,self.origin)
else:
@ -403,7 +549,7 @@ class Result:
Parameters
----------
x : str
Label of scalar, vector, or tensor dataset to take absolute value of.
Name of scalar, vector, or tensor dataset to take absolute value of.
"""
self._add_generic_pointwise(self._add_absolute,{'x':x})
@ -424,24 +570,51 @@ class Result:
'creator': 'add_calculation'
}
}
def add_calculation(self,label,formula,unit='n/a',description=None):
def add_calculation(self,name,formula,unit='n/a',description=None):
"""
Add result of a general formula.
Parameters
----------
label : str
Label of resulting dataset.
name : str
Name of resulting dataset.
formula : str
Formula to calculate resulting dataset. Existing datasets are referenced by '#TheirLabel#'.
Formula to calculate resulting dataset. Existing datasets are referenced by '#TheirName#'.
unit : str, optional
Physical unit of the result.
description : str, optional
Human-readable description of the result.
Examples
--------
Add total dislocation density, i.e. the sum of mobile dislocation
density 'rho_mob' and dislocation dipole density 'rho_dip' over
all slip systems:
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_calculation('rho_mob_total','np.sum(#rho_mob#,axis=1)',
... '1/m²','total mobile dislocation density')
>>> r.add_calculation('rho_dip_total','np.sum(#rho_dip#,axis=1)',
... '1/m²','total dislocation dipole density')
>>> r.add_calculation('rho_total','#rho_dip_total#+#rho_mob_total',
... '1/m²','total dislocation density')
Add Mises equivalent of the Cauchy stress without storage of
intermediate results. Define a user function for better readability:
>>> import damask
>>> def equivalent_stress(F,P):
... sigma = damask.mechanics.stress_Cauchy(F=F,P=P)
... return damask.mechanics.equivalent_stress_Mises(sigma)
>>> r = damask.Result('my_file.hdf5')
>>> r.enable_user_function(equivalent_stress)
>>> r.add_calculation('sigma_vM','equivalent_stress(#F#,#P#)','Pa',
... 'Mises equivalent of the Cauchy stress')
"""
dataset_mapping = {d:d for d in set(re.findall(r'#(.*?)#',formula))} # datasets used in the formula
args = {'formula':formula,'label':label,'unit':unit,'description':description}
args = {'formula':formula,'label':name,'unit':unit,'description':description}
self._add_generic_pointwise(self._add_calculation,dataset_mapping,args)
@ -465,9 +638,9 @@ class Result:
Parameters
----------
P : str, optional
Label of the dataset containing the first Piola-Kirchhoff stress. Defaults to 'P'.
Name of the dataset containing the first Piola-Kirchhoff stress. Defaults to 'P'.
F : str, optional
Label of the dataset containing the deformation gradient. Defaults to 'F'.
Name of the dataset containing the deformation gradient. Defaults to 'F'.
"""
self._add_generic_pointwise(self._add_stress_Cauchy,{'P':P,'F':F})
@ -491,7 +664,15 @@ class Result:
Parameters
----------
T : str
Label of tensor dataset.
Name of tensor dataset.
Examples
--------
Add the determinant of plastic deformation gradient 'F_p':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_determinant('F_p')
"""
self._add_generic_pointwise(self._add_determinant,{'T':T})
@ -515,7 +696,7 @@ class Result:
Parameters
----------
T : str
Label of tensor dataset.
Name of tensor dataset.
"""
self._add_generic_pointwise(self._add_deviator,{'T':T})
@ -546,7 +727,7 @@ class Result:
Parameters
----------
T_sym : str
Label of symmetric tensor dataset.
Name of symmetric tensor dataset.
eigenvalue : str, optional
Eigenvalue. Select from 'max', 'mid', 'min'. Defaults to 'max'.
@ -579,7 +760,7 @@ class Result:
Parameters
----------
T_sym : str
Label of symmetric tensor dataset.
Name of symmetric tensor dataset.
eigenvalue : str, optional
Eigenvalue to which the eigenvector corresponds.
Select from 'max', 'mid', 'min'. Defaults to 'max'.
@ -619,9 +800,17 @@ class Result:
l : numpy.array of shape (3)
Lab frame direction for inverse pole figure.
q : str
Label of the dataset containing the crystallographic orientation as quaternions.
Name of the dataset containing the crystallographic orientation as quaternions.
Defaults to 'O'.
Examples
--------
Add the IPF color along [0,1,1] for orientation 'O':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_IPF_color(np.array([0,1,1]))
"""
self._add_generic_pointwise(self._add_IPF_color,{'q':q},{'l':l})
@ -644,7 +833,7 @@ class Result:
Parameters
----------
T_sym : str
Label of symmetric tensor dataset.
Name of symmetric tensor dataset.
"""
self._add_generic_pointwise(self._add_maximum_shear,{'T_sym':T_sym})
@ -678,11 +867,25 @@ class Result:
Parameters
----------
T_sym : str
Label of symmetric tensorial stress or strain dataset.
Name of symmetric tensorial stress or strain dataset.
kind : {'stress', 'strain', None}, optional
Kind of the von Mises equivalent. Defaults to None, in which case
it is selected based on the unit of the dataset ('1' -> strain, 'Pa' -> stress).
Examples
--------
Add the Mises equivalent of the Cauchy stress 'sigma':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_equivalent_Mises('sigma')
Add the Mises equivalent of the spatial logarithmic strain 'epsilon_V^0.0(F)':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_equivalent_Mises('epsilon_V^0.0(F)')
"""
self._add_generic_pointwise(self._add_equivalent_Mises,{'T_sym':T_sym},{'kind':kind})
@ -717,7 +920,7 @@ class Result:
Parameters
----------
x : str
Label of vector or tensor dataset.
Name of vector or tensor dataset.
ord : {non-zero int, inf, -inf, 'fro', 'nuc'}, optional
Order of the norm. inf means NumPys inf object. For details refer to numpy.linalg.norm.
@ -745,9 +948,9 @@ class Result:
Parameters
----------
P : str, optional
Label of first Piola-Kirchhoff stress dataset. Defaults to 'P'.
Name of first Piola-Kirchhoff stress dataset. Defaults to 'P'.
F : str, optional
Label of deformation gradient dataset. Defaults to 'F'.
Name of deformation gradient dataset. Defaults to 'F'.
"""
self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F})
@ -786,7 +989,7 @@ class Result:
# Parameters
# ----------
# q : str
# Label of the dataset containing the crystallographic orientation as quaternions.
# Name of the dataset containing the crystallographic orientation as quaternions.
# p : numpy.array of shape (3)
# Crystallographic direction or plane.
# polar : bool, optional
@ -813,8 +1016,16 @@ class Result:
Parameters
----------
F : str, optional
Label of deformation gradient dataset.
F : str
Name of deformation gradient dataset.
Examples
--------
Add the rotational part of deformation gradient 'F':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_rotation('F')
"""
self._add_generic_pointwise(self._add_rotation,{'F':F})
@ -838,7 +1049,15 @@ class Result:
Parameters
----------
T : str
Label of tensor dataset.
Name of tensor dataset.
Examples
--------
Add the hydrostatic part of the Cauchy stress 'sigma':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.add_spherical('sigma')
"""
self._add_generic_pointwise(self._add_spherical,{'T':T})
@ -864,13 +1083,28 @@ class Result:
Parameters
----------
F : str, optional
Label of deformation gradient dataset. Defaults to 'F'.
Name of deformation gradient dataset. Defaults to 'F'.
t : {'V', 'U'}, optional
Type of the polar decomposition, 'V' for left stretch tensor and 'U' for right stretch tensor.
Defaults to 'V'.
m : float, optional
Order of the strain calculation. Defaults to 0.0.
Examples
--------
Add the Biot strain based on the deformation gradient 'F':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.strain(t='U',m=0.5)
Add the plastic Euler-Almansi strain based on the
plastic deformation gradient 'F_p':
>>> import damask
>>> r = damask.Result('my_file.hdf5')
>>> r.strain('F_p','V',-1)
"""
self._add_generic_pointwise(self._add_strain,{'F':F},{'t':t,'m':m})
@ -894,7 +1128,7 @@ class Result:
Parameters
----------
F : str, optional
Label of deformation gradient dataset. Defaults to 'F'.
Name of deformation gradient dataset. Defaults to 'F'.
t : {'V', 'U'}, optional
Type of the polar decomposition, 'V' for left stretch tensor and 'U' for right stretch tensor.
Defaults to 'V'.
@ -947,7 +1181,7 @@ class Result:
for inc in self.visible['increments']:
for ty in ['phase','homogenization']:
for label in self.visible[ty+'s']:
for field in self.visible['fields']:
for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()):
group = '/'.join([inc,ty,label,field])
if set(datasets.values()).issubset(f[group].keys()): groups.append(group)
@ -1001,10 +1235,14 @@ class Result:
"""
Write XDMF file to directly visualize data in DADF5 file.
The XDMF format is only supported for structured grids
with single phase and single constituent.
For other cases use `save_VTK`.
Parameters
----------
output : (list of) str
Labels of the datasets to read.
Names of the datasets included in the XDMF file.
Defaults to '*', in which case all datasets are considered.
"""
@ -1082,7 +1320,7 @@ class Result:
for ty in ['phase','homogenization']:
for label in self.visible[ty+'s']:
for field in self.visible['fields']:
for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()):
for out in _match(output,f['/'.join([inc,ty,label,field])].keys()):
name = '/'.join([inc,ty,label,field,out])
shape = f[name].shape[1:]
@ -1131,10 +1369,16 @@ class Result:
"""
Export to VTK cell/point data.
One VTK file per visible increment is created.
For cell data, the VTK format is a rectilinear grid (.vtr) for
grid-based simulations and an unstructured grid (.vtu) for
mesh-baed simulations. For point data, the VTK format is poly
data (.vtp).
Parameters
----------
output : (list of) str, optional
Labels of the datasets to place.
Names of the datasets included in the VTK file.
Defaults to '*', in which case all datasets are exported.
mode : {'cell', 'point'}
Export in cell format or point format.
@ -1212,7 +1456,7 @@ class Result:
Parameters
----------
output : (list of) str
Labels of the datasets to read.
Names of the datasets to read.
Defaults to '*', in which case all datasets are read.
flatten : bool
Remove singular levels of the folder hierarchy.
@ -1264,7 +1508,7 @@ class Result:
Parameters
----------
output : (list of) str, optional
Labels of the datasets to place.
Names of the datasets to read.
Defaults to '*', in which case all datasets are placed.
flatten : bool
Remove singular levels of the folder hierarchy.

View File

@ -24,7 +24,7 @@ class Rotation:
- Euler angle triplets are implemented using the Bunge convention,
with angular ranges of [0,2π], [0,π], [0,2π].
- The rotation angle ω is limited to the interval [0,π].
- The real part of a quaternion is positive, Re(q) > 0
- The real part of a quaternion is positive, Re(q) 0
- P = -1 (as default).
Examples
@ -357,7 +357,7 @@ class Rotation:
Parameters
----------
other : Rotation or list of Rotations.
other : damask.Rotation
"""
return self.copy(rotation=np.vstack(tuple(map(lambda x:x.quaternion,
@ -365,12 +365,28 @@ class Rotation:
def flatten(self,order = 'C'):
"""Flatten array."""
"""
Flatten array.
Returns
-------
flattened : damask.Rotation
Rotation flattened to single dimension.
"""
return self.copy(rotation=self.quaternion.reshape((-1,4),order=order))
def reshape(self,shape,order = 'C'):
"""Reshape array."""
"""
Reshape array.
Returns
-------
reshaped : damask.Rotation
Rotation of given shape.
"""
if isinstance(shape,(int,np.integer)): shape = (shape,)
return self.copy(rotation=self.quaternion.reshape(tuple(shape)+(4,),order=order))
@ -387,6 +403,11 @@ class Rotation:
Where to preferentially locate missing dimensions.
Either 'left' or 'right' (default).
Returns
-------
broadcasted : damask.Rotation
Rotation broadcasted to given shape.
"""
if isinstance(shape,(int,np.integer)): shape = (shape,)
return self.copy(rotation=np.broadcast_to(self.quaternion.reshape(util.shapeshifter(self.shape,shape,mode)+(4,)),
@ -404,7 +425,7 @@ class Rotation:
Returns
-------
average : Rotation
average : damask.Rotation
Weighted average of original Rotation field.
References
@ -438,9 +459,14 @@ class Rotation:
Parameters
----------
other : Rotation
other : damask.Rotation
Rotation to which the misorientation is computed.
Returns
-------
g : damask.Rotation
Misorientation.
"""
return other*~self

View File

@ -49,7 +49,7 @@ class Table:
Returns
-------
slice : Table
slice : damask.Table
Sliced part of the Table.
Examples
@ -157,7 +157,7 @@ class Table:
Parameters
----------
other : Table
other : damask.Table
Table to compare against.
rtol : float, optional
Relative tolerance of equality.
@ -185,7 +185,7 @@ class Table:
Parameters
----------
other : Table
other : damask.Table
Table to compare against.
rtol : float, optional
Relative tolerance of equality.
@ -223,6 +223,11 @@ class Table:
fname : file, str, or pathlib.Path
Filename or file for reading.
Returns
-------
loaded : damask.Table
Table data from file.
"""
try:
f = open(fname)
@ -275,6 +280,11 @@ class Table:
fname : file, str, or pathlib.Path
Filename or file for reading.
Returns
-------
loaded : damask.Table
Table data from file.
"""
try:
f = open(fname)
@ -334,14 +344,14 @@ class Table:
----------
label : str
Column label.
data : np.ndarray
data : numpy.ndarray
New data.
info : str, optional
Human-readable information about the new data.
Returns
-------
table : Table
updated : damask.Table
Updated table.
"""
@ -367,14 +377,14 @@ class Table:
----------
label : str
Column label.
data : np.ndarray
data : numpy.ndarray
Modified data.
info : str, optional
Human-readable information about the modified data.
Returns
-------
table : Table
udated : damask.Table
Updated table.
"""
@ -402,7 +412,7 @@ class Table:
Returns
-------
table : Table
udated : damask.Table
Updated table.
"""
@ -425,7 +435,7 @@ class Table:
Returns
-------
table : Table
udated : damask.Table
Updated table.
"""
@ -451,7 +461,7 @@ class Table:
Returns
-------
table : Table
udated : damask.Table
Updated table.
"""
@ -479,13 +489,13 @@ class Table:
Parameters
----------
other : Table
other : damask.Table
Table to append.
Returns
-------
table : Table
Concatenated table.
udated : damask.Table
Updated table.
"""
if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns):
@ -504,13 +514,13 @@ class Table:
Parameters
----------
other : Table
other : damask.Table
Table to join.
Returns
-------
table : Table
Joined table.
udated : damask.Table
Updated table.
"""
if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:

View File

@ -3,7 +3,6 @@ import multiprocessing as mp
from pathlib import Path
import numpy as np
import numpy.ma as ma
import vtk
from vtk.util.numpy_support import numpy_to_vtk as np_to_vtk
from vtk.util.numpy_support import numpy_to_vtkIdTypeArray as np_to_vtkIdTypeArray
@ -51,6 +50,11 @@ class VTK:
origin : iterable of float, len (3), optional
Spatial origin coordinates.
Returns
-------
new : damask.VTK
VTK-based geometry without nodal or cell data.
"""
vtk_data = vtk.vtkRectilinearGrid()
vtk_data.SetDimensions(*(np.array(grid)+1))
@ -68,7 +72,7 @@ class VTK:
"""
Create VTK of type vtk.vtkUnstructuredGrid.
This is the common type for FEM solver results.
This is the common type for mesh solver results.
Parameters
----------
@ -80,6 +84,11 @@ class VTK:
cell_type : str
Name of the vtk.vtkCell subclass. Tested for TRIANGLE, QUAD, TETRA, and HEXAHEDRON.
Returns
-------
new : damask.VTK
VTK-based geometry without nodal or cell data.
"""
vtk_nodes = vtk.vtkPoints()
vtk_nodes.SetData(np_to_vtk(nodes))
@ -110,6 +119,11 @@ class VTK:
points : numpy.ndarray of shape (:,3)
Spatial position of the points.
Returns
-------
new : damask.VTK
VTK-based geometry without nodal or cell data.
"""
N = points.shape[0]
vtk_points = vtk.vtkPoints()
@ -137,8 +151,13 @@ class VTK:
fname : str or pathlib.Path
Filename for reading. Valid extensions are .vtr, .vtu, .vtp, and .vtk.
dataset_type : str, optional
Name of the vtk.vtkDataSet subclass when opening a .vtk file. Valid types are vtkRectilinearGrid,
vtkUnstructuredGrid, and vtkPolyData.
Name of the vtk.vtkDataSet subclass when opening a .vtk file.
Valid types are vtkRectilinearGrid, vtkUnstructuredGrid, and vtkPolyData.
Returns
-------
loaded : damask.VTK
VTK-based geometry from file.
"""
if not os.path.isfile(fname): # vtk has a strange error handling
@ -149,13 +168,13 @@ class VTK:
reader.SetFileName(str(fname))
if dataset_type is None:
raise TypeError('Dataset type for *.vtk file not given.')
elif dataset_type.lower().endswith('rectilineargrid'):
elif dataset_type.lower().endswith(('rectilineargrid','rectilinear_grid')):
reader.Update()
vtk_data = reader.GetRectilinearGridOutput()
elif dataset_type.lower().endswith('unstructuredgrid'):
elif dataset_type.lower().endswith(('unstructuredgrid','unstructured_grid')):
reader.Update()
vtk_data = reader.GetUnstructuredGridOutput()
elif dataset_type.lower().endswith('polydata'):
elif dataset_type.lower().endswith(('polydata','poly_data')):
reader.Update()
vtk_data = reader.GetPolyDataOutput()
else:
@ -246,7 +265,7 @@ class VTK:
raise ValueError('No label defined for numpy.ndarray')
N_data = data.shape[0]
data_ = np.where(data.mask,data.fill_value,data) if isinstance(data,ma.MaskedArray) else\
data_ = np.where(data.mask,data.fill_value,data) if isinstance(data,np.ma.MaskedArray) else\
data
d = np_to_vtk((data_.astype(np.single) if data_.dtype in [np.double, np.longdouble] else
data_).reshape(N_data,-1),deep=True) # avoid large files
@ -277,6 +296,11 @@ class VTK:
label : str
Data label.
Returns
-------
data : numpy.ndarray
Data stored under the given label.
"""
cell_data = self.vtk_data.GetCellData()
for a in range(cell_data.GetNumberOfArrays()):

View File

@ -1,15 +1,14 @@
"""
Filters for operations on regular grids.
Notes
-----
The grids are defined as (x,y,z,...) where x is fastest and z is slowest.
This convention is consistent with the layout in grid vtr files.
When converting to/from a plain list (e.g. storage in ASCII table),
the following operations are required for tensorial data:
D3 = D1.reshape(cells+(-1,),order='F').reshape(cells+(3,3))
D1 = D3.reshape(cells+(-1,)).reshape(-1,9,order='F')
- D3 = D1.reshape(cells+(-1,),order='F').reshape(cells+(3,3))
- D1 = D3.reshape(cells+(-1,)).reshape(-1,9,order='F')
"""
from scipy import spatial as _spatial
@ -29,10 +28,12 @@ def _ks(size,cells,first_order=False):
Correction for first order derivatives, defaults to False.
"""
k_sk = _np.where(_np.arange(cells[0])>cells[0]//2,_np.arange(cells[0])-cells[0],_np.arange(cells[0]))/size[0]
k_sk = _np.where(_np.arange(cells[0])>cells[0]//2,
_np.arange(cells[0])-cells[0],_np.arange(cells[0]))/size[0]
if cells[0]%2 == 0 and first_order: k_sk[cells[0]//2] = 0 # Nyquist freq=0 for even cells (Johnson, MIT, 2011)
k_sj = _np.where(_np.arange(cells[1])>cells[1]//2,_np.arange(cells[1])-cells[1],_np.arange(cells[1]))/size[1]
k_sj = _np.where(_np.arange(cells[1])>cells[1]//2,
_np.arange(cells[1])-cells[1],_np.arange(cells[1]))/size[1]
if cells[1]%2 == 0 and first_order: k_sj[cells[1]//2] = 0 # Nyquist freq=0 for even cells (Johnson, MIT, 2011)
k_si = _np.arange(cells[2]//2+1)/size[2]
@ -40,74 +41,89 @@ def _ks(size,cells,first_order=False):
return _np.stack(_np.meshgrid(k_sk,k_sj,k_si,indexing = 'ij'), axis=-1)
def curl(size,field):
"""
def curl(size,f):
u"""
Calculate curl of a vector or tensor field in Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
f : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
Periodic field of which the curl is calculated.
Returns
-------
× f : numpy.ndarray
Curl of f.
"""
n = _np.prod(field.shape[3:])
k_s = _ks(size,field.shape[:3],True)
n = _np.prod(f.shape[3:])
k_s = _ks(size,f.shape[:3],True)
e = _np.zeros((3, 3, 3))
e[0, 1, 2] = e[1, 2, 0] = e[2, 0, 1] = +1.0 # Levi-Civita symbol
e[0, 2, 1] = e[2, 1, 0] = e[1, 0, 2] = -1.0
field_fourier = _np.fft.rfftn(field,axes=(0,1,2))
curl_ = (_np.einsum('slm,ijkl,ijkm ->ijks', e,k_s,field_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 3
_np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,field_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3x3
f_fourier = _np.fft.rfftn(f,axes=(0,1,2))
curl_ = (_np.einsum('slm,ijkl,ijkm ->ijks', e,k_s,f_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 3
_np.einsum('slm,ijkl,ijknm->ijksn',e,k_s,f_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3x3
return _np.fft.irfftn(curl_,axes=(0,1,2),s=field.shape[:3])
return _np.fft.irfftn(curl_,axes=(0,1,2),s=f.shape[:3])
def divergence(size,field):
"""
def divergence(size,f):
u"""
Calculate divergence of a vector or tensor field in Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
f : numpy.ndarray of shape (:,:,:,3) or (:,:,:,3,3)
Periodic field of which the divergence is calculated.
Returns
-------
· f : numpy.ndarray
Divergence of f.
"""
n = _np.prod(field.shape[3:])
k_s = _ks(size,field.shape[:3],True)
n = _np.prod(f.shape[3:])
k_s = _ks(size,f.shape[:3],True)
field_fourier = _np.fft.rfftn(field,axes=(0,1,2))
div_ = (_np.einsum('ijkl,ijkl ->ijk', k_s,field_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 1
_np.einsum('ijkm,ijklm->ijkl',k_s,field_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3
f_fourier = _np.fft.rfftn(f,axes=(0,1,2))
div_ = (_np.einsum('ijkl,ijkl ->ijk', k_s,f_fourier)*2.0j*_np.pi if n == 3 else # vector, 3 -> 1
_np.einsum('ijkm,ijklm->ijkl',k_s,f_fourier)*2.0j*_np.pi) # tensor, 3x3 -> 3
return _np.fft.irfftn(div_,axes=(0,1,2),s=field.shape[:3])
return _np.fft.irfftn(div_,axes=(0,1,2),s=f.shape[:3])
def gradient(size,field):
"""
Calculate gradient of a scalar or vector field in Fourier space.
def gradient(size,f):
u"""
Calculate gradient of a scalar or vector fieldin Fourier space.
Parameters
----------
size : numpy.ndarray of shape (3)
Physical size of the periodic field.
field : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
f : numpy.ndarray of shape (:,:,:,1) or (:,:,:,3)
Periodic field of which the gradient is calculated.
Returns
-------
f : numpy.ndarray
Divergence of f.
"""
n = _np.prod(field.shape[3:])
k_s = _ks(size,field.shape[:3],True)
n = _np.prod(f.shape[3:])
k_s = _ks(size,f.shape[:3],True)
field_fourier = _np.fft.rfftn(field,axes=(0,1,2))
grad_ = (_np.einsum('ijkl,ijkm->ijkm', field_fourier,k_s)*2.0j*_np.pi if n == 1 else # scalar, 1 -> 3
_np.einsum('ijkl,ijkm->ijklm',field_fourier,k_s)*2.0j*_np.pi) # vector, 3 -> 3x3
f_fourier = _np.fft.rfftn(f,axes=(0,1,2))
grad_ = (_np.einsum('ijkl,ijkm->ijkm', f_fourier,k_s)*2.0j*_np.pi if n == 1 else # scalar, 1 -> 3
_np.einsum('ijkl,ijkm->ijklm',f_fourier,k_s)*2.0j*_np.pi) # vector, 3 -> 3x3
return _np.fft.irfftn(grad_,axes=(0,1,2),s=field.shape[:3])
return _np.fft.irfftn(grad_,axes=(0,1,2),s=f.shape[:3])
def coordinates0_point(cells,size,origin=_np.zeros(3)):
@ -123,6 +139,11 @@ def coordinates0_point(cells,size,origin=_np.zeros(3)):
origin : numpy.ndarray, optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_p_0 : numpy.ndarray
Undeformed cell center coordinates.
"""
start = origin + size/cells*.5
end = origin + size - size/cells*.5
@ -144,6 +165,11 @@ def displacement_fluct_point(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_p_fluct : numpy.ndarray
Fluctuating part of the cell center displacements.
"""
integrator = 0.5j*size/_np.pi
@ -171,6 +197,11 @@ def displacement_avg_point(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_p_avg : numpy.ndarray
Average part of the cell center displacements.
"""
F_avg = _np.average(F,axis=(0,1,2))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_point(F.shape[:3],size))
@ -187,6 +218,11 @@ def displacement_point(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_p : numpy.ndarray
Cell center displacements.
"""
return displacement_avg_point(size,F) + displacement_fluct_point(size,F)
@ -204,6 +240,11 @@ def coordinates_point(size,F,origin=_np.zeros(3)):
origin : numpy.ndarray of shape (3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_p : numpy.ndarray
Cell center coordinates.
"""
return coordinates0_point(F.shape[:3],size,origin) + displacement_point(size,F)
@ -220,6 +261,11 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True):
Expect coordinates0 data to be ordered (x fast, z slow).
Defaults to True.
Returns
-------
cells, size, origin : Three numpy.ndarray, each of shape (3)
Information to reconstruct grid.
"""
coords = [_np.unique(coordinates0[:,i]) for i in range(3)]
mincorner = _np.array(list(map(min,coords)))
@ -252,19 +298,6 @@ def cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True):
return (cells,size,origin)
def coordinates0_check(coordinates0):
"""
Check whether coordinates lie on a regular grid.
Parameters
----------
coordinates0 : numpy.ndarray
Array of undeformed cell coordinates.
"""
cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True)
def coordinates0_node(cells,size,origin=_np.zeros(3)):
"""
Nodal positions (undeformed).
@ -278,6 +311,11 @@ def coordinates0_node(cells,size,origin=_np.zeros(3)):
origin : numpy.ndarray of shape (3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_n_0 : numpy.ndarray
Undeformed nodal coordinates.
"""
return _np.stack(_np.meshgrid(_np.linspace(origin[0],size[0]+origin[0],cells[0]+1),
_np.linspace(origin[1],size[1]+origin[1],cells[1]+1),
@ -296,6 +334,11 @@ def displacement_fluct_node(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_n_fluct : numpy.ndarray
Fluctuating part of the nodal displacements.
"""
return point_to_node(displacement_fluct_point(size,F))
@ -311,6 +354,11 @@ def displacement_avg_node(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_n_avg : numpy.ndarray
Average part of the nodal displacements.
"""
F_avg = _np.average(F,axis=(0,1,2))
return _np.einsum('ml,ijkl->ijkm',F_avg - _np.eye(3),coordinates0_node(F.shape[:3],size))
@ -327,6 +375,11 @@ def displacement_node(size,F):
F : numpy.ndarray
Deformation gradient field.
Returns
-------
u_p : numpy.ndarray
Nodal displacements.
"""
return displacement_avg_node(size,F) + displacement_fluct_node(size,F)
@ -344,28 +397,15 @@ def coordinates_node(size,F,origin=_np.zeros(3)):
origin : numpy.ndarray of shape (3), optional
Physical origin of the periodic field. Defaults to [0.0,0.0,0.0].
Returns
-------
x_n : numpy.ndarray
Nodal coordinates.
"""
return coordinates0_node(F.shape[:3],size,origin) + displacement_node(size,F)
def point_to_node(cell_data):
"""Interpolate periodic point data to nodal data."""
n = ( cell_data + _np.roll(cell_data,1,(0,1,2))
+ _np.roll(cell_data,1,(0,)) + _np.roll(cell_data,1,(1,)) + _np.roll(cell_data,1,(2,))
+ _np.roll(cell_data,1,(0,1)) + _np.roll(cell_data,1,(1,2)) + _np.roll(cell_data,1,(2,0)))*0.125
return _np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
def node_2_point(node_data):
"""Interpolate periodic nodal data to point data."""
c = ( node_data + _np.roll(node_data,1,(0,1,2))
+ _np.roll(node_data,1,(0,)) + _np.roll(node_data,1,(1,)) + _np.roll(node_data,1,(2,))
+ _np.roll(node_data,1,(0,1)) + _np.roll(node_data,1,(1,2)) + _np.roll(node_data,1,(2,0)))*0.125
return c[1:,1:,1:]
def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
"""
Return grid 'DNA', i.e. cells, size, and origin from 1D array of nodal positions.
@ -378,6 +418,11 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
Expect coordinates0 data to be ordered (x fast, z slow).
Defaults to True.
Returns
-------
cells, size, origin : Three numpy.ndarray, each of shape (3)
Information to reconstruct grid.
"""
coords = [_np.unique(coordinates0[:,i]) for i in range(3)]
mincorner = _np.array(list(map(min,coords)))
@ -402,6 +447,72 @@ def cellsSizeOrigin_coordinates0_node(coordinates0,ordered=True):
return (cells,size,origin)
def point_to_node(cell_data):
"""
Interpolate periodic point data to nodal data.
Parameters
----------
cell_data : numpy.ndarray of shape (:,:,:,...)
Data defined on the cell centers of a periodic grid.
Returns
-------
node_data : numpy.ndarray of shape (:,:,:,...)
Data defined on the nodes of a periodic grid.
"""
n = ( cell_data + _np.roll(cell_data,1,(0,1,2))
+ _np.roll(cell_data,1,(0,)) + _np.roll(cell_data,1,(1,)) + _np.roll(cell_data,1,(2,))
+ _np.roll(cell_data,1,(0,1)) + _np.roll(cell_data,1,(1,2)) + _np.roll(cell_data,1,(2,0)))*0.125
return _np.pad(n,((0,1),(0,1),(0,1))+((0,0),)*len(cell_data.shape[3:]),mode='wrap')
def node_to_point(node_data):
"""
Interpolate periodic nodal data to point data.
Parameters
----------
node_data : numpy.ndarray of shape (:,:,:,...)
Data defined on the nodes of a periodic grid.
Returns
-------
cell_data : numpy.ndarray of shape (:,:,:,...)
Data defined on the cell centers of a periodic grid.
"""
c = ( node_data + _np.roll(node_data,1,(0,1,2))
+ _np.roll(node_data,1,(0,)) + _np.roll(node_data,1,(1,)) + _np.roll(node_data,1,(2,))
+ _np.roll(node_data,1,(0,1)) + _np.roll(node_data,1,(1,2)) + _np.roll(node_data,1,(2,0)))*0.125
return c[1:,1:,1:]
def coordinates0_valid(coordinates0):
"""
Check whether coordinates lie on a regular grid.
Parameters
----------
coordinates0 : numpy.ndarray
Array of undeformed cell coordinates.
Returns
-------
valid : bool
Wheter the coordinates lie on a regular grid.
"""
try:
cellsSizeOrigin_coordinates0_point(coordinates0,ordered=True)
return True
except ValueError:
return False
def regrid(size,F,cells):
"""
Return mapping from coordinates in deformed configuration to a regular grid.

View File

@ -1,4 +1,9 @@
"""Finite-strain continuum mechanics."""
"""
Finite-strain continuum mechanics.
All routines operate on numpy.ndarrays of shape (...,3,3).
"""
from . import tensor as _tensor
from . import _rotation
@ -154,7 +159,6 @@ def strain(F,t,m):
return eps
def stress_Cauchy(P,F):
"""
Calculate the Cauchy stress (true stress).

View File

@ -24,6 +24,11 @@ def from_random(size,N_seeds,cells=None,rng_seed=None):
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
Returns
-------
coords : numpy.ndarray of shape (N_seeds,3)
Seed coordinates in 3D space.
"""
rng = _np.random.default_rng(rng_seed)
if cells is None:
@ -56,6 +61,11 @@ def from_Poisson_disc(size,N_seeds,N_candidates,distance,periodic=True,rng_seed=
A seed to initialize the BitGenerator. Defaults to None.
If None, then fresh, unpredictable entropy will be pulled from the OS.
Returns
-------
coords : numpy.ndarray of shape (N_seeds,3)
Seed coordinates in 3D space.
"""
rng = _np.random.default_rng(rng_seed)
coords = _np.empty((N_seeds,3))
@ -94,6 +104,11 @@ def from_grid(grid,selection=None,invert=False,average=False,periodic=True):
periodic : boolean, optional
Center of gravity with periodic boundaries.
Returns
-------
coords, materials : numpy.ndarray of shape (:,3), numpy.ndarray of shape (:)
Seed coordinates in 3D space, material IDs.
"""
material = grid.material.reshape((-1,1),order='F')
mask = _np.full(grid.cells.prod(),True,dtype=bool) if selection is None else \

View File

@ -1,3 +1,3 @@
"""Tools to control the various solvers."""
"""Run simulations directly from python."""
from ._marc import Marc # noqa

View File

@ -1,10 +1,7 @@
"""
Tensor operations.
Tensor mathematics.
Notes
-----
This is not a tensor class, but a collection of routines
to operate on numpy.ndarrays of shape (...,3,3).
All routines operate on numpy.ndarrays of shape (...,3,3).
"""

View File

@ -49,7 +49,7 @@ _colors = {
####################################################################################################
def srepr(arg,glue = '\n'):
r"""
Join arguments with glue string.
Join items with glue string.
Parameters
----------
@ -58,6 +58,11 @@ def srepr(arg,glue = '\n'):
glue : str, optional
Glue used for joining operation. Defaults to \n.
Returns
-------
joined : str
String representation of the joined items.
"""
if (not hasattr(arg, 'strip') and
(hasattr(arg, '__getitem__') or
@ -68,19 +73,71 @@ def srepr(arg,glue = '\n'):
def emph(what):
"""Formats string with emphasis."""
"""
Format with emphasis.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
return _colors['bold']+srepr(what)+_colors['end_color']
def deemph(what):
"""Formats string with deemphasis."""
"""
Format with deemphasis.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
return _colors['dim']+srepr(what)+_colors['end_color']
def warn(what):
"""Formats string for warning."""
"""
Format for warning.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
return _colors['warning']+emph(what)+_colors['end_color']
def strikeout(what):
"""Formats string as strikeout."""
"""
Format as strikeout.
Parameters
----------
what : object with __repr__ or iterable of objects with __repr__.
Message to format.
Returns
-------
formatted : str
Formatted string representation of the joined items.
"""
return _colors['crossout']+srepr(what)+_colors['end_color']
@ -97,6 +154,11 @@ def execute(cmd,wd='./',env=None):
env : dict, optional
Environment for execution.
Returns
-------
stdout, stderr : str
Output of the executed command.
"""
print(f"executing '{cmd}' in '{wd}'")
process = subprocess.run(shlex.split(cmd),
@ -157,6 +219,11 @@ def scale_to_coprime(v):
v : numpy.ndarray of shape (:)
Vector to scale.
Returns
-------
m : numpy.ndarray of shape (:)
Vector scaled to co-prime numbers.
"""
MAX_DENOMINATOR = 1000000
@ -371,9 +438,14 @@ def DREAM3D_base_group(fname):
Parameters
----------
fname : str
fname : str or pathlib.Path
Filename of the DREAM.3D (HDF5) file.
Returns
-------
path : str
Path to the base group.
"""
with h5py.File(fname,'r') as f:
base_group = f.visit(lambda path: path.rsplit('/',2)[0] if '_SIMPL_GEOMETRY/SPACING' in path else None)
@ -393,9 +465,14 @@ def DREAM3D_cell_data_group(fname):
Parameters
----------
fname : str
fname : str or pathlib.Path
Filename of the DREAM.3D (HDF5) file.
Returns
-------
path : str
Path to the cell data group.
"""
base_group = DREAM3D_base_group(fname)
with h5py.File(fname,'r') as f:

View File

@ -271,7 +271,7 @@ class TestResult:
@pytest.mark.parametrize('overwrite',['off','on'])
def test_add_overwrite(self,default,overwrite):
last = default.view('times',default.times_in_range(0,np.inf)[-1])
last = default.view('increments',-1)
last.add_stress_Cauchy()
@ -279,9 +279,9 @@ class TestResult:
created_first = datetime.strptime(created_first,'%Y-%m-%d %H:%M:%S%z')
if overwrite == 'on':
last = last.allow_modification()
last = last.modification_enable()
else:
last = last.disallow_modification()
last = last.modification_disable()
time.sleep(2.)
try:
@ -301,14 +301,24 @@ class TestResult:
def test_rename(self,default,allowed):
if allowed == 'on':
F = default.place('F')
default = default.allow_modification()
default = default.modification_enable()
default.rename('F','new_name')
assert np.all(F == default.place('new_name'))
default = default.disallow_modification()
default = default.modification_disable()
with pytest.raises(PermissionError):
default.rename('P','another_new_name')
@pytest.mark.parametrize('allowed',['off','on'])
def test_remove(self,default,allowed):
if allowed == 'on':
unsafe = default.modification_enable()
unsafe.remove('F')
assert unsafe.get('F') is None
else:
with pytest.raises(PermissionError):
default.remove('F')
@pytest.mark.parametrize('mode',['cell','node'])
def test_coordinates(self,default,mode):
if mode == 'cell':

View File

@ -60,7 +60,7 @@ class TestGridFilters:
cell_field_x = np.interp(coordinates0_point_x,coordinates_node_x,node_field_x,period=np.pi*2.)
cell_field = np.broadcast_to(cell_field_x.reshape(-1,1,1),cells)
assert np.allclose(cell_field,grid_filters.node_2_point(node_field))
assert np.allclose(cell_field,grid_filters.node_to_point(node_field))
@pytest.mark.parametrize('mode',['point','node'])
def test_coordinates0_origin(self,mode):
@ -93,14 +93,24 @@ class TestGridFilters:
F = np.broadcast_to(np.random.random((3,3)), tuple(cells)+(3,3))
assert np.allclose(function(size,F),0.0)
@pytest.mark.parametrize('function',[grid_filters.coordinates0_check,
grid_filters.cellsSizeOrigin_coordinates0_node,
grid_filters.cellsSizeOrigin_coordinates0_point])
@pytest.mark.parametrize('function',[grid_filters.cellsSizeOrigin_coordinates0_point,
grid_filters.cellsSizeOrigin_coordinates0_node])
def test_invalid_coordinates(self,function):
invalid_coordinates = np.random.random((np.random.randint(12,52),3))
with pytest.raises(ValueError):
function(invalid_coordinates)
@pytest.mark.parametrize('function',[grid_filters.coordinates0_point,
grid_filters.coordinates0_node])
def test_valid_coordinates_check(self,function):
valid_coordinates = function(np.random.randint(4,10,(3)),np.random.rand(3))
assert grid_filters.coordinates0_valid(valid_coordinates.reshape(-1,3,order='F'))
def test_invalid_coordinates_check(self):
invalid_coordinates = np.random.random((np.random.randint(12,52),3))
assert not grid_filters.coordinates0_valid(invalid_coordinates)
@pytest.mark.parametrize('function',[grid_filters.cellsSizeOrigin_coordinates0_node,
grid_filters.cellsSizeOrigin_coordinates0_point])
def test_uneven_spaced_coordinates(self,function):

View File

@ -40,7 +40,6 @@
#include "phase_damage_isobrittle.f90"
#include "phase_damage_isoductile.f90"
#include "phase_damage_anisobrittle.f90"
#include "phase_damage_anisoductile.f90"
#include "homogenization.f90"
#include "homogenization_mechanical.f90"
#include "homogenization_mechanical_pass.f90"

View File

@ -35,7 +35,7 @@ module material
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt, & !< phase ID of each element
material_phaseID, & !< per (constituent,cell)
material_phaseEntry !< per (constituent,cell
material_phaseEntry !< per (constituent,cell)
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem)
material_phaseMemberAt !< position of the element within its phase instance

View File

@ -12,8 +12,7 @@ submodule(phase) damage
DAMAGE_UNDEFINED_ID, &
DAMAGE_ISOBRITTLE_ID, &
DAMAGE_ISODUCTILE_ID, &
DAMAGE_ANISOBRITTLE_ID, &
DAMAGE_ANISODUCTILE_ID
DAMAGE_ANISOBRITTLE_ID
end enum
@ -34,10 +33,6 @@ submodule(phase) damage
logical, dimension(:), allocatable :: mySources
end function anisobrittle_init
module function anisoductile_init() result(mySources)
logical, dimension(:), allocatable :: mySources
end function anisoductile_init
module function isobrittle_init() result(mySources)
logical, dimension(:), allocatable :: mySources
end function isobrittle_init
@ -62,10 +57,6 @@ submodule(phase) damage
S
end subroutine anisobrittle_dotState
module subroutine anisoductile_dotState(ph,me)
integer, intent(in) :: ph,me
end subroutine anisoductile_dotState
module subroutine isoductile_dotState(ph,me)
integer, intent(in) :: ph,me
end subroutine isoductile_dotState
@ -75,11 +66,6 @@ submodule(phase) damage
character(len=*), intent(in) :: group
end subroutine anisobrittle_results
module subroutine anisoductile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine anisoductile_results
module subroutine isobrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
@ -146,7 +132,6 @@ module subroutine damage_init
where(isobrittle_init() ) phase_source = DAMAGE_ISOBRITTLE_ID
where(isoductile_init() ) phase_source = DAMAGE_ISODUCTILE_ID
where(anisobrittle_init()) phase_source = DAMAGE_ANISOBRITTLE_ID
where(anisoductile_init()) phase_source = DAMAGE_ANISODUCTILE_ID
endif
end subroutine damage_init
@ -171,7 +156,7 @@ module function phase_f_phi(phi,co,ce) result(f)
en = material_phaseEntry(co,ce)
select case(phase_source(ph))
case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID,DAMAGE_ANISODUCTILE_ID)
case(DAMAGE_ISOBRITTLE_ID,DAMAGE_ISODUCTILE_ID,DAMAGE_ANISOBRITTLE_ID)
f = 1.0_pReal &
- phi*damageState(ph)%state(1,en)
case default
@ -304,9 +289,6 @@ module subroutine damage_results(group,ph)
case (DAMAGE_ANISOBRITTLE_ID) sourceType
call anisobrittle_results(ph,group//'damage/')
case (DAMAGE_ANISODUCTILE_ID) sourceType
call anisoductile_results(ph,group//'damage/')
end select sourceType
end subroutine damage_results
@ -332,9 +314,6 @@ function phase_damage_collectDotState(ph,me) result(broken)
case (DAMAGE_ISODUCTILE_ID) sourceType
call isoductile_dotState(ph,me)
case (DAMAGE_ANISODUCTILE_ID) sourceType
call anisoductile_dotState(ph,me)
case (DAMAGE_ANISOBRITTLE_ID) sourceType
call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress?

View File

@ -1,138 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incorporating anisotropic ductile damage source mechanism
!> @details to be done
!--------------------------------------------------------------------------------------------------
submodule(phase:damage) anisoductile
type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: &
q !< damage rate sensitivity
real(pReal), dimension(:), allocatable :: &
gamma_crit !< critical plastic strain per slip system
character(len=pStringLen), allocatable, dimension(:) :: &
output
end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
module function anisoductile_init() result(mySources)
logical, dimension(:), allocatable :: mySources
class(tNode), pointer :: &
phases, &
phase, &
mech, &
pl, &
sources, &
src
integer :: Ninstances,Nmembers,ph
integer, dimension(:), allocatable :: N_sl
character(len=pStringLen) :: extmsg = ''
mySources = source_active('anisoductile')
if(count(mySources) == 0) return
print'(/,a)', ' <<<+- phase:damage:anisoductile init -+>>>'
print'(a,i0)', ' # phases: ',count(mySources); flush(IO_STDOUT)
phases => config_material%get('phase')
allocate(param(phases%length))
do ph = 1, phases%length
if(mySources(ph)) then
phase => phases%get(ph)
mech => phase%get('mechanical')
pl => mech%get('plastic')
sources => phase%get('damage')
associate(prm => param(ph))
src => sources%get(1)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
prm%q = src%get_asFloat('q')
prm%gamma_crit = src%get_as1dFloat('gamma_crit',requiredSize=size(N_sl))
! expand: family => system
prm%gamma_crit = math_expand(prm%gamma_crit,N_sl)
#if defined (__GFORTRAN__)
prm%output = output_as1dString(src)
#else
prm%output = src%get_as1dString('output',defaultVal=emptyStringArray)
#endif
! sanity checks
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
Nmembers=count(material_phaseID==ph)
call phase_allocateState(damageState(ph),Nmembers,1,1,0)
damageState(ph)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
end associate
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoDuctile)')
endif
enddo
end function anisoductile_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!--------------------------------------------------------------------------------------------------
module subroutine anisoductile_dotState(ph,me)
integer, intent(in) :: &
ph, &
me
associate(prm => param(ph))
damageState(ph)%dotState(1,me) = sum(plasticState(ph)%slipRate(:,me)/(damage_phi(ph,me)**prm%q)/prm%gamma_crit)
end associate
end subroutine anisoductile_dotState
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
module subroutine anisoductile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
integer :: o
associate(prm => param(phase), stt => damageState(phase)%state)
outputsLoop: do o = 1,size(prm%output)
select case(trim(prm%output(o)))
case ('f_phi')
call results_writeDataset(group,stt,trim(prm%output(o)),'driving force','J/m³')
end select
enddo outputsLoop
end associate
end subroutine anisoductile_results
end submodule anisoductile

View File

@ -190,7 +190,7 @@ end function phase_thermal_collectDotState
!--------------------------------------------------------------------------------------------------
!> @brief Damage viscosity.
!> @brief Thermal viscosity.
!--------------------------------------------------------------------------------------------------
module function phase_mu_T(co,ce) result(mu)
@ -205,7 +205,7 @@ end function phase_mu_T
!--------------------------------------------------------------------------------------------------
!> @brief Damage conductivity/diffusivity in reference configuration.
!> @brief Thermal conductivity/diffusivity in reference configuration.
!--------------------------------------------------------------------------------------------------
module function phase_K_T(co,ce) result(K)