223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

outsourcing tensor math to mechanics class 

strain calculation is generalize to arbitrary order and simplified: No
need for svd, F^T F/F F^T does the job.
This commit is contained in:
Martin Diehl 2019-10-19 12:54:16 +02:00
parent e51f6cee72
commit b31de5d0f6
2 changed files with 324 additions and 195 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

217
python/damask/dadf5.py
View File

@ -360,7 +360,7 @@ class DADF5():
def cell_coordinates(self): def cell_coordinates(self):
"""Initial coordinates of the cell centers.""" """Return initial coordinates of the cell centers."""
if self.structured: if self.structured:
delta = self.size/self.grid*0.5 delta = self.size/self.grid*0.5
z, y, x = np.meshgrid(np.linspace(delta[2],self.size[2]-delta[2],self.grid[2]), z, y, x = np.meshgrid(np.linspace(delta[2],self.size[2]-delta[2],self.grid[2]),
@ -375,11 +375,17 @@ class DADF5():
def add_Cauchy(self,P='P',F='F'): def add_Cauchy(self,P='P',F='F'):
""" """
Adds Cauchy stress calculated from 1st Piola-Kirchhoff stress and deformation gradient. Add Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
Resulting tensor is symmetrized as the Cauchy stress should be symmetric. Parameters
----------
P : str, optional
Label of the dataset containing the 1. Piola-Kirchhoff stress. Default value is ‘P’.
F : str, optional
Label of the dataset containing the deformation gradient. Default value is ‘F’.
""" """
def Cauchy(F,P): def __add_Cauchy(F,P):
return { return {
'data': mechanics.Cauchy(F['data'],P['data']), 'data': mechanics.Cauchy(F['data'],P['data']),
'label': 'sigma', 'label': 'sigma',
@ -394,43 +400,47 @@ class DADF5():
requested = [{'label':F,'arg':'F'}, requested = [{'label':F,'arg':'F'},
{'label':P,'arg':'P'} ] {'label':P,'arg':'P'} ]
self.__add_generic_pointwise(Cauchy,requested) self.__add_generic_pointwise(__add_Cauchy,requested)
def add_Mises(self,x): def add_Mises(self,x):
"""Adds the equivalent Mises stress or strain of a tensor.""" """
def Mises(x): Add the equivalent Mises stress or strain of a symmetric tensor.
if x['meta']['Unit'] == 'Pa': #ToDo: Should we use this? Then add_Cauchy and add_strain_tensors also should perform sanity checks Parameters
t = 'stress' ----------
elif x['meta']['Unit'] == '1': x : str
t = 'strain' Label of the dataset containing a symmetric stress or strain tensor
else: """
print(x['meta']['Unit']) def __add_Mises(x):
raise ValueError
return { return {
'data' : mechanics.Mises_strain(x['data']) if t=='strain' else mechanics.Mises_stress(x['data']), 'data': mechanics.Mises_strain(x) if t=='strain' else mechanics.Mises_stress(x),
'label': '{}_vM'.format(x['label']), 'label': '{}_vM'.format(x['label']),
'meta': { 'meta': {
'Unit': x['meta']['Unit'], 'Unit': x['meta']['Unit'],
'Description': 'Mises equivalent {} of {} ({})'.format(t,x['label'],x['meta']['Description']), 'Description': 'Mises equivalent {} of {} ({})'.format(t,x['label'],x['meta']['Description']),
'Creator' : 'dadf5.py:add_Mises_stress v{}'.format(version) 'Creator': 'dadf5.py:add_Mises v{}'.format(version)
} }
} }
requested = [{'label':x,'arg':'x'}] requested = [{'label':x,'arg':'x'}]
self.__add_generic_pointwise(Mises,requested) self.__add_generic_pointwise(__add_Mises,requested)
def add_norm(self,x,ord=None): def add_norm(self,x,ord=None):
""" """
Adds norm of vector or tensor. Add the norm of vector or tensor.
See numpy.linalg.norm manual for details. Parameters
----------
x : str
Label of the dataset containing a vector or tensor.
ord : {non-zero int, inf, -inf, ‘fro’, ‘nuc’}, optional
Order of the norm. inf means numpy’s inf object. For details refer to numpy.linalg.norm.
""" """
def norm(x,ord): def __add_norm(x,ord):
o = ord o = ord
if len(x['data'].shape) == 2: if len(x['data'].shape) == 2:
@ -456,12 +466,19 @@ class DADF5():
requested = [{'label':x,'arg':'x'}] requested = [{'label':x,'arg':'x'}]
self.__add_generic_pointwise(norm,requested,{'ord':ord}) self.__add_generic_pointwise(__add_norm,requested,{'ord':ord})
def add_absolute(self,x): def add_absolute(self,x):
"""Adds absolute value.""" """
def absolute(x): Add absolute value.
Parameters
----------
x : str
Label of the dataset containing a scalar, vector, or tensor.
"""
def __add_absolute(x):
return { return {
'data': np.abs(x['data']), 'data': np.abs(x['data']),
@ -475,12 +492,19 @@ class DADF5():
requested = [{'label':x,'arg':'x'}] requested = [{'label':x,'arg':'x'}]
self.__add_generic_pointwise(absolute,requested) self.__add_generic_pointwise(__add_absolute,requested)
def add_determinant(self,x): def add_determinant(self,x):
"""Adds the determinant component of a tensor.""" """
def determinant(x): Add the determinant of a tensor.
Parameters
----------
x : str
Label of the dataset containing a tensor.
"""
def __add_determinant(x):
return { return {
'data': np.linalg.det(x['data']), 'data': np.linalg.det(x['data']),
@ -494,19 +518,26 @@ class DADF5():
requested = [{'label':x,'arg':'x'}] requested = [{'label':x,'arg':'x'}]
self.__add_generic_pointwise(determinant,requested) self.__add_generic_pointwise(__add_determinant,requested)
def add_spherical(self,x): def add_spherical(self,x):
"""Adds the spherical component of a tensor.""" """
def spherical(x): Add the spherical (hydrostatic) part of a tensor.
Parameters
----------
x : str
Label of the dataset containing a tensor.
"""
def __add_spherical(x):
if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])): if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])):
raise ValueError raise ValueError
return { return {
'data' : np.trace(x['data'],axis1=1,axis2=2)/3.0, 'data': mechanics.spherical_part(x),
'label' : 'sph({})'.format(x['label']), 'label': 'p_{}'.format(x['label']),
'meta': { 'meta': {
'Unit': x['meta']['Unit'], 'Unit': x['meta']['Unit'],
'Description': 'Spherical component of tensor {} ({})'.format(x['label'],x['meta']['Description']), 'Description': 'Spherical component of tensor {} ({})'.format(x['label'],x['meta']['Description']),
@ -516,19 +547,26 @@ class DADF5():
requested = [{'label':x,'arg':'x'}] requested = [{'label':x,'arg':'x'}]
self.__add_generic_pointwise(spherical,requested) self.__add_generic_pointwise(__add_spherical,requested)
def add_deviator(self,x): def add_deviator(self,x):
"""Adds the deviator of a tensor.""" """
def deviator(x): Add the deviatoric part of a tensor.
Parameters
----------
x : str
Label of the dataset containing a tensor.
"""
def __add_deviator(x):
if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])): if not np.all(np.array(x['data'].shape[1:]) == np.array([3,3])):
raise ValueError raise ValueError
return { return {
'data': mechanics.deviator(x['data']), 'data': mechanics.deviator(x['data']),
'label' : 'dev({})'.format(x['label']), 'label': 's_{}'.format(x['label']),
'meta': { 'meta': {
'Unit': x['meta']['Unit'], 'Unit': x['meta']['Unit'],
'Description': 'Deviator of tensor {} ({})'.format(x['label'],x['meta']['Description']), 'Description': 'Deviator of tensor {} ({})'.format(x['label'],x['meta']['Description']),
@ -538,20 +576,30 @@ class DADF5():
requested = [{'label':x,'arg':'x'}] requested = [{'label':x,'arg':'x'}]
self.__add_generic_pointwise(deviator,requested) self.__add_generic_pointwise(__add_deviator,requested)
def add_calculation(self,formula,label,unit='n/a',description=None,vectorized=True): def add_calculation(self,formula,label,unit='n/a',description=None,vectorized=True):
""" """
General formula. Add result of a general formula.
Works currently only for vectorized expressions
Parameters
----------
formula : str
Formula, refer to datasets by ‘#Label#‘.
label : str
Label of the dataset containing the result of the calculation.
unit : str, optional
Physical unit of the result.
description : str, optional
Human readable description of the result.
vectorized : bool, optional
Indicate whether the formula is written in vectorized form.
""" """
if vectorized is not True: if vectorized is not True:
raise NotImplementedError raise NotImplementedError
def calculation(**kwargs): def __add_calculation(**kwargs):
formula = kwargs['formula'] formula = kwargs['formula']
for d in re.findall(r'#(.*?)#',formula): for d in re.findall(r'#(.*?)#',formula):
@ -570,62 +618,57 @@ class DADF5():
requested = [{'label':d,'arg':d} for d in set(re.findall(r'#(.*?)#',formula))] # datasets used in the formula requested = [{'label':d,'arg':d} for d in set(re.findall(r'#(.*?)#',formula))] # datasets used in the formula
pass_through = {'formula':formula,'label':label,'unit':unit,'description':description} pass_through = {'formula':formula,'label':label,'unit':unit,'description':description}
self.__add_generic_pointwise(calculation,requested,pass_through) self.__add_generic_pointwise(__add_calculation,requested,pass_through)
def add_strain_tensor(self,t,ord,defgrad='F'): def add_strain_tensor(self,F='F',t='U',ord=0):
""" """
Adds the a strain tensor. Add strain tensor calculated from a deformation gradient.
Albrecht Bertram: Elasticity and Plasticity of Large Deformations An Introduction (3rd Edition, 2012), p. 102. For details refer to damask.mechanics.strain_tensor
Parameters
----------
F : str, optional
Label of the dataset containing the deformation gradient. Default value is ‘F’.
t : {‘V’, ‘U’}, optional
Type of the polar decomposition, ‘V’ for right stretch tensor and ‘U’ for left stretch tensor.
Defaults value is ‘U’.
ord : float, optional
Order of the strain calculation. Default value is ‘0.0’.
""" """
def strain_tensor(defgrad,t,ord): def __add_strain_tensor(F,t,ord):
operator = {
'V#ln': lambda V: np.log(V),
'U#ln': lambda U: np.log(U),
'V#Biot': lambda V: np.broadcast_to(np.ones(3),[V.shape[0],3]) - 1.0/V,
'U#Biot': lambda U: U - np.broadcast_to(np.ones(3),[U.shape[0],3]),
'V#Green':lambda V: np.broadcast_to(np.ones(3),[V.shape[0],3]) - 1.0/V**2,
'U#Green':lambda U: U**2 - np.broadcast_to(np.ones(3),[U.shape[0],3]),
}
if t.lower() in ['l','left']:
stretch = 'V'
elif t.lower() in ['r','right']:
stretch = 'U'
else:
raise KeyError
(U,S,Vh) = np.linalg.svd(defgrad['data']) # singular value decomposition
R_inv = np.transpose(np.matmul(U,Vh),(0,2,1)) # transposed rotation of polar decomposition
s = np.matmul(R_inv,defgrad['data']) if stretch == 'U' else \
np.matmul(defgrad['data'],R_inv) # compute either left or right stretch
(D,V) = np.linalg.eigh((s+np.transpose(s,(0,2,1)))*.5) # eigen decomposition (of symmetric(ed) matrix)
d = operator[stretch+'#'+{0:'ln',1:'Biot',2:'Green'}[ord]](D)
a = np.matmul(V,np.einsum('ij,ikj->ijk',d,V))
return { return {
'data' : a, 'data': mechanics.strain_tensor(F['data'],t,ord),
'label' : 'epsilon_{}^{}({})'.format(stretch,ord,defgrad['label']), 'label': 'epsilon_{}^{}({})'.format(t,ord,F['label']),
'meta': { 'meta': {
'Unit' : defgrad['meta']['Unit'], 'Unit': F['meta']['Unit'],
'Description' : 'Strain tensor of {} ({})'.format(defgrad['label'],defgrad['meta']['Description']), 'Description': 'Strain tensor of {} ({})'.format(F['label'],F['meta']['Description']),
'Creator': 'dadf5.py:add_strain_tensor v{}'.format(version) 'Creator': 'dadf5.py:add_strain_tensor v{}'.format(version)
} }
} }
requested = [{'label':defgrad,'arg':'defgrad'}] requested = [{'label':F,'arg':'F'}]
self.__add_generic_pointwise(strain_tensor,requested,{'t':t,'ord':ord}) self.__add_generic_pointwise(__add_strain_tensor,requested,{'t':t,'ord':ord})
def add_principal_components(self,x): def add_principal_components(self,x):
"""Adds principal components of symmetric tensor.""" """
def principal_components(x): Add principal components of symmetric tensor.
The principal components are sorted in descending order, each repeated according to its multiplicity.
Parameters
----------
x : str
Label of the dataset containing a symmetric tensor.
"""
def __add_principal_components(x):
return { return {
'data' : np.linalg.eigvalsh((x['data']+np.transpose(x['data'],(0,2,1)))*.5)[:,::-1], # eigenvalues (of symmetric(ed) matrix) 'data': mechanics.principal_components(x),
'label': 'lambda_{}'.format(x['label']), 'label': 'lambda_{}'.format(x['label']),
'meta': { 'meta': {
'Unit': x['meta']['Unit'], 'Unit': x['meta']['Unit'],
@ -636,17 +679,22 @@ class DADF5():
requested = [{'label':x,'arg':'x'}] requested = [{'label':x,'arg':'x'}]
self.__add_generic_pointwise(principal_components,requested) self.__add_generic_pointwise(__add_principal_components,requested)
def add_maximum_shear(self,x): def add_maximum_shear(self,x):
"""Adds maximum shear components of symmetric tensor.""" """
def maximum_shear(x): Add maximum shear components of symmetric tensor.
w = np.linalg.eigvalsh((x['data']+np.transpose(x['data'],(0,2,1)))*.5) # eigenvalues (of symmetric(ed) matrix) Parameters
----------
x : str
Label of the dataset containing a symmetric tensor.
"""
def __add_maximum_shear(x):
return { return {
'data' : (w[:,2] - w[:,0])*0.5, 'data': mechanics.maximum_shear(x),
'label': 'max_shear({})'.format(x['label']), 'label': 'max_shear({})'.format(x['label']),
'meta': { 'meta': {
'Unit': x['meta']['Unit'], 'Unit': x['meta']['Unit'],
@ -657,7 +705,7 @@ class DADF5():
requested = [{'label':x,'arg':'x'}] requested = [{'label':x,'arg':'x'}]
self.__add_generic_pointwise(maximum_shear,requested) self.__add_generic_pointwise(__add_maximum_shear,requested)
def __add_generic_pointwise(self,func,datasets_requested,extra_args={}): def __add_generic_pointwise(self,func,datasets_requested,extra_args={}):
@ -672,7 +720,6 @@ class DADF5():
Details of the datasets to be used: label (in HDF5 file) and arg (argument to which the data is parsed in func). Details of the datasets to be used: label (in HDF5 file) and arg (argument to which the data is parsed in func).
extra_args : dictionary, optional extra_args : dictionary, optional
Any extra arguments parsed to func. Any extra arguments parsed to func.
""" """
def job(args): def job(args):
"""Call function with input data + extra arguments, returns results + group.""" """Call function with input data + extra arguments, returns results + group."""

134
python/damask/mechanics.py
View File

@ -2,7 +2,7 @@ import numpy as np
def Cauchy(F,P): def Cauchy(F,P):
""" """
Calculate Cauchy stress from 1st Piola-Kirchhoff stress and deformation gradient. Return Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric.
@ -11,51 +11,81 @@ def Cauchy(F,P):
F : numpy.array of shape (x,3,3) or (3,3) F : numpy.array of shape (x,3,3) or (3,3)
Deformation gradient. Deformation gradient.
P : numpy.array of shape (x,3,3) or (3,3) P : numpy.array of shape (x,3,3) or (3,3)
1st Piola-Kirchhoff. 1. Piola-Kirchhoff stress.
""" """
if np.shape(F) == np.shape(P) == (3,3): if np.shape(F) == np.shape(P) == (3,3):
sigma = 1.0/np.linalg.det(F) * np.dot(F,P) sigma = 1.0/np.linalg.det(F) * np.dot(F,P)
return (sigma+sigma.T)*0.5
else: else:
sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F) sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F)
return (sigma + np.transpose(sigma,(0,2,1)))*0.5 return symmetric(sigma)
def strain_tensor(F,t,ord):
"""
Return strain tensor calculated from deformation gradient.
For details refer to Albrecht Bertram: Elasticity and Plasticity of Large Deformations:
An Introduction (3rd Edition, 2012), p. 102.
Parameters
----------
F : numpy.array of shape (x,3,3) or (3,3)
Deformation gradient.
t : {‘V’, ‘U’}
Type of the polar decomposition, ‘V’ for right stretch tensor and ‘U’ for left stretch tensor.
ord : float
Order of the strain
"""
F_expanded = F if len(F.shape) == 3 else F.reshape(1,3,3)
if t == 'U':
B = np.matmul(F_expanded,transpose(F_expanded))
U,n = np.linalg.eigh(symmetric(B))
l = np.log(U) if ord == 0 else U**ord - np.broadcast_to(np.ones(3),[U.shape[0],3])
elif t == 'V':
C = np.matmul(transpose(F_expanded),F_expanded)
V,n = np.linalg.eigh(symmetric(C))
l = np.log(V) if ord == 0 else np.broadcast_to(np.ones(3),[V.shape[0],3]) - 1.0/V**ord
epsilon = np.matmul(n,np.einsum('ij,ikj->ijk',l,n))
return epsilon.reshape((3,3)) if np.shape(F) == (3,3) else \
epsilon
def deviatoric_part(x): def deviatoric_part(x):
""" """
Calculate deviatoric part of a tensor. Return deviatoric part of a tensor.
Parameters Parameters
---------- ----------
x : numpy.array of shape (x,3,3) or (3,3) x : numpy.array of shape (x,3,3) or (3,3)
Tensor. Tensor.
""" """
if np.shape(x) == (3,3): return x - np.eye(3)*spherical_part(x) if np.shape(x) == (3,3) else \
return x - np.eye(3)*np.trace(x)/3.0 x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x))
else:
return x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),np.trace(x,axis1=1,axis2=2)/3.0)
def spherical_part(x): def spherical_part(x):
""" """
Calculate spherical(hydrostatic) part of a tensor. Return spherical (hydrostatic) part of a tensor.
A single scalar is returned, i.e. the hydrostatic part is not mapped on the 3rd order identity matrix. A single scalar is returned, i.e. the hydrostatic part is not mapped on the 3rd order identity
matrix.
Parameters Parameters
---------- ----------
x : numpy.array of shape (x,3,3) or (3,3) x : numpy.array of shape (x,3,3) or (3,3)
Tensor. Tensor.
""" """
if np.shape(x) == (3,3):
return np.trace(x)/3.0 return np.trace(x)/3.0 if np.shape(x) == (3,3) else \
else: np.trace(x,axis1=1,axis2=2)/3.0
return np.trace(x,axis1=1,axis2=2)/3.0
def Mises_stress(sigma): def Mises_stress(sigma):
""" """
Calculate the Mises equivalent of a stress tensor. Return the Mises equivalent of a stress tensor.
Parameters Parameters
---------- ----------
@ -63,23 +93,75 @@ def Mises_stress(sigma):
Symmetric stress tensor. Symmetric stress tensor.
""" """
s = deviatoric_part(sigma) s = deviatoric_part(sigma)
if np.shape(sigma) == (3,3): return np.sqrt(3.0/2.0*np.trace(s)) if np.shape(sigma) == (3,3) else \
return np.sqrt(3.0/2.0*np.trace(s)) np.sqrt(3.0/2.0*np.einsum('ijk->i',s))
else:
return np.sqrt(np.einsum('ijk->i',s)*3.0/2.0)
def Mises_strain(epsilon): def Mises_strain(epsilon):
""" """
Calculate the Mises equivalent of a strain tensor. Return the Mises equivalent of a strain tensor.
Parameters Parameters
---------- ----------
sigma : numpy.array of shape (x,3,3) or (3,3) epsilon : numpy.array of shape (x,3,3) or (3,3)
Symmetric strain tensor. Symmetric strain tensor.
""" """
s = deviatoric_part(epsilon) s = deviatoric_part(epsilon)
if np.shape(epsilon) == (3,3): return np.sqrt(2.0/3.0*np.trace(s)) if np.shape(epsilon) == (3,3) else \
return np.sqrt(2.0/3.0*np.trace(s)) np.sqrt(2.0/3.0*np.einsum('ijk->i',s))
else:
return np.sqrt(2.0/3.0*np.einsum('ijk->i',s))
def symmetric(x):
"""
Return the symmetrized tensor.
Parameters
----------
x : numpy.array of shape (x,3,3) or (3,3)
Tensor.
"""
return (x+transpose(x))*0.5
def maximum_shear(x):
"""
Return the maximum shear component of a symmetric tensor.
Parameters
----------
x : numpy.array of shape (x,3,3) or (3,3)
Symmetric tensor.
"""
w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order
return (w[2] - w[0])*0.5 if np.shape(epsilon) == (3,3) else \
(w[:,2] - w[:,0])*0.5
def principal_components(x):
"""
Return the principal components of a symmetric tensor.
The principal components (eigenvalues) are sorted in descending order, each repeated according to
its multiplicity.
Parameters
----------
x : numpy.array of shape (x,3,3) or (3,3)
Symmetric tensor.
"""
w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order
return w[::-1] if np.shape(epsilon) == (3,3) else \
w[:,::-1]
def transpose(x):
"""
Return the transpose of a tensor.
Parameters
----------
x : numpy.array of shape (x,3,3) or (3,3)
Tensor.
"""
return x.T if np.shape(x) == (3,3) else \
np.transpose(x,(0,2,1))