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

Merge branch 'development' into report-systems

This commit is contained in:
Vitesh Shah 2021-08-10 09:16:13 +02:00
parent cf66b81e3a 7f67ca243c
commit b29b0ecd71
23 changed files with 1563 additions and 1577 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit 9699f20f21f8a5f532c735a1aa9daeba395da94d Subproject commit 4ce625b4ac0da9d490620f8cf1694d0a057cfa47

117
processing/legacy/addSchmidfactors.py vendored
View File

@ -1,117 +0,0 @@
#!/usr/bin/env python3
import os
import sys
from io import StringIO
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
slipSystems = {
'fcc': damask.lattice.kinematics['cF']['slip'][:12],
'bcc': damask.lattice.kinematics['cI']['slip'],
'hex': damask.lattice.kinematics['hP']['slip'],
}
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(usage='%prog options [ASCIItable(s)]', description = """
Add columns listing Schmid factors (and optional trace vector of selected system) for given Euler angles.
""", version = scriptID)
lattice_choices = list(slipSystems.keys())
parser.add_option('-l',
'--lattice',
dest = 'lattice', type = 'choice', choices = lattice_choices, metavar='string',
help = 'type of lattice structure [%default] {}'.format(lattice_choices))
parser.add_option('--covera',
dest = 'CoverA', type = 'float', metavar = 'float',
help = 'C over A ratio for hexagonal systems [%default]')
parser.add_option('-f',
'--force',
dest = 'force',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'force direction in lab frame [%default]')
parser.add_option('-n',
'--normal',
dest = 'normal',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'stress plane normal in lab frame, per default perpendicular to the force')
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.set_defaults(force = (0.0,0.0,1.0),
quaternion='orientation',
normal = None,
lattice = lattice_choices[0],
CoverA = np.sqrt(8./3.),
)
(options, filenames) = parser.parse_args()
if filenames == []: filenames = [None]
force = np.array(options.force)/np.linalg.norm(options.force)
if options.normal is not None:
normal = np.array(options.normal)/np.linalg.norm(options.ormal)
if abs(np.dot(force,normal)) > 1e-3:
parser.error('stress plane normal not orthogonal to force direction')
else:
normal = force
if options.lattice in ['bcc','fcc']:
slip_direction = slipSystems[options.lattice][:,:3]
slip_normal = slipSystems[options.lattice][:,3:]
elif options.lattice == 'hex':
slip_direction = np.zeros((len(slipSystems['hex']),3),'d')
slip_normal = np.zeros_like(slip_direction)
# convert 4 Miller index notation of hex to orthogonal 3 Miller index notation
for i in range(len(slip_direction)):
slip_direction[i] = np.array([slipSystems['hex'][i,0]*1.5,
(slipSystems['hex'][i,0] + 2.*slipSystems['hex'][i,1])*0.5*np.sqrt(3),
slipSystems['hex'][i,3]*options.CoverA,
])
slip_normal[i] = np.array([slipSystems['hex'][i,4],
(slipSystems['hex'][i,4] + 2.*slipSystems['hex'][i,5])/np.sqrt(3),
slipSystems['hex'][i,7]/options.CoverA,
])
slip_direction /= np.linalg.norm(slip_direction,axis=1,keepdims=True)
slip_normal /= np.linalg.norm(slip_normal, axis=1,keepdims=True)
labels = ['S[{direction[0]:.1g}_{direction[1]:.1g}_{direction[2]:.1g}]'
'({normal[0]:.1g}_{normal[1]:.1g}_{normal[2]:.1g})'\
.format(normal = theNormal, direction = theDirection,
) for theNormal,theDirection in zip(slip_normal,slip_direction)]
for name in filenames:
damask.util.report(scriptName,name)
table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
o = damask.Rotation.from_quaternion(table.get(options.quaternion))
force = np.broadcast_to(force, o.shape+(3,))
normal = np.broadcast_to(normal,o.shape+(3,))
slip_direction = np.broadcast_to(slip_direction,o.shape+slip_direction.shape)
slip_normal = np.broadcast_to(slip_normal, o.shape+slip_normal.shape)
S = np.abs(np.einsum('ijk,ik->ij',slip_direction,(o@force))*
np.einsum('ijk,ik->ij',slip_normal, (o@normal)))
for i,label in enumerate(labels):
table = table.add(label,S[:,i],scriptID+' '+' '.join(sys.argv[1:]))
table.save((sys.stdout if name is None else name))

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha4-221-g4a8c83611 v3.0.0-alpha4-276-gf75235f6a

2
python/damask/__init__.py
View File

@ -14,10 +14,10 @@ from . import tensor # noqa
from . import mechanics # noqa from . import mechanics # noqa
from . import solver # noqa from . import solver # noqa
from . import grid_filters # noqa from . import grid_filters # noqa
from . import lattice # noqa
#Modules that contain only one class (of the same name), are prefixed by a '_'. #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'. #For example, '_colormap' containsa class called 'Colormap' which is imported as 'damask.Colormap'.
from ._rotation import Rotation # noqa from ._rotation import Rotation # noqa
from ._crystal import Crystal # noqa
from ._orientation import Orientation # noqa from ._orientation import Orientation # noqa
from ._table import Table # noqa from ._table import Table # noqa
from ._vtk import VTK # noqa from ._vtk import VTK # noqa

849
python/damask/_crystal.py Normal file
View File

@ -0,0 +1,849 @@
import numpy as np
from . import util
from . import Rotation
lattice_symmetries = {
'aP': 'triclinic',
'mP': 'monoclinic',
'mS': 'monoclinic',
'oP': 'orthorhombic',
'oS': 'orthorhombic',
'oI': 'orthorhombic',
'oF': 'orthorhombic',
'tP': 'tetragonal',
'tI': 'tetragonal',
'hP': 'hexagonal',
'cP': 'cubic',
'cI': 'cubic',
'cF': 'cubic',
}
class Crystal():
"""Crystal lattice."""
def __init__(self,*,
family = None,
lattice = None,
a = None,b = None,c = None,
alpha = None,beta = None,gamma = None,
degrees = False):
"""
Representation of crystal in terms of crystal family or Bravais lattice.
Parameters
----------
family : {'triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic'}, optional.
Name of the crystal family.
Will be inferred if 'lattice' is given.
lattice : {'aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF'}, optional.
Name of the Bravais lattice in Pearson notation.
a : float, optional
Length of lattice parameter 'a'.
b : float, optional
Length of lattice parameter 'b'.
c : float, optional
Length of lattice parameter 'c'.
alpha : float, optional
Angle between b and c lattice basis.
beta : float, optional
Angle between c and a lattice basis.
gamma : float, optional
Angle between a and b lattice basis.
degrees : bool, optional
Angles are given in degrees. Defaults to False.
"""
if family not in [None] + list(lattice_symmetries.values()):
raise KeyError(f'invalid crystal family "{family}"')
if lattice is not None and family is not None and family != lattice_symmetries[lattice]:
raise KeyError(f'incompatible family "{family}" for lattice "{lattice}"')
self.family = lattice_symmetries[lattice] if family is None else family
self.lattice = lattice
if self.lattice is not None:
self.a = 1 if a is None else a
self.b = b
self.c = c
self.a = float(self.a) if self.a is not None else \
(self.b / self.ratio['b'] if self.b is not None and self.ratio['b'] is not None else
self.c / self.ratio['c'] if self.c is not None and self.ratio['c'] is not None else None)
self.b = float(self.b) if self.b is not None else \
(self.a * self.ratio['b'] if self.a is not None and self.ratio['b'] is not None else
self.c / self.ratio['c'] * self.ratio['b']
if self.c is not None and self.ratio['b'] is not None and self.ratio['c'] is not None else None)
self.c = float(self.c) if self.c is not None else \
(self.a * self.ratio['c'] if self.a is not None and self.ratio['c'] is not None else
self.b / self.ratio['b'] * self.ratio['c']
if self.c is not None and self.ratio['b'] is not None and self.ratio['c'] is not None else None)
self.alpha = np.radians(alpha) if degrees and alpha is not None else alpha
self.beta = np.radians(beta) if degrees and beta is not None else beta
self.gamma = np.radians(gamma) if degrees and gamma is not None else gamma
if self.alpha is None and 'alpha' in self.immutable: self.alpha = self.immutable['alpha']
if self.beta is None and 'beta' in self.immutable: self.beta = self.immutable['beta']
if self.gamma is None and 'gamma' in self.immutable: self.gamma = self.immutable['gamma']
if \
(self.a is None) \
or (self.b is None or ('b' in self.immutable and self.b != self.immutable['b'] * self.a)) \
or (self.c is None or ('c' in self.immutable and self.c != self.immutable['c'] * self.b)) \
or (self.alpha is None or ('alpha' in self.immutable and self.alpha != self.immutable['alpha'])) \
or (self.beta is None or ('beta' in self.immutable and self.beta != self.immutable['beta'])) \
or (self.gamma is None or ('gamma' in self.immutable and self.gamma != self.immutable['gamma'])):
raise ValueError (f'Incompatible parameters {self.parameters} for crystal family {self.family}')
if np.any(np.array([self.alpha,self.beta,self.gamma]) <= 0):
raise ValueError ('Lattice angles must be positive')
if np.any([np.roll([self.alpha,self.beta,self.gamma],r)[0]
> np.sum(np.roll([self.alpha,self.beta,self.gamma],r)[1:]) for r in range(3)]):
raise ValueError ('Each lattice angle must be less than sum of others')
else:
self.a = self.b = self.c = None
self.alpha = self.beta = self.gamma = None
def __repr__(self):
"""Represent."""
return '\n'.join([f'Crystal family {self.family}']
+ ([] if self.lattice is None else [f'Bravais lattice {self.lattice}']+
list(map(lambda x:f'{x[0]}:{x[1]:.5g}',
zip(['a','b','c','alpha','beta','gamma',],
self.parameters))))
)
def __eq__(self,other):
"""
Equal to other.
Parameters
----------
other : Crystal
Crystal to check for equality.
"""
return self.lattice == other.lattice and \
self.parameters == other.parameters and \
self.family == other.family
@property
def parameters(self):
"""Return lattice parameters a, b, c, alpha, beta, gamma."""
return (self.a,self.b,self.c,self.alpha,self.beta,self.gamma)
@property
def immutable(self):
"""Return immutable lattice parameters."""
_immutable = {
'cubic': {
'b': 1.0,
'c': 1.0,
'alpha': np.pi/2.,
'beta': np.pi/2.,
'gamma': np.pi/2.,
},
'hexagonal': {
'b': 1.0,
'alpha': np.pi/2.,
'beta': np.pi/2.,
'gamma': 2.*np.pi/3.,
},
'tetragonal': {
'b': 1.0,
'alpha': np.pi/2.,
'beta': np.pi/2.,
'gamma': np.pi/2.,
},
'orthorhombic': {
'alpha': np.pi/2.,
'beta': np.pi/2.,
'gamma': np.pi/2.,
},
'monoclinic': {
'alpha': np.pi/2.,
'gamma': np.pi/2.,
},
'triclinic': {}
}
return _immutable[self.family]
@property
def standard_triangle(self):
"""
Corners of the standard triangle.
Notes
-----
Not yet defined for monoclinic.
References
----------
Bases are computed from
>>> basis = {
... 'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,1.]/np.sqrt(2.), # green
... [1.,1.,1.]/np.sqrt(3.)]).T), # blue
... 'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T), # blue
... 'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [1.,1.,0.]/np.sqrt(2.)]).T), # blue
... 'orthorhombic': np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [0.,1.,0.]]).T), # blue
... }
"""
_basis = {
'cubic': {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
'proper':np.array([ [ 0. , -1. , 1. ],
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
},
'hexagonal':
{'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
},
'tetragonal':
{'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
},
'orthorhombic':
{'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
'proper':np.array([ [ 0., 0., 1.],
[-1., 0., 0.],
[ 0., 1., 0.] ]),
}}
return _basis.get(self.family,None)
@property
def ratio(self):
"""Return axes ratios of own lattice."""
_ratio = { 'hexagonal': {'c': np.sqrt(8./3.)}}
return dict(b = self.immutable['b']
if 'b' in self.immutable else
_ratio[self.family]['b'] if self.family in _ratio and 'b' in _ratio[self.family] else None,
c = self.immutable['c']
if 'c' in self.immutable else
_ratio[self.family]['c'] if self.family in _ratio and 'c' in _ratio[self.family] else None,
)
@property
def basis_real(self):
"""
Return orthogonal real space crystal basis.
References
----------
C.T. Young and J.L. Lytton, Journal of Applied Physics 43:1408–1417, 1972
https://doi.org/10.1063/1.1661333
"""
if None in self.parameters:
raise KeyError('missing crystal lattice parameters')
return np.array([
[1,0,0],
[np.cos(self.gamma),np.sin(self.gamma),0],
[np.cos(self.beta),
(np.cos(self.alpha)-np.cos(self.beta)*np.cos(self.gamma)) /np.sin(self.gamma),
np.sqrt(1 - np.cos(self.alpha)**2 - np.cos(self.beta)**2 - np.cos(self.gamma)**2
+ 2 * np.cos(self.alpha) * np.cos(self.beta) * np.cos(self.gamma))/np.sin(self.gamma)],
],dtype=float).T \
* np.array([self.a,self.b,self.c])
@property
def basis_reciprocal(self):
"""Return reciprocal (dual) crystal basis."""
return np.linalg.inv(self.basis_real.T)
def to_lattice(self,*,direction=None,plane=None):
"""
Calculate lattice vector corresponding to crystal frame direction or plane normal.
Parameters
----------
direction|plane : numpy.ndarray of shape (...,3)
Vector along direction or plane normal.
Returns
-------
Miller : numpy.ndarray of shape (...,3)
Lattice vector of direction or plane.
Use util.scale_to_coprime to convert to (integer) Miller indices.
"""
if (direction is not None) ^ (plane is None):
raise KeyError('Specify either "direction" or "plane"')
axis,basis = (np.array(direction),self.basis_reciprocal.T) \
if plane is None else \
(np.array(plane),self.basis_real.T)
return np.einsum('il,...l',basis,axis)
def to_frame(self,*,uvw=None,hkl=None):
"""
Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
Parameters
----------
uvw|hkl : numpy.ndarray of shape (...,3)
Miller indices of crystallographic direction or plane normal.
Returns
-------
vector : numpy.ndarray of shape (...,3) or (N,...,3)
Crystal frame vector (or vectors if with_symmetry) along [uvw] direction or (hkl) plane normal.
"""
if (uvw is not None) ^ (hkl is None):
raise KeyError('Specify either "uvw" or "hkl"')
axis,basis = (np.array(uvw),self.basis_real) \
if hkl is None else \
(np.array(hkl),self.basis_reciprocal)
return np.einsum('il,...l',basis,axis)
def kinematics(self,mode):
"""
Return crystal kinematics systems.
Parameters
----------
mode : {'slip','twin'}
Deformation mode.
Returns
-------
direction_plane : dictionary
Directions and planes of deformation mode families.
"""
_kinematics = {
'cF': {
'slip' :[np.array([
[+0,+1,-1, +1,+1,+1],
[-1,+0,+1, +1,+1,+1],
[+1,-1,+0, +1,+1,+1],
[+0,-1,-1, -1,-1,+1],
[+1,+0,+1, -1,-1,+1],
[-1,+1,+0, -1,-1,+1],
[+0,-1,+1, +1,-1,-1],
[-1,+0,-1, +1,-1,-1],
[+1,+1,+0, +1,-1,-1],
[+0,+1,+1, -1,+1,-1],
[+1,+0,-1, -1,+1,-1],
[-1,-1,+0, -1,+1,-1]]),
np.array([
[+1,+1,+0, +1,-1,+0],
[+1,-1,+0, +1,+1,+0],
[+1,+0,+1, +1,+0,-1],
[+1,+0,-1, +1,+0,+1],
[+0,+1,+1, +0,+1,-1],
[+0,+1,-1, +0,+1,+1]])],
'twin' :[np.array([
[-2, 1, 1, 1, 1, 1],
[ 1,-2, 1, 1, 1, 1],
[ 1, 1,-2, 1, 1, 1],
[ 2,-1, 1, -1,-1, 1],
[-1, 2, 1, -1,-1, 1],
[-1,-1,-2, -1,-1, 1],
[-2,-1,-1, 1,-1,-1],
[ 1, 2,-1, 1,-1,-1],
[ 1,-1, 2, 1,-1,-1],
[ 2, 1,-1, -1, 1,-1],
[-1,-2,-1, -1, 1,-1],
[-1, 1, 2, -1, 1,-1]])]
},
'cI': {
'slip' :[np.array([
[+1,-1,+1, +0,+1,+1],
[-1,-1,+1, +0,+1,+1],
[+1,+1,+1, +0,-1,+1],
[-1,+1,+1, +0,-1,+1],
[-1,+1,+1, +1,+0,+1],
[-1,-1,+1, +1,+0,+1],
[+1,+1,+1, -1,+0,+1],
[+1,-1,+1, -1,+0,+1],
[-1,+1,+1, +1,+1,+0],
[-1,+1,-1, +1,+1,+0],
[+1,+1,+1, -1,+1,+0],
[+1,+1,-1, -1,+1,+0]]),
np.array([
[-1,+1,+1, +2,+1,+1],
[+1,+1,+1, -2,+1,+1],
[+1,+1,-1, +2,-1,+1],
[+1,-1,+1, +2,+1,-1],
[+1,-1,+1, +1,+2,+1],
[+1,+1,-1, -1,+2,+1],
[+1,+1,+1, +1,-2,+1],
[-1,+1,+1, +1,+2,-1],
[+1,+1,-1, +1,+1,+2],
[+1,-1,+1, -1,+1,+2],
[-1,+1,+1, +1,-1,+2],
[+1,+1,+1, +1,+1,-2]]),
np.array([
[+1,+1,-1, +1,+2,+3],
[+1,-1,+1, -1,+2,+3],
[-1,+1,+1, +1,-2,+3],
[+1,+1,+1, +1,+2,-3],
[+1,-1,+1, +1,+3,+2],
[+1,+1,-1, -1,+3,+2],
[+1,+1,+1, +1,-3,+2],
[-1,+1,+1, +1,+3,-2],
[+1,+1,-1, +2,+1,+3],
[+1,-1,+1, -2,+1,+3],
[-1,+1,+1, +2,-1,+3],
[+1,+1,+1, +2,+1,-3],
[+1,-1,+1, +2,+3,+1],
[+1,+1,-1, -2,+3,+1],
[+1,+1,+1, +2,-3,+1],
[-1,+1,+1, +2,+3,-1],
[-1,+1,+1, +3,+1,+2],
[+1,+1,+1, -3,+1,+2],
[+1,+1,-1, +3,-1,+2],
[+1,-1,+1, +3,+1,-2],
[-1,+1,+1, +3,+2,+1],
[+1,+1,+1, -3,+2,+1],
[+1,+1,-1, +3,-2,+1],
[+1,-1,+1, +3,+2,-1]])],
'twin' :[np.array([
[-1, 1, 1, 2, 1, 1],
[ 1, 1, 1, -2, 1, 1],
[ 1, 1,-1, 2,-1, 1],
[ 1,-1, 1, 2, 1,-1],
[ 1,-1, 1, 1, 2, 1],
[ 1, 1,-1, -1, 2, 1],
[ 1, 1, 1, 1,-2, 1],
[-1, 1, 1, 1, 2,-1],
[ 1, 1,-1, 1, 1, 2],
[ 1,-1, 1, -1, 1, 2],
[-1, 1, 1, 1,-1, 2],
[ 1, 1, 1, 1, 1,-2]])]
},
'hP': {
'slip' :[np.array([
[+2,-1,-1,+0, +0,+0,+0,+1],
[-1,+2,-1,+0, +0,+0,+0,+1],
[-1,-1,+2,+0, +0,+0,+0,+1]]),
np.array([
[+2,-1,-1,+0, +0,+1,-1,+0],
[-1,+2,-1,+0, -1,+0,+1,+0],
[-1,-1,+2,+0, +1,-1,+0,+0]]),
np.array([
[-1,+1,+0,+0, +1,+1,-2,+0],
[+0,-1,+1,+0, -2,+1,+1,+0],
[+1,+0,-1,+0, +1,-2,+1,+0]]),
np.array([
[-1,+2,-1,+0, +1,+0,-1,+1],
[-2,+1,+1,+0, +0,+1,-1,+1],
[-1,-1,+2,+0, -1,+1,+0,+1],
[+1,-2,+1,+0, -1,+0,+1,+1],
[+2,-1,-1,+0, +0,-1,+1,+1],
[+1,+1,-2,+0, +1,-1,+0,+1]]),
np.array([
[-2,+1,+1,+3, +1,+0,-1,+1],
[-1,-1,+2,+3, +1,+0,-1,+1],
[-1,-1,+2,+3, +0,+1,-1,+1],
[+1,-2,+1,+3, +0,+1,-1,+1],
[+1,-2,+1,+3, -1,+1,+0,+1],
[+2,-1,-1,+3, -1,+1,+0,+1],
[+2,-1,-1,+3, -1,+0,+1,+1],
[+1,+1,-2,+3, -1,+0,+1,+1],
[+1,+1,-2,+3, +0,-1,+1,+1],
[-1,+2,-1,+3, +0,-1,+1,+1],
[-1,+2,-1,+3, +1,-1,+0,+1],
[-2,+1,+1,+3, +1,-1,+0,+1]]),
np.array([
[-1,-1,+2,+3, +1,+1,-2,+2],
[+1,-2,+1,+3, -1,+2,-1,+2],
[+2,-1,-1,+3, -2,+1,+1,+2],
[+1,+1,-2,+3, -1,-1,+2,+2],
[-1,+2,-1,+3, +1,-2,+1,+2],
[-2,+1,+1,+3, +2,-1,-1,+2]])],
'twin' :[np.array([
[-1, 0, 1, 1, 1, 0,-1, 2], # shear = (3-(c/a)^2)/(sqrt(3) c/a) <-10.1>{10.2}
[ 0,-1, 1, 1, 0, 1,-1, 2],
[ 1,-1, 0, 1, -1, 1, 0, 2],
[ 1, 0,-1, 1, -1, 0, 1, 2],
[ 0, 1,-1, 1, 0,-1, 1, 2],
[-1, 1, 0, 1, 1,-1, 0, 2]]),
np.array([
[-1,-1, 2, 6, 1, 1,-2, 1], # shear = 1/(c/a) <11.6>{-1-1.1}
[ 1,-2, 1, 6, -1, 2,-1, 1],
[ 2,-1,-1, 6, -2, 1, 1, 1],
[ 1, 1,-2, 6, -1,-1, 2, 1],
[-1, 2,-1, 6, 1,-2, 1, 1],
[-2, 1, 1, 6, 2,-1,-1, 1]]),
np.array([
[ 1, 0,-1,-2, 1, 0,-1, 1], # shear = (4(c/a)^2-9)/(4 sqrt(3) c/a) <10.-2>{10.1}
[ 0, 1,-1,-2, 0, 1,-1, 1],
[-1, 1, 0,-2, -1, 1, 0, 1],
[-1, 0, 1,-2, -1, 0, 1, 1],
[ 0,-1, 1,-2, 0,-1, 1, 1],
[ 1,-1, 0,-2, 1,-1, 0, 1]]),
np.array([
[ 1, 1,-2,-3, 1, 1,-2, 2], # shear = 2((c/a)^2-2)/(3 c/a) <11.-3>{11.2}
[-1, 2,-1,-3, -1, 2,-1, 2],
[-2, 1, 1,-3, -2, 1, 1, 2],
[-1,-1, 2,-3, -1,-1, 2, 2],
[ 1,-2, 1,-3, 1,-2, 1, 2],
[ 2,-1,-1,-3, 2,-1,-1, 2]])]
},
}
master = _kinematics[self.lattice][mode]
if self.lattice == 'hP':
return {'direction':[util.Bravais_to_Miller(uvtw=m[:,0:4]) for m in master],
'plane': [util.Bravais_to_Miller(hkil=m[:,4:8]) for m in master]}
else:
return {'direction':[m[:,0:3] for m in master],
'plane': [m[:,3:6] for m in master]}
def relation_operations(self,model):
"""
Crystallographic orientation relationships for phase transformations.
Parameters
----------
model : str
Name of orientation relationship.
Returns
-------
operations : (string, damask.Rotation)
Resulting lattice and rotations characterizing the orientation relationship.
References
----------
S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
https://doi.org/10.1016/j.jallcom.2012.02.004
K. Kitahara et al., Acta Materialia 54(5):1279-1288, 2006
https://doi.org/10.1016/j.actamat.2005.11.001
Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
https://doi.org/10.1107/S0021889805038276
H. Kitahara et al., Materials Characterization 54(4-5):378-386, 2005
https://doi.org/10.1016/j.matchar.2004.12.015
Y. He et al., Acta Materialia 53(4):1179-1190, 2005
https://doi.org/10.1016/j.actamat.2004.11.021
"""
_orientation_relationships = {
'KS': {
'cF' : np.array([
[[-1, 0, 1],[ 1, 1, 1]],
[[-1, 0, 1],[ 1, 1, 1]],
[[ 0, 1,-1],[ 1, 1, 1]],
[[ 0, 1,-1],[ 1, 1, 1]],
[[ 1,-1, 0],[ 1, 1, 1]],
[[ 1,-1, 0],[ 1, 1, 1]],
[[ 1, 0,-1],[ 1,-1, 1]],
[[ 1, 0,-1],[ 1,-1, 1]],
[[-1,-1, 0],[ 1,-1, 1]],
[[-1,-1, 0],[ 1,-1, 1]],
[[ 0, 1, 1],[ 1,-1, 1]],
[[ 0, 1, 1],[ 1,-1, 1]],
[[ 0,-1, 1],[-1, 1, 1]],
[[ 0,-1, 1],[-1, 1, 1]],
[[-1, 0,-1],[-1, 1, 1]],
[[-1, 0,-1],[-1, 1, 1]],
[[ 1, 1, 0],[-1, 1, 1]],
[[ 1, 1, 0],[-1, 1, 1]],
[[-1, 1, 0],[ 1, 1,-1]],
[[-1, 1, 0],[ 1, 1,-1]],
[[ 0,-1,-1],[ 1, 1,-1]],
[[ 0,-1,-1],[ 1, 1,-1]],
[[ 1, 0, 1],[ 1, 1,-1]],
[[ 1, 0, 1],[ 1, 1,-1]],
],dtype=float),
'cI' : np.array([
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
[[-1,-1, 1],[ 0, 1, 1]],
[[-1, 1,-1],[ 0, 1, 1]],
],dtype=float),
},
'GT': {
'cF' : np.array([
[[ -5,-12, 17],[ 1, 1, 1]],
[[ 17, -5,-12],[ 1, 1, 1]],
[[-12, 17, -5],[ 1, 1, 1]],
[[ 5, 12, 17],[ -1, -1, 1]],
[[-17, 5,-12],[ -1, -1, 1]],
[[ 12,-17, -5],[ -1, -1, 1]],
[[ -5, 12,-17],[ -1, 1, 1]],
[[ 17, 5, 12],[ -1, 1, 1]],
[[-12,-17, 5],[ -1, 1, 1]],
[[ 5,-12,-17],[ 1, -1, 1]],
[[-17, -5, 12],[ 1, -1, 1]],
[[ 12, 17, 5],[ 1, -1, 1]],
[[ -5, 17,-12],[ 1, 1, 1]],
[[-12, -5, 17],[ 1, 1, 1]],
[[ 17,-12, -5],[ 1, 1, 1]],
[[ 5,-17,-12],[ -1, -1, 1]],
[[ 12, 5, 17],[ -1, -1, 1]],
[[-17, 12, -5],[ -1, -1, 1]],
[[ -5,-17, 12],[ -1, 1, 1]],
[[-12, 5,-17],[ -1, 1, 1]],
[[ 17, 12, 5],[ -1, 1, 1]],
[[ 5, 17, 12],[ 1, -1, 1]],
[[ 12, -5,-17],[ 1, -1, 1]],
[[-17,-12, 5],[ 1, -1, 1]],
],dtype=float),
'cI' : np.array([
[[-17, -7, 17],[ 1, 0, 1]],
[[ 17,-17, -7],[ 1, 1, 0]],
[[ -7, 17,-17],[ 0, 1, 1]],
[[ 17, 7, 17],[ -1, 0, 1]],
[[-17, 17, -7],[ -1, -1, 0]],
[[ 7,-17,-17],[ 0, -1, 1]],
[[-17, 7,-17],[ -1, 0, 1]],
[[ 17, 17, 7],[ -1, 1, 0]],
[[ -7,-17, 17],[ 0, 1, 1]],
[[ 17, -7,-17],[ 1, 0, 1]],
[[-17,-17, 7],[ 1, -1, 0]],
[[ 7, 17, 17],[ 0, -1, 1]],
[[-17, 17, -7],[ 1, 1, 0]],
[[ -7,-17, 17],[ 0, 1, 1]],
[[ 17, -7,-17],[ 1, 0, 1]],
[[ 17,-17, -7],[ -1, -1, 0]],
[[ 7, 17, 17],[ 0, -1, 1]],
[[-17, 7,-17],[ -1, 0, 1]],
[[-17,-17, 7],[ -1, 1, 0]],
[[ -7, 17,-17],[ 0, 1, 1]],
[[ 17, 7, 17],[ -1, 0, 1]],
[[ 17, 17, 7],[ 1, -1, 0]],
[[ 7,-17,-17],[ 0, -1, 1]],
[[-17, -7, 17],[ 1, 0, 1]],
],dtype=float),
},
'GT_prime': {
'cF' : np.array([
[[ 0, 1, -1],[ 7, 17, 17]],
[[ -1, 0, 1],[ 17, 7, 17]],
[[ 1, -1, 0],[ 17, 17, 7]],
[[ 0, -1, -1],[ -7,-17, 17]],
[[ 1, 0, 1],[-17, -7, 17]],
[[ 1, -1, 0],[-17,-17, 7]],
[[ 0, 1, -1],[ 7,-17,-17]],
[[ 1, 0, 1],[ 17, -7,-17]],
[[ -1, -1, 0],[ 17,-17, -7]],
[[ 0, -1, -1],[ -7, 17,-17]],
[[ -1, 0, 1],[-17, 7,-17]],
[[ -1, -1, 0],[-17, 17, -7]],
[[ 0, -1, 1],[ 7, 17, 17]],
[[ 1, 0, -1],[ 17, 7, 17]],
[[ -1, 1, 0],[ 17, 17, 7]],
[[ 0, 1, 1],[ -7,-17, 17]],
[[ -1, 0, -1],[-17, -7, 17]],
[[ -1, 1, 0],[-17,-17, 7]],
[[ 0, -1, 1],[ 7,-17,-17]],
[[ -1, 0, -1],[ 17, -7,-17]],
[[ 1, 1, 0],[ 17,-17, -7]],
[[ 0, 1, 1],[ -7, 17,-17]],
[[ 1, 0, -1],[-17, 7,-17]],
[[ 1, 1, 0],[-17, 17, -7]],
],dtype=float),
'cI' : np.array([
[[ 1, 1, -1],[ 12, 5, 17]],
[[ -1, 1, 1],[ 17, 12, 5]],
[[ 1, -1, 1],[ 5, 17, 12]],
[[ -1, -1, -1],[-12, -5, 17]],
[[ 1, -1, 1],[-17,-12, 5]],
[[ 1, -1, -1],[ -5,-17, 12]],
[[ -1, 1, -1],[ 12, -5,-17]],
[[ 1, 1, 1],[ 17,-12, -5]],
[[ -1, -1, 1],[ 5,-17,-12]],
[[ 1, -1, -1],[-12, 5,-17]],
[[ -1, -1, 1],[-17, 12, -5]],
[[ -1, -1, -1],[ -5, 17,-12]],
[[ 1, -1, 1],[ 12, 17, 5]],
[[ 1, 1, -1],[ 5, 12, 17]],
[[ -1, 1, 1],[ 17, 5, 12]],
[[ -1, 1, 1],[-12,-17, 5]],
[[ -1, -1, -1],[ -5,-12, 17]],
[[ -1, 1, -1],[-17, -5, 12]],
[[ -1, -1, 1],[ 12,-17, -5]],
[[ -1, 1, -1],[ 5,-12,-17]],
[[ 1, 1, 1],[ 17, -5,-12]],
[[ 1, 1, 1],[-12, 17, -5]],
[[ 1, -1, -1],[ -5, 12,-17]],
[[ 1, 1, -1],[-17, 5,-12]],
],dtype=float),
},
'NW': {
'cF' : np.array([
[[ 2, -1, -1],[ 1, 1, 1]],
[[ -1, 2, -1],[ 1, 1, 1]],
[[ -1, -1, 2],[ 1, 1, 1]],
[[ -2, -1, -1],[ -1, 1, 1]],
[[ 1, 2, -1],[ -1, 1, 1]],
[[ 1, -1, 2],[ -1, 1, 1]],
[[ 2, 1, -1],[ 1, -1, 1]],
[[ -1, -2, -1],[ 1, -1, 1]],
[[ -1, 1, 2],[ 1, -1, 1]],
[[ 2, -1, 1],[ -1, -1, 1]],
[[ -1, 2, 1],[ -1, -1, 1]],
[[ -1, -1, -2],[ -1, -1, 1]],
],dtype=float),
'cI' : np.array([
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
],dtype=float),
},
'Pitsch': {
'cF' : np.array([
[[ 1, 0, 1],[ 0, 1, 0]],
[[ 1, 1, 0],[ 0, 0, 1]],
[[ 0, 1, 1],[ 1, 0, 0]],
[[ 0, 1, -1],[ 1, 0, 0]],
[[ -1, 0, 1],[ 0, 1, 0]],
[[ 1, -1, 0],[ 0, 0, 1]],
[[ 1, 0, -1],[ 0, 1, 0]],
[[ -1, 1, 0],[ 0, 0, 1]],
[[ 0, -1, 1],[ 1, 0, 0]],
[[ 0, 1, 1],[ 1, 0, 0]],
[[ 1, 0, 1],[ 0, 1, 0]],
[[ 1, 1, 0],[ 0, 0, 1]],
],dtype=float),
'cI' : np.array([
[[ 1, -1, 1],[ -1, 0, 1]],
[[ 1, 1, -1],[ 1, -1, 0]],
[[ -1, 1, 1],[ 0, 1, -1]],
[[ -1, 1, -1],[ 0, -1, -1]],
[[ -1, -1, 1],[ -1, 0, -1]],
[[ 1, -1, -1],[ -1, -1, 0]],
[[ 1, -1, -1],[ -1, 0, -1]],
[[ -1, 1, -1],[ -1, -1, 0]],
[[ -1, -1, 1],[ 0, -1, -1]],
[[ -1, 1, 1],[ 0, -1, 1]],
[[ 1, -1, 1],[ 1, 0, -1]],
[[ 1, 1, -1],[ -1, 1, 0]],
],dtype=float),
},
'Bain': {
'cF' : np.array([
[[ 0, 1, 0],[ 1, 0, 0]],
[[ 0, 0, 1],[ 0, 1, 0]],
[[ 1, 0, 0],[ 0, 0, 1]],
],dtype=float),
'cI' : np.array([
[[ 0, 1, 1],[ 1, 0, 0]],
[[ 1, 0, 1],[ 0, 1, 0]],
[[ 1, 1, 0],[ 0, 0, 1]],
],dtype=float),
},
'Burgers' : {
'cI' : np.array([
[[ -1, 1, 1],[ 1, 1, 0]],
[[ -1, 1, -1],[ 1, 1, 0]],
[[ 1, 1, 1],[ 1, -1, 0]],
[[ 1, 1, -1],[ 1, -1, 0]],
[[ 1, 1, -1],[ 1, 0, 1]],
[[ -1, 1, 1],[ 1, 0, 1]],
[[ 1, 1, 1],[ -1, 0, 1]],
[[ 1, -1, 1],[ -1, 0, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, -1, 1]],
[[ 1, 1, 1],[ 0, -1, 1]],
],dtype=float),
'hP' : np.array([
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
],dtype=float),
},
}
orientation_relationships = {k:v for k,v in _orientation_relationships.items() if self.lattice in v}
if model not in orientation_relationships:
raise KeyError(f'unknown orientation relationship "{model}"')
r = orientation_relationships[model]
sl = self.lattice
ol = (set(r)-{sl}).pop()
m = r[sl]
o = r[ol]
p_,_p = np.zeros(m.shape[:-1]+(3,)),np.zeros(o.shape[:-1]+(3,))
p_[...,0,:] = m[...,0,:] if m.shape[-1] == 3 else util.Bravais_to_Miller(uvtw=m[...,0,0:4])
p_[...,1,:] = m[...,1,:] if m.shape[-1] == 3 else util.Bravais_to_Miller(hkil=m[...,1,0:4])
_p[...,0,:] = o[...,0,:] if o.shape[-1] == 3 else util.Bravais_to_Miller(uvtw=o[...,0,0:4])
_p[...,1,:] = o[...,1,:] if o.shape[-1] == 3 else util.Bravais_to_Miller(hkil=o[...,1,0:4])
return (ol,Rotation.from_parallel(p_,_p))

824
python/damask/_orientation.py
View File

@ -1,11 +1,12 @@
import inspect import inspect
import copy
import numpy as np import numpy as np
from . import Rotation from . import Rotation
from . import Crystal
from . import util from . import util
from . import tensor from . import tensor
from . import lattice as lattice_
lattice_symmetries = { lattice_symmetries = {
@ -30,10 +31,12 @@ lattice_symmetries = {
} }
_parameter_doc = \ _parameter_doc = \
"""lattice : str """
Either a crystal family out of {triclinic, monoclinic, orthorhombic, tetragonal, hexagonal, cubic} family : {'triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic'}, optional.
or a Bravais lattice out of {aP, mP, mS, oP, oS, oI, oF, tP, tI, hP, cP, cI, cF}. Name of the crystal family.
When specifying a Bravais lattice, additional lattice parameters might be required. Will be infered if 'lattice' is given.
lattice : {'aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF'}, optional.
Name of the Bravais lattice in Pearson notation.
a : float, optional a : float, optional
Length of lattice parameter 'a'. Length of lattice parameter 'a'.
b : float, optional b : float, optional
@ -52,7 +55,7 @@ _parameter_doc = \
""" """
class Orientation(Rotation): class Orientation(Rotation,Crystal):
""" """
Representation of crystallographic orientation as combination of rotation and either crystal family or Bravais lattice. Representation of crystallographic orientation as combination of rotation and either crystal family or Bravais lattice.
@ -73,7 +76,7 @@ class Orientation(Rotation):
- triclinic - triclinic
- aP : primitive - aP : primitive
- monoclininic - monoclinic
- mP : primitive - mP : primitive
- mS : base-centered - mS : base-centered
@ -104,13 +107,15 @@ class Orientation(Rotation):
-------- --------
An array of 3 x 5 random orientations reduced to the fundamental zone of tetragonal symmetry: An array of 3 x 5 random orientations reduced to the fundamental zone of tetragonal symmetry:
>>> damask.Orientation.from_random(shape=(3,5),lattice='tetragonal').reduced >>> import damask
>>> o=damask.Orientation.from_random(shape=(3,5),family='tetragonal').reduced
""" """
@util.extend_docstring(_parameter_doc) @util.extend_docstring(_parameter_doc)
def __init__(self, def __init__(self,
rotation = None, rotation = np.array([1.0,0.0,0.0,0.0]), *,
family = None,
lattice = None, lattice = None,
a = None,b = None,c = None, a = None,b = None,c = None,
alpha = None,beta = None,gamma = None, alpha = None,beta = None,gamma = None,
@ -126,79 +131,23 @@ class Orientation(Rotation):
Defaults to no rotation. Defaults to no rotation.
""" """
Rotation.__init__(self) if rotation is None else Rotation.__init__(self,rotation=rotation) Rotation.__init__(self,rotation)
Crystal.__init__(self,family=family, lattice=lattice,
if lattice in lattice_symmetries: a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma, degrees=degrees)
self.family = lattice_symmetries[lattice]
self.lattice = lattice
self.a = 1 if a is None else a
self.b = b
self.c = c
self.a = float(self.a) if self.a is not None else \
(self.b / self.ratio['b'] if self.b is not None and self.ratio['b'] is not None else
self.c / self.ratio['c'] if self.c is not None and self.ratio['c'] is not None else None)
self.b = float(self.b) if self.b is not None else \
(self.a * self.ratio['b'] if self.a is not None and self.ratio['b'] is not None else
self.c / self.ratio['c'] * self.ratio['b']
if self.c is not None and self.ratio['b'] is not None and self.ratio['c'] is not None else None)
self.c = float(self.c) if self.c is not None else \
(self.a * self.ratio['c'] if self.a is not None and self.ratio['c'] is not None else
self.b / self.ratio['b'] * self.ratio['c']
if self.c is not None and self.ratio['b'] is not None and self.ratio['c'] is not None else None)
self.alpha = np.radians(alpha) if degrees and alpha is not None else alpha
self.beta = np.radians(beta) if degrees and beta is not None else beta
self.gamma = np.radians(gamma) if degrees and gamma is not None else gamma
if self.alpha is None and 'alpha' in self.immutable: self.alpha = self.immutable['alpha']
if self.beta is None and 'beta' in self.immutable: self.beta = self.immutable['beta']
if self.gamma is None and 'gamma' in self.immutable: self.gamma = self.immutable['gamma']
if \
(self.a is None) \
or (self.b is None or ('b' in self.immutable and self.b != self.immutable['b'] * self.a)) \
or (self.c is None or ('c' in self.immutable and self.c != self.immutable['c'] * self.b)) \
or (self.alpha is None or ('alpha' in self.immutable and self.alpha != self.immutable['alpha'])) \
or (self.beta is None or ( 'beta' in self.immutable and self.beta != self.immutable['beta'])) \
or (self.gamma is None or ('gamma' in self.immutable and self.gamma != self.immutable['gamma'])):
raise ValueError (f'Incompatible parameters {self.parameters} for crystal family {self.family}')
if np.any(np.array([self.alpha,self.beta,self.gamma]) <= 0):
raise ValueError ('Lattice angles must be positive')
if np.any([np.roll([self.alpha,self.beta,self.gamma],r)[0]
> np.sum(np.roll([self.alpha,self.beta,self.gamma],r)[1:]) for r in range(3)]):
raise ValueError ('Each lattice angle must be less than sum of others')
elif lattice in set(lattice_symmetries.values()):
self.family = lattice
self.lattice = None
self.a = self.b = self.c = None
self.alpha = self.beta = self.gamma = None
else:
raise KeyError(f'Lattice "{lattice}" is unknown')
def __repr__(self): def __repr__(self):
"""Represent.""" """Represent."""
return '\n'.join(([] if self.lattice is None else [f'Bravais lattice {self.lattice}']) return '\n'.join([Crystal.__repr__(self),
+ ([f'Crystal family {self.family}']) Rotation.__repr__(self)])
+ [super().__repr__()])
def __copy__(self,**kwargs): def __copy__(self,rotation=None):
"""Create deep copy.""" """Create deep copy."""
return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion, dup = copy.deepcopy(self)
lattice =kwargs['lattice'] if 'lattice' in kwargs else self.lattice if rotation is not None:
if self.lattice is not None else self.family, dup.quaternion = Orientation(rotation,family='cubic').quaternion
a =kwargs['a'] if 'a' in kwargs else self.a, return dup
b =kwargs['b'] if 'b' in kwargs else self.b,
c =kwargs['c'] if 'c' in kwargs else self.c,
alpha =kwargs['alpha'] if 'alpha' in kwargs else self.alpha,
beta =kwargs['beta'] if 'beta' in kwargs else self.beta,
gamma =kwargs['gamma'] if 'gamma' in kwargs else self.gamma,
degrees =kwargs['degrees'] if 'degrees' in kwargs else None,
)
copy = __copy__ copy = __copy__
@ -213,8 +162,9 @@ class Orientation(Rotation):
Orientation to check for equality. Orientation to check for equality.
""" """
matching_type = all([hasattr(other,attr) and getattr(self,attr) == getattr(other,attr) matching_type = self.family == other.family and \
for attr in ['family','lattice','parameters']]) self.lattice == other.lattice and \
self.parameters == other.parameters
return np.logical_and(matching_type,super(self.__class__,self.reduced).__eq__(other.reduced)) return np.logical_and(matching_type,super(self.__class__,self.reduced).__eq__(other.reduced))
def __ne__(self,other): def __ne__(self,other):
@ -251,8 +201,9 @@ class Orientation(Rotation):
Mask indicating where corresponding orientations are close. Mask indicating where corresponding orientations are close.
""" """
matching_type = all([hasattr(other,attr) and getattr(self,attr) == getattr(other,attr) matching_type = self.family == other.family and \
for attr in ['family','lattice','parameters']]) self.lattice == other.lattice and \
self.parameters == other.parameters
return np.logical_and(matching_type,super(self.__class__,self.reduced).isclose(other.reduced)) return np.logical_and(matching_type,super(self.__class__,self.reduced).isclose(other.reduced))
@ -431,84 +382,6 @@ class Orientation(Rotation):
return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True)))) return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
@property
def symmetry_operations(self):
"""Symmetry operations as Rotations."""
if self.family == 'cubic':
sym_quats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ],
[ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, -0.5 ],
[-0.5, 0.5, -0.5, 0.5 ],
[-0.5, -0.5, 0.5, 0.5 ],
[-0.5, -0.5, 0.5, -0.5 ],
[-0.5, -0.5, -0.5, 0.5 ],
[-0.5, 0.5, -0.5, -0.5 ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ],
[-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ],
]
elif self.family == 'hexagonal':
sym_quats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ],
[ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[ 0.0, 0.0, 0.0, 1.0 ],
[-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ],
[ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
]
elif self.family == 'tetragonal':
sym_quats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
]
elif self.family == 'orthorhombic':
sym_quats = [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
]
elif self.family == 'monoclinic':
sym_quats = [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
]
elif self.family == 'triclinic':
sym_quats = [
[ 1.0,0.0,0.0,0.0 ],
]
else:
raise KeyError(f'unknown crystal family "{self.family}"')
return Rotation.from_quaternion(sym_quats,accept_homomorph=True)
@property @property
def equivalent(self): def equivalent(self):
""" """
@ -518,7 +391,8 @@ class Orientation(Rotation):
is added to the left of the Rotation array. is added to the left of the Rotation array.
""" """
o = self.symmetry_operations.broadcast_to(self.symmetry_operations.shape+self.shape,mode='right') sym_ops = self.symmetry_operations
o = sym_ops.broadcast_to(sym_ops.shape+self.shape,mode='right')
return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left')) return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
@ -619,381 +493,6 @@ class Orientation(Rotation):
return np.ones_like(rho[...,0],dtype=bool) return np.ones_like(rho[...,0],dtype=bool)
def relation_operations(self,model,return_lattice=False):
"""
Crystallographic orientation relationships for phase transformations.
Parameters
----------
model : str
Name of orientation relationship.
return_lattice : bool, optional
Return the target lattice in addition.
Returns
-------
operations : Rotations
Rotations characterizing the orientation relationship.
References
----------
S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
https://doi.org/10.1016/j.jallcom.2012.02.004
K. Kitahara et al., Acta Materialia 54(5):1279-1288, 2006
https://doi.org/10.1016/j.actamat.2005.11.001
Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
https://doi.org/10.1107/S0021889805038276
H. Kitahara et al., Materials Characterization 54(4-5):378-386, 2005
https://doi.org/10.1016/j.matchar.2004.12.015
Y. He et al., Acta Materialia 53(4):1179-1190, 2005
https://doi.org/10.1016/j.actamat.2004.11.021
"""
if model not in lattice_.relations:
raise KeyError(f'unknown orientation relationship "{model}"')
r = lattice_.relations[model]
if self.lattice not in r:
raise KeyError(f'relationship "{model}" not supported for lattice "{self.lattice}"')
sl = self.lattice
ol = (set(r)-{sl}).pop()
m = r[sl]
o = r[ol]
p_,_p = np.zeros(m.shape[:-1]+(3,)),np.zeros(o.shape[:-1]+(3,))
p_[...,0,:] = m[...,0,:] if m.shape[-1] == 3 else lattice_.Bravais_to_Miller(uvtw=m[...,0,0:4])
p_[...,1,:] = m[...,1,:] if m.shape[-1] == 3 else lattice_.Bravais_to_Miller(hkil=m[...,1,0:4])
_p[...,0,:] = o[...,0,:] if o.shape[-1] == 3 else lattice_.Bravais_to_Miller(uvtw=o[...,0,0:4])
_p[...,1,:] = o[...,1,:] if o.shape[-1] == 3 else lattice_.Bravais_to_Miller(hkil=o[...,1,0:4])
return (Rotation.from_parallel(p_,_p),ol) \
if return_lattice else \
Rotation.from_parallel(p_,_p)
def related(self,model):
"""
Orientations derived from the given relationship.
One dimension (length according to number of related orientations)
is added to the left of the Rotation array.
"""
o,lattice = self.relation_operations(model,return_lattice=True)
target = Orientation(lattice=lattice)
o = o.broadcast_to(o.shape+self.shape,mode='right')
return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'),
lattice=lattice,
b = self.b if target.ratio['b'] is None else self.a*target.ratio['b'],
c = self.c if target.ratio['c'] is None else self.a*target.ratio['c'],
alpha = None if 'alpha' in target.immutable else self.alpha,
beta = None if 'beta' in target.immutable else self.beta,
gamma = None if 'gamma' in target.immutable else self.gamma,
)
@property
def parameters(self):
"""Return lattice parameters a, b, c, alpha, beta, gamma."""
return (self.a,self.b,self.c,self.alpha,self.beta,self.gamma)
@property
def immutable(self):
"""Return immutable parameters of own lattice."""
if self.family == 'triclinic':
return {}
if self.family == 'monoclinic':
return {
'alpha': np.pi/2.,
'gamma': np.pi/2.,
}
if self.family == 'orthorhombic':
return {
'alpha': np.pi/2.,
'beta': np.pi/2.,
'gamma': np.pi/2.,
}
if self.family == 'tetragonal':
return {
'b': 1.0,
'alpha': np.pi/2.,
'beta': np.pi/2.,
'gamma': np.pi/2.,
}
if self.family == 'hexagonal':
return {
'b': 1.0,
'alpha': np.pi/2.,
'beta': np.pi/2.,
'gamma': 2.*np.pi/3.,
}
if self.family == 'cubic':
return {
'b': 1.0,
'c': 1.0,
'alpha': np.pi/2.,
'beta': np.pi/2.,
'gamma': np.pi/2.,
}
@property
def ratio(self):
"""Return axes ratios of own lattice."""
_ratio = { 'hexagonal': {'c': np.sqrt(8./3.)}}
return dict(b = self.immutable['b']
if 'b' in self.immutable else
_ratio[self.family]['b'] if self.family in _ratio and 'b' in _ratio[self.family] else None,
c = self.immutable['c']
if 'c' in self.immutable else
_ratio[self.family]['c'] if self.family in _ratio and 'c' in _ratio[self.family] else None,
)
@property
def basis_real(self):
"""
Calculate orthogonal real space crystal basis.
References
----------
C.T. Young and J.L. Lytton, Journal of Applied Physics 43:1408–1417, 1972
https://doi.org/10.1063/1.1661333
"""
if None in self.parameters:
raise KeyError('missing crystal lattice parameters')
return np.array([
[1,0,0],
[np.cos(self.gamma),np.sin(self.gamma),0],
[np.cos(self.beta),
(np.cos(self.alpha)-np.cos(self.beta)*np.cos(self.gamma)) /np.sin(self.gamma),
np.sqrt(1 - np.cos(self.alpha)**2 - np.cos(self.beta)**2 - np.cos(self.gamma)**2
+ 2 * np.cos(self.alpha) * np.cos(self.beta) * np.cos(self.gamma))/np.sin(self.gamma)],
],dtype=float).T \
* np.array([self.a,self.b,self.c])
@property
def basis_reciprocal(self):
"""Calculate reciprocal (dual) crystal basis."""
return np.linalg.inv(self.basis_real.T)
def in_SST(self,vector,proper=False):
"""
Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry.
Parameters
----------
vector : numpy.ndarray of shape (...,3)
Vector to check.
proper : bool, optional
Consider only vectors with z >= 0, hence combine two neighboring SSTs.
Defaults to False.
Returns
-------
in : numpy.ndarray of shape (...)
Boolean array indicating whether vector falls into SST.
References
----------
Bases are computed from
>>> basis = {
... 'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,1.]/np.sqrt(2.), # green
... [1.,1.,1.]/np.sqrt(3.)]).T), # blue
... 'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T), # blue
... 'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [1.,1.,0.]/np.sqrt(2.)]).T), # blue
... 'orthorhombic': np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [0.,1.,0.]]).T), # blue
... }
"""
if not isinstance(vector,np.ndarray) or vector.shape[-1] != 3:
raise ValueError('input is not a field of three-dimensional vectors')
if self.family == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
'proper':np.array([ [ 0. , -1. , 1. ],
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
}
elif self.family == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
}
elif self.family == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
}
elif self.family == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
'proper':np.array([ [ 0., 0., 1.],
[-1., 0., 0.],
[ 0., 1., 0.] ]),
}
else: # direct exit for unspecified symmetry
return np.ones_like(vector[...,0],bool)
if proper:
components_proper = np.around(np.einsum('...ji,...i',
np.broadcast_to(basis['proper'], vector.shape+(3,)),
vector), 12)
components_improper = np.around(np.einsum('...ji,...i',
np.broadcast_to(basis['improper'], vector.shape+(3,)),
vector), 12)
return np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 0.0,axis=-1)
else:
components = np.around(np.einsum('...ji,...i',
np.broadcast_to(basis['improper'], vector.shape+(3,)),
np.block([vector[...,:2],np.abs(vector[...,2:3])])), 12)
return np.all(components >= 0.0,axis=-1)
def IPF_color(self,vector,in_SST=True,proper=False):
"""
Map vector to RGB color within standard stereographic triangle of own symmetry.
Parameters
----------
vector : numpy.ndarray of shape (...,3)
Vector to colorize.
in_SST : bool, optional
Consider symmetrically equivalent orientations such that poles are located in SST.
Defaults to True.
proper : bool, optional
Consider only vectors with z >= 0, hence combine two neighboring SSTs (with mirrored colors).
Defaults to False.
Returns
-------
rgb : numpy.ndarray of shape (...,3)
RGB array of IPF colors.
Examples
--------
Inverse pole figure color of the e_3 direction for a crystal in "Cube" orientation with cubic symmetry:
>>> o = damask.Orientation(lattice='cubic')
>>> o.IPF_color([0,0,1])
array([1., 0., 0.])
References
----------
Bases are computed from
>>> basis = {
... 'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,1.]/np.sqrt(2.), # green
... [1.,1.,1.]/np.sqrt(3.)]).T), # blue
... 'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T), # blue
... 'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [1.,1.,0.]/np.sqrt(2.)]).T), # blue
... 'orthorhombic': np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # green
... [0.,1.,0.]]).T), # blue
... }
"""
if np.array(vector).shape[-1] != 3:
raise ValueError('input is not a field of three-dimensional vectors')
vector_ = self.to_SST(vector,proper) if in_SST else \
self @ np.broadcast_to(vector,self.shape+(3,))
if self.family == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
'proper':np.array([ [ 0. , -1. , 1. ],
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
}
elif self.family == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
}
elif self.family == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
}
elif self.family == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
'proper':np.array([ [ 0., 0., 1.],
[-1., 0., 0.],
[ 0., 1., 0.] ]),
}
else: # direct exit for unspecified symmetry
return np.zeros_like(vector_)
if proper:
components_proper = np.around(np.einsum('...ji,...i',
np.broadcast_to(basis['proper'], vector_.shape+(3,)),
vector_), 12)
components_improper = np.around(np.einsum('...ji,...i',
np.broadcast_to(basis['improper'], vector_.shape+(3,)),
vector_), 12)
in_SST = np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 0.0,axis=-1)
components = np.where((in_SST & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis],
components_proper,components_improper)
else:
components = np.around(np.einsum('...ji,...i',
np.broadcast_to(basis['improper'], vector_.shape+(3,)),
np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
in_SST = np.all(components >= 0.0,axis=-1)
with np.errstate(invalid='ignore',divide='ignore'):
rgb = (components/np.linalg.norm(components,axis=-1,keepdims=True))**0.5 # smoothen color ramps
rgb = np.clip(rgb,0.,1.) # clip intensity
rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1
rgb[np.broadcast_to(~in_SST[...,np.newaxis],rgb.shape)] = 0.0
return rgb
def disorientation(self,other,return_operators=False): def disorientation(self,other,return_operators=False):
""" """
Calculate disorientation between myself and given other orientation. Calculate disorientation between myself and given other orientation.
@ -1148,57 +647,185 @@ class Orientation(Rotation):
) )
def to_lattice(self,*,direction=None,plane=None): def in_SST(self,vector,proper=False):
""" """
Calculate lattice vector corresponding to crystal frame direction or plane normal. Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry.
Parameters Parameters
---------- ----------
direction|normal : numpy.ndarray of shape (...,3) vector : numpy.ndarray of shape (...,3)
Vector along direction or plane normal. Vector to check.
proper : bool, optional
Consider only vectors with z >= 0, hence combine two neighboring SSTs.
Defaults to False.
Returns Returns
------- -------
Miller : numpy.ndarray of shape (...,3) in : numpy.ndarray of shape (...)
lattice vector of direction or plane. Boolean array indicating whether vector falls into SST.
Use util.scale_to_coprime to convert to (integer) Miller indices.
""" """
if (direction is not None) ^ (plane is None): if not isinstance(vector,np.ndarray) or vector.shape[-1] != 3:
raise KeyError('specify either "direction" or "plane"') raise ValueError('input is not a field of three-dimensional vectors')
axis,basis = (np.array(direction),self.basis_reciprocal.T) \
if plane is None else \ if self.standard_triangle is None: # direct exit for no symmetry
(np.array(plane),self.basis_real.T) return np.ones_like(vector[...,0],bool)
return np.einsum('il,...l',basis,axis)
if proper:
components_proper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['proper'], vector.shape+(3,)),
vector), 12)
components_improper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)),
vector), 12)
return np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 0.0,axis=-1)
else:
components = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)),
np.block([vector[...,:2],np.abs(vector[...,2:3])])), 12)
return np.all(components >= 0.0,axis=-1)
def to_frame(self,*,uvw=None,hkl=None,with_symmetry=False): def IPF_color(self,vector,in_SST=True,proper=False):
""" """
Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl). Map vector to RGB color within standard stereographic triangle of own symmetry.
Parameters Parameters
---------- ----------
uvw|hkl : numpy.ndarray of shape (...,3) vector : numpy.ndarray of shape (...,3)
Miller indices of crystallographic direction or plane normal. Vector to colorize.
with_symmetry : bool, optional in_SST : bool, optional
Calculate all N symmetrically equivalent vectors. Consider symmetrically equivalent orientations such that poles are located in SST.
Defaults to True.
proper : bool, optional
Consider only vectors with z >= 0, hence combine two neighboring SSTs (with mirrored colors).
Defaults to False.
Returns Returns
------- -------
vector : numpy.ndarray of shape (...,3) or (N,...,3) rgb : numpy.ndarray of shape (...,3)
Crystal frame vector (or vectors if with_symmetry) along [uvw] direction or (hkl) plane normal. RGB array of IPF colors.
Examples
--------
Inverse pole figure color of the e_3 direction for a crystal in "Cube" orientation with cubic symmetry:
>>> import damask
>>> o = damask.Orientation(lattice='cubic')
>>> o.IPF_color([0,0,1])
array([1., 0., 0.])
""" """
if (uvw is not None) ^ (hkl is None): if np.array(vector).shape[-1] != 3:
raise KeyError('specify either "uvw" or "hkl"') raise ValueError('input is not a field of three-dimensional vectors')
axis,basis = (np.array(uvw),self.basis_real) \
if hkl is None else \
(np.array(hkl),self.basis_reciprocal)
return (self.symmetry_operations.broadcast_to(self.symmetry_operations.shape+axis.shape[:-1],mode='right')
@ np.broadcast_to(np.einsum('il,...l',basis,axis),self.symmetry_operations.shape+axis.shape)
if with_symmetry else
np.einsum('il,...l',basis,axis))
vector_ = self.to_SST(vector,proper) if in_SST else \
self @ np.broadcast_to(vector,self.shape+(3,))
if self.standard_triangle is None: # direct exit for no symmetry
return np.zeros_like(vector_)
if proper:
components_proper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)),
vector_), 12)
components_improper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
vector_), 12)
in_SST = np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 0.0,axis=-1)
components = np.where((in_SST & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis],
components_proper,components_improper)
else:
components = np.around(np.einsum('...ji,...i',
np.broadcast_to(self .standard_triangle['improper'], vector_.shape+(3,)),
np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
in_SST = np.all(components >= 0.0,axis=-1)
with np.errstate(invalid='ignore',divide='ignore'):
rgb = (components/np.linalg.norm(components,axis=-1,keepdims=True))**0.5 # smoothen color ramps
rgb = np.clip(rgb,0.,1.) # clip intensity
rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1
rgb[np.broadcast_to(~in_SST[...,np.newaxis],rgb.shape)] = 0.0
return rgb
@property
def symmetry_operations(self):
"""Symmetry operations as Rotations."""
_symmetry_operations = {
'cubic': [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ],
[ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, -0.5 ],
[-0.5, 0.5, -0.5, 0.5 ],
[-0.5, -0.5, 0.5, 0.5 ],
[-0.5, -0.5, 0.5, -0.5 ],
[-0.5, -0.5, -0.5, 0.5 ],
[-0.5, 0.5, -0.5, -0.5 ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ],
[-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ],
],
'hexagonal': [
[ 1.0, 0.0, 0.0, 0.0 ],
[-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ],
[ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[ 0.0, 0.0, 0.0, 1.0 ],
[-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ],
[ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
],
'tetragonal': [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
],
'orthorhombic': [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
],
'monoclinic': [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
],
'triclinic': [
[ 1.0,0.0,0.0,0.0 ],
]}
return Rotation.from_quaternion(_symmetry_operations[self.family],accept_homomorph=True)
####################################################################################################
# functions that require lattice, not just family
def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False): def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False):
""" """
@ -1217,51 +844,78 @@ class Orientation(Rotation):
Lab frame vector (or vectors if with_symmetry) along [uvw] direction or (hkl) plane normal. Lab frame vector (or vectors if with_symmetry) along [uvw] direction or (hkl) plane normal.
""" """
v = self.to_frame(uvw=uvw,hkl=hkl,with_symmetry=with_symmetry) v = self.to_frame(uvw=uvw,hkl=hkl)
if with_symmetry:
sym_ops = self.symmetry_operations
v = sym_ops.broadcast_to(sym_ops.shape+v.shape[:-1],mode='right') \
@ np.broadcast_to(v,sym_ops.shape+v.shape)
return ~(self if self.shape+v.shape[:-1] == () else self.broadcast_to(self.shape+v.shape[:-1],mode='right')) \ return ~(self if self.shape+v.shape[:-1] == () else self.broadcast_to(self.shape+v.shape[:-1],mode='right')) \
@ np.broadcast_to(v,self.shape+v.shape) @ np.broadcast_to(v,self.shape+v.shape)
def Schmid(self,mode): def Schmid(self,*,N_slip=None,N_twin=None):
u""" u"""
Calculate Schmid matrix P = d ⨂ n in the lab frame for given lattice shear kinematics. Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems.
Parameters Parameters
---------- ----------
mode : {'slip', 'twin'} N_slip|N_twin : iterable of int
Type of kinematics. Number of deformation systems per family of the deformation system.
Use '*' to select all.
Returns Returns
------- -------
P : numpy.ndarray of shape (...,N,3,3) P : numpy.ndarray of shape (N,...,3,3)
Schmid matrix for each of the N deformation systems. Schmid matrix for each of the N deformation systems.
Examples Examples
-------- --------
Schmid matrix (in lab frame) of slip systems of a face-centered Schmid matrix (in lab frame) of first octahedral slip system of a face-centered
cubic crystal in "Goss" orientation. cubic crystal in "Goss" orientation.
>>> import damask >>> import damask
>>> import numpy as np >>> import numpy as np
>>> np.set_printoptions(3,suppress=True,floatmode='fixed') >>> np.set_printoptions(3,suppress=True,floatmode='fixed')
>>> damask.Orientation.from_Eulers(phi=[0,45,0],degrees=True,lattice='cF').Schmid('slip')[0] >>> O = damask.Orientation.from_Euler_angles(phi=[0,45,0],degrees=True,lattice='cF')
>>> O.Schmid(N_slip=[1])
array([[ 0.000, 0.000, 0.000], array([[ 0.000, 0.000, 0.000],
[ 0.577, -0.000, 0.816], [ 0.577, -0.000, 0.816],
[ 0.000, 0.000, 0.000]]) [ 0.000, 0.000, 0.000]])
""" """
try: if (N_slip is not None) ^ (N_twin is None):
master = lattice_.kinematics[self.lattice][mode] raise KeyError('Specify either "N_slip" or "N_twin"')
kinematics = {'direction':master[:,0:3],'plane':master[:,3:6]} \
if master.shape[-1] == 6 else \
{'direction':lattice_.Bravais_to_Miller(uvtw=master[:,0:4]),
'plane': lattice_.Bravais_to_Miller(hkil=master[:,4:8])}
except KeyError:
raise (f'"{mode}" not defined for lattice "{self.lattice}"')
d = self.to_frame(uvw=kinematics['direction'],with_symmetry=False)
p = self.to_frame(hkl=kinematics['plane'] ,with_symmetry=False)
P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=-1,keepdims=True),
p/np.linalg.norm(p,axis=-1,keepdims=True))
return ~self.broadcast_to( self.shape+P.shape[:-2],mode='right') \ kinematics,active = (self.kinematics('slip'),N_slip) if N_twin is None else \
@ np.broadcast_to(P,self.shape+P.shape) (self.kinematics('twin'),N_twin)
if active == '*': active = [len(a) for a in kinematics['direction']]
d = self.to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)]))
p = self.to_frame(hkl=np.vstack([kinematics['plane'][i][:n] for i,n in enumerate(active)]))
P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=1,keepdims=True),
p/np.linalg.norm(p,axis=1,keepdims=True))
shape = P.shape[0:1]+self.shape+(3,3)
return ~self.broadcast_to(shape[:-2]) \
@ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape)
def related(self,model):
"""
Orientations derived from the given relationship.
One dimension (length according to number of related orientations)
is added to the left of the Rotation array.
"""
lattice,o = self.relation_operations(model)
target = Crystal(lattice=lattice)
o = o.broadcast_to(o.shape+self.shape,mode='right')
return Orientation(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'),
lattice=lattice,
b = self.b if target.ratio['b'] is None else self.a*target.ratio['b'],
c = self.c if target.ratio['c'] is None else self.a*target.ratio['c'],
alpha = None if 'alpha' in target.immutable else self.alpha,
beta = None if 'beta' in target.immutable else self.beta,
gamma = None if 'gamma' in target.immutable else self.gamma,
)

68
python/damask/_result.py
View File

@ -981,47 +981,35 @@ class Result:
self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F}) self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F})
# The add_pole functionality needs discussion.
# The new Crystal object can perform such a calculation but the outcome depends on the lattice parameters
# as well as on whether a direction or plane is concerned (see the DAMASK_examples/pole_figure notebook).
# Below code appears to be too simplistic.
# @staticmethod @staticmethod
# def _add_pole(q,p,polar): def _add_pole(q,uvw,hkl):
# pole = np.array(p) c = q['meta']['c/a'] if 'c/a' in q['meta'] else 1
# unit_pole = pole/np.linalg.norm(pole) pole = Orientation(q['data'], lattice=q['meta']['lattice'], a=1, c=c).to_pole(uvw=uvw,hkl=hkl)
# m = util.scale_to_coprime(pole)
# rot = Rotation(q['data'].view(np.double).reshape(-1,4)) return {
# 'data': pole,
# rotatedPole = rot @ np.broadcast_to(unit_pole,rot.shape+(3,)) # rotate pole according to crystal orientation 'label': 'p^[{} {} {}]'.format(*uvw) if uvw else 'p^({} {} {})'.format(*hkl),
# xy = rotatedPole[:,0:2]/(1.+abs(unit_pole[2])) # stereographic projection 'meta' : {
# coords = xy if not polar else \ 'unit': '1',
# np.block([np.sqrt(xy[:,0:1]*xy[:,0:1]+xy[:,1:2]*xy[:,1:2]),np.arctan2(xy[:,1:2],xy[:,0:1])]) 'description': f"lab frame vector along lattice {'direction' if uvw else 'plane'}",
# return { 'creator': 'add_pole'
# 'data': coords, }
# 'label': 'p^{}_[{} {} {})'.format(u'rφ' if polar else 'xy',*m), }
# 'meta' : { def add_pole(self,q='O',*,uvw=None,hkl=None):
# 'unit': '1', """
# 'description': '{} coordinates of stereographic projection of pole (direction/plane) in crystal frame'\ Add lab frame vector along lattice direction [uvw] or plane normal (hkl).
# .format('Polar' if polar else 'Cartesian'),
# 'creator': 'add_pole' Parameters
# } ----------
# } q : str
# def add_pole(self,q,p,polar=False): Name of the dataset containing the crystallographic orientation as quaternions.
# """ Defaults to 'O'.
# Add coordinates of stereographic projection of given pole in crystal frame. uvw|hkl : numpy.ndarray of shape (...,3)
# Miller indices of crystallographic direction or plane normal.
# Parameters
# ---------- """
# q : str self._add_generic_pointwise(self._add_pole,{'q':q},{'uvw':uvw,'hkl':hkl})
# Name of the dataset containing the crystallographic orientation as quaternions.
# p : numpy.array of shape (3)
# Crystallographic direction or plane.
# polar : bool, optional
# Give pole in polar coordinates. Defaults to False.
#
# """
# self._add_generic_pointwise(self._add_pole,{'q':q},{'p':p,'polar':polar})
@staticmethod @staticmethod

501
python/damask/lattice.py
View File

@ -1,501 +0,0 @@
import numpy as _np
def Bravais_to_Miller(*,uvtw=None,hkil=None):
"""
Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
Parameters
----------
uvtw|hkil : numpy.ndarray of shape (...,4)
Miller–Bravais indices of crystallographic direction [uvtw] or plane normal (hkil).
Returns
-------
uvw|hkl : numpy.ndarray of shape (...,3)
Miller indices of [uvw] direction or (hkl) plane normal.
"""
if (uvtw is not None) ^ (hkil is None):
raise KeyError('Specify either "uvtw" or "hkil"')
axis,basis = (_np.array(uvtw),_np.array([[1,0,-1,0],
[0,1,-1,0],
[0,0, 0,1]])) \
if hkil is None else \
(_np.array(hkil),_np.array([[1,0,0,0],
[0,1,0,0],
[0,0,0,1]]))
return _np.einsum('il,...l',basis,axis)
def Miller_to_Bravais(*,uvw=None,hkl=None):
"""
Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil).
Parameters
----------
uvw|hkl : numpy.ndarray of shape (...,3)
Miller indices of crystallographic direction [uvw] or plane normal (hkl).
Returns
-------
uvtw|hkil : numpy.ndarray of shape (...,4)
Miller–Bravais indices of [uvtw] direction or (hkil) plane normal.
"""
if (uvw is not None) ^ (hkl is None):
raise KeyError('Specify either "uvw" or "hkl"')
axis,basis = (_np.array(uvw),_np.array([[ 2,-1, 0],
[-1, 2, 0],
[-1,-1, 0],
[ 0, 0, 3]])/3) \
if hkl is None else \
(_np.array(hkl),_np.array([[ 1, 0, 0],
[ 0, 1, 0],
[-1,-1, 0],
[ 0, 0, 1]]))
return _np.einsum('il,...l',basis,axis)
kinematics = {
'cF': {
'slip' : _np.array([
[+0,+1,-1 , +1,+1,+1],
[-1,+0,+1 , +1,+1,+1],
[+1,-1,+0 , +1,+1,+1],
[+0,-1,-1 , -1,-1,+1],
[+1,+0,+1 , -1,-1,+1],
[-1,+1,+0 , -1,-1,+1],
[+0,-1,+1 , +1,-1,-1],
[-1,+0,-1 , +1,-1,-1],
[+1,+1,+0 , +1,-1,-1],
[+0,+1,+1 , -1,+1,-1],
[+1,+0,-1 , -1,+1,-1],
[-1,-1,+0 , -1,+1,-1],
[+1,+1,+0 , +1,-1,+0],
[+1,-1,+0 , +1,+1,+0],
[+1,+0,+1 , +1,+0,-1],
[+1,+0,-1 , +1,+0,+1],
[+0,+1,+1 , +0,+1,-1],
[+0,+1,-1 , +0,+1,+1],
],'d'),
'twin' : _np.array([
[-2, 1, 1, 1, 1, 1],
[ 1,-2, 1, 1, 1, 1],
[ 1, 1,-2, 1, 1, 1],
[ 2,-1, 1, -1,-1, 1],
[-1, 2, 1, -1,-1, 1],
[-1,-1,-2, -1,-1, 1],
[-2,-1,-1, 1,-1,-1],
[ 1, 2,-1, 1,-1,-1],
[ 1,-1, 2, 1,-1,-1],
[ 2, 1,-1, -1, 1,-1],
[-1,-2,-1, -1, 1,-1],
[-1, 1, 2, -1, 1,-1],
],dtype=float),
},
'cI': {
'slip' : _np.array([
[+1,-1,+1 , +0,+1,+1],
[-1,-1,+1 , +0,+1,+1],
[+1,+1,+1 , +0,-1,+1],
[-1,+1,+1 , +0,-1,+1],
[-1,+1,+1 , +1,+0,+1],
[-1,-1,+1 , +1,+0,+1],
[+1,+1,+1 , -1,+0,+1],
[+1,-1,+1 , -1,+0,+1],
[-1,+1,+1 , +1,+1,+0],
[-1,+1,-1 , +1,+1,+0],
[+1,+1,+1 , -1,+1,+0],
[+1,+1,-1 , -1,+1,+0],
[-1,+1,+1 , +2,+1,+1],
[+1,+1,+1 , -2,+1,+1],
[+1,+1,-1 , +2,-1,+1],
[+1,-1,+1 , +2,+1,-1],
[+1,-1,+1 , +1,+2,+1],
[+1,+1,-1 , -1,+2,+1],
[+1,+1,+1 , +1,-2,+1],
[-1,+1,+1 , +1,+2,-1],
[+1,+1,-1 , +1,+1,+2],
[+1,-1,+1 , -1,+1,+2],
[-1,+1,+1 , +1,-1,+2],
[+1,+1,+1 , +1,+1,-2],
[+1,+1,-1 , +1,+2,+3],
[+1,-1,+1 , -1,+2,+3],
[-1,+1,+1 , +1,-2,+3],
[+1,+1,+1 , +1,+2,-3],
[+1,-1,+1 , +1,+3,+2],
[+1,+1,-1 , -1,+3,+2],
[+1,+1,+1 , +1,-3,+2],
[-1,+1,+1 , +1,+3,-2],
[+1,+1,-1 , +2,+1,+3],
[+1,-1,+1 , -2,+1,+3],
[-1,+1,+1 , +2,-1,+3],
[+1,+1,+1 , +2,+1,-3],
[+1,-1,+1 , +2,+3,+1],
[+1,+1,-1 , -2,+3,+1],
[+1,+1,+1 , +2,-3,+1],
[-1,+1,+1 , +2,+3,-1],
[-1,+1,+1 , +3,+1,+2],
[+1,+1,+1 , -3,+1,+2],
[+1,+1,-1 , +3,-1,+2],
[+1,-1,+1 , +3,+1,-2],
[-1,+1,+1 , +3,+2,+1],
[+1,+1,+1 , -3,+2,+1],
[+1,+1,-1 , +3,-2,+1],
[+1,-1,+1 , +3,+2,-1],
],'d'),
'twin' : _np.array([
[-1, 1, 1, 2, 1, 1],
[ 1, 1, 1, -2, 1, 1],
[ 1, 1,-1, 2,-1, 1],
[ 1,-1, 1, 2, 1,-1],
[ 1,-1, 1, 1, 2, 1],
[ 1, 1,-1, -1, 2, 1],
[ 1, 1, 1, 1,-2, 1],
[-1, 1, 1, 1, 2,-1],
[ 1, 1,-1, 1, 1, 2],
[ 1,-1, 1, -1, 1, 2],
[-1, 1, 1, 1,-1, 2],
[ 1, 1, 1, 1, 1,-2],
],dtype=float),
},
'hP': {
'slip' : _np.array([
[+2,-1,-1,+0 , +0,+0,+0,+1],
[-1,+2,-1,+0 , +0,+0,+0,+1],
[-1,-1,+2,+0 , +0,+0,+0,+1],
[+2,-1,-1,+0 , +0,+1,-1,+0],
[-1,+2,-1,+0 , -1,+0,+1,+0],
[-1,-1,+2,+0 , +1,-1,+0,+0],
[-1,+1,+0,+0 , +1,+1,-2,+0],
[+0,-1,+1,+0 , -2,+1,+1,+0],
[+1,+0,-1,+0 , +1,-2,+1,+0],
[-1,+2,-1,+0 , +1,+0,-1,+1],
[-2,+1,+1,+0 , +0,+1,-1,+1],
[-1,-1,+2,+0 , -1,+1,+0,+1],
[+1,-2,+1,+0 , -1,+0,+1,+1],
[+2,-1,-1,+0 , +0,-1,+1,+1],
[+1,+1,-2,+0 , +1,-1,+0,+1],
[-2,+1,+1,+3 , +1,+0,-1,+1],
[-1,-1,+2,+3 , +1,+0,-1,+1],
[-1,-1,+2,+3 , +0,+1,-1,+1],
[+1,-2,+1,+3 , +0,+1,-1,+1],
[+1,-2,+1,+3 , -1,+1,+0,+1],
[+2,-1,-1,+3 , -1,+1,+0,+1],
[+2,-1,-1,+3 , -1,+0,+1,+1],
[+1,+1,-2,+3 , -1,+0,+1,+1],
[+1,+1,-2,+3 , +0,-1,+1,+1],
[-1,+2,-1,+3 , +0,-1,+1,+1],
[-1,+2,-1,+3 , +1,-1,+0,+1],
[-2,+1,+1,+3 , +1,-1,+0,+1],
[-1,-1,+2,+3 , +1,+1,-2,+2],
[+1,-2,+1,+3 , -1,+2,-1,+2],
[+2,-1,-1,+3 , -2,+1,+1,+2],
[+1,+1,-2,+3 , -1,-1,+2,+2],
[-1,+2,-1,+3 , +1,-2,+1,+2],
[-2,+1,+1,+3 , +2,-1,-1,+2],
],'d'),
'twin' : _np.array([
[-1, 0, 1, 1, 1, 0, -1, 2], # shear = (3-(c/a)^2)/(sqrt(3) c/a) <-10.1>{10.2}
[ 0, -1, 1, 1, 0, 1, -1, 2],
[ 1, -1, 0, 1, -1, 1, 0, 2],
[ 1, 0, -1, 1, -1, 0, 1, 2],
[ 0, 1, -1, 1, 0, -1, 1, 2],
[-1, 1, 0, 1, 1, -1, 0, 2],
[-1, -1, 2, 6, 1, 1, -2, 1], # shear = 1/(c/a) <11.6>{-1-1.1}
[ 1, -2, 1, 6, -1, 2, -1, 1],
[ 2, -1, -1, 6, -2, 1, 1, 1],
[ 1, 1, -2, 6, -1, -1, 2, 1],
[-1, 2, -1, 6, 1, -2, 1, 1],
[-2, 1, 1, 6, 2, -1, -1, 1],
[ 1, 0, -1, -2, 1, 0, -1, 1], # shear = (4(c/a)^2-9)/(4 sqrt(3) c/a) <10.-2>{10.1}
[ 0, 1, -1, -2, 0, 1, -1, 1],
[-1, 1, 0, -2, -1, 1, 0, 1],
[-1, 0, 1, -2, -1, 0, 1, 1],
[ 0, -1, 1, -2, 0, -1, 1, 1],
[ 1, -1, 0, -2, 1, -1, 0, 1],
[ 1, 1, -2, -3, 1, 1, -2, 2], # shear = 2((c/a)^2-2)/(3 c/a) <11.-3>{11.2}
[-1, 2, -1, -3, -1, 2, -1, 2],
[-2, 1, 1, -3, -2, 1, 1, 2],
[-1, -1, 2, -3, -1, -1, 2, 2],
[ 1, -2, 1, -3, 1, -2, 1, 2],
[ 2, -1, -1, -3, 2, -1, -1, 2],
],dtype=float),
},
}
# Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation
# from S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
# also see K. Kitahara et al., Acta Materialia 54:1279-1288, 2006
relations = {
'KS': {
'cF' : _np.array([
[[ -1, 0, 1],[ 1, 1, 1]],
[[ -1, 0, 1],[ 1, 1, 1]],
[[ 0, 1, -1],[ 1, 1, 1]],
[[ 0, 1, -1],[ 1, 1, 1]],
[[ 1, -1, 0],[ 1, 1, 1]],
[[ 1, -1, 0],[ 1, 1, 1]],
[[ 1, 0, -1],[ 1, -1, 1]],
[[ 1, 0, -1],[ 1, -1, 1]],
[[ -1, -1, 0],[ 1, -1, 1]],
[[ -1, -1, 0],[ 1, -1, 1]],
[[ 0, 1, 1],[ 1, -1, 1]],
[[ 0, 1, 1],[ 1, -1, 1]],
[[ 0, -1, 1],[ -1, 1, 1]],
[[ 0, -1, 1],[ -1, 1, 1]],
[[ -1, 0, -1],[ -1, 1, 1]],
[[ -1, 0, -1],[ -1, 1, 1]],
[[ 1, 1, 0],[ -1, 1, 1]],
[[ 1, 1, 0],[ -1, 1, 1]],
[[ -1, 1, 0],[ 1, 1, -1]],
[[ -1, 1, 0],[ 1, 1, -1]],
[[ 0, -1, -1],[ 1, 1, -1]],
[[ 0, -1, -1],[ 1, 1, -1]],
[[ 1, 0, 1],[ 1, 1, -1]],
[[ 1, 0, 1],[ 1, 1, -1]],
],dtype=float),
'cI' : _np.array([
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ -1, -1, 1],[ 0, 1, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
],dtype=float),
},
'GT': {
'cF' : _np.array([
[[ -5,-12, 17],[ 1, 1, 1]],
[[ 17, -5,-12],[ 1, 1, 1]],
[[-12, 17, -5],[ 1, 1, 1]],
[[ 5, 12, 17],[ -1, -1, 1]],
[[-17, 5,-12],[ -1, -1, 1]],
[[ 12,-17, -5],[ -1, -1, 1]],
[[ -5, 12,-17],[ -1, 1, 1]],
[[ 17, 5, 12],[ -1, 1, 1]],
[[-12,-17, 5],[ -1, 1, 1]],
[[ 5,-12,-17],[ 1, -1, 1]],
[[-17, -5, 12],[ 1, -1, 1]],
[[ 12, 17, 5],[ 1, -1, 1]],
[[ -5, 17,-12],[ 1, 1, 1]],
[[-12, -5, 17],[ 1, 1, 1]],
[[ 17,-12, -5],[ 1, 1, 1]],
[[ 5,-17,-12],[ -1, -1, 1]],
[[ 12, 5, 17],[ -1, -1, 1]],
[[-17, 12, -5],[ -1, -1, 1]],
[[ -5,-17, 12],[ -1, 1, 1]],
[[-12, 5,-17],[ -1, 1, 1]],
[[ 17, 12, 5],[ -1, 1, 1]],
[[ 5, 17, 12],[ 1, -1, 1]],
[[ 12, -5,-17],[ 1, -1, 1]],
[[-17,-12, 5],[ 1, -1, 1]],
],dtype=float),
'cI' : _np.array([
[[-17, -7, 17],[ 1, 0, 1]],
[[ 17,-17, -7],[ 1, 1, 0]],
[[ -7, 17,-17],[ 0, 1, 1]],
[[ 17, 7, 17],[ -1, 0, 1]],
[[-17, 17, -7],[ -1, -1, 0]],
[[ 7,-17,-17],[ 0, -1, 1]],
[[-17, 7,-17],[ -1, 0, 1]],
[[ 17, 17, 7],[ -1, 1, 0]],
[[ -7,-17, 17],[ 0, 1, 1]],
[[ 17, -7,-17],[ 1, 0, 1]],
[[-17,-17, 7],[ 1, -1, 0]],
[[ 7, 17, 17],[ 0, -1, 1]],
[[-17, 17, -7],[ 1, 1, 0]],
[[ -7,-17, 17],[ 0, 1, 1]],
[[ 17, -7,-17],[ 1, 0, 1]],
[[ 17,-17, -7],[ -1, -1, 0]],
[[ 7, 17, 17],[ 0, -1, 1]],
[[-17, 7,-17],[ -1, 0, 1]],
[[-17,-17, 7],[ -1, 1, 0]],
[[ -7, 17,-17],[ 0, 1, 1]],
[[ 17, 7, 17],[ -1, 0, 1]],
[[ 17, 17, 7],[ 1, -1, 0]],
[[ 7,-17,-17],[ 0, -1, 1]],
[[-17, -7, 17],[ 1, 0, 1]],
],dtype=float),
},
'GT_prime': {
'cF' : _np.array([
[[ 0, 1, -1],[ 7, 17, 17]],
[[ -1, 0, 1],[ 17, 7, 17]],
[[ 1, -1, 0],[ 17, 17, 7]],
[[ 0, -1, -1],[ -7,-17, 17]],
[[ 1, 0, 1],[-17, -7, 17]],
[[ 1, -1, 0],[-17,-17, 7]],
[[ 0, 1, -1],[ 7,-17,-17]],
[[ 1, 0, 1],[ 17, -7,-17]],
[[ -1, -1, 0],[ 17,-17, -7]],
[[ 0, -1, -1],[ -7, 17,-17]],
[[ -1, 0, 1],[-17, 7,-17]],
[[ -1, -1, 0],[-17, 17, -7]],
[[ 0, -1, 1],[ 7, 17, 17]],
[[ 1, 0, -1],[ 17, 7, 17]],
[[ -1, 1, 0],[ 17, 17, 7]],
[[ 0, 1, 1],[ -7,-17, 17]],
[[ -1, 0, -1],[-17, -7, 17]],
[[ -1, 1, 0],[-17,-17, 7]],
[[ 0, -1, 1],[ 7,-17,-17]],
[[ -1, 0, -1],[ 17, -7,-17]],
[[ 1, 1, 0],[ 17,-17, -7]],
[[ 0, 1, 1],[ -7, 17,-17]],
[[ 1, 0, -1],[-17, 7,-17]],
[[ 1, 1, 0],[-17, 17, -7]],
],dtype=float),
'cI' : _np.array([
[[ 1, 1, -1],[ 12, 5, 17]],
[[ -1, 1, 1],[ 17, 12, 5]],
[[ 1, -1, 1],[ 5, 17, 12]],
[[ -1, -1, -1],[-12, -5, 17]],
[[ 1, -1, 1],[-17,-12, 5]],
[[ 1, -1, -1],[ -5,-17, 12]],
[[ -1, 1, -1],[ 12, -5,-17]],
[[ 1, 1, 1],[ 17,-12, -5]],
[[ -1, -1, 1],[ 5,-17,-12]],
[[ 1, -1, -1],[-12, 5,-17]],
[[ -1, -1, 1],[-17, 12, -5]],
[[ -1, -1, -1],[ -5, 17,-12]],
[[ 1, -1, 1],[ 12, 17, 5]],
[[ 1, 1, -1],[ 5, 12, 17]],
[[ -1, 1, 1],[ 17, 5, 12]],
[[ -1, 1, 1],[-12,-17, 5]],
[[ -1, -1, -1],[ -5,-12, 17]],
[[ -1, 1, -1],[-17, -5, 12]],
[[ -1, -1, 1],[ 12,-17, -5]],
[[ -1, 1, -1],[ 5,-12,-17]],
[[ 1, 1, 1],[ 17, -5,-12]],
[[ 1, 1, 1],[-12, 17, -5]],
[[ 1, -1, -1],[ -5, 12,-17]],
[[ 1, 1, -1],[-17, 5,-12]],
],dtype=float),
},
'NW': {
'cF' : _np.array([
[[ 2, -1, -1],[ 1, 1, 1]],
[[ -1, 2, -1],[ 1, 1, 1]],
[[ -1, -1, 2],[ 1, 1, 1]],
[[ -2, -1, -1],[ -1, 1, 1]],
[[ 1, 2, -1],[ -1, 1, 1]],
[[ 1, -1, 2],[ -1, 1, 1]],
[[ 2, 1, -1],[ 1, -1, 1]],
[[ -1, -2, -1],[ 1, -1, 1]],
[[ -1, 1, 2],[ 1, -1, 1]],
[[ 2, -1, 1],[ -1, -1, 1]],
[[ -1, 2, 1],[ -1, -1, 1]],
[[ -1, -1, -2],[ -1, -1, 1]],
],dtype=float),
'cI' : _np.array([
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
[[ 0, -1, 1],[ 0, 1, 1]],
],dtype=float),
},
'Pitsch': {
'cF' : _np.array([
[[ 1, 0, 1],[ 0, 1, 0]],
[[ 1, 1, 0],[ 0, 0, 1]],
[[ 0, 1, 1],[ 1, 0, 0]],
[[ 0, 1, -1],[ 1, 0, 0]],
[[ -1, 0, 1],[ 0, 1, 0]],
[[ 1, -1, 0],[ 0, 0, 1]],
[[ 1, 0, -1],[ 0, 1, 0]],
[[ -1, 1, 0],[ 0, 0, 1]],
[[ 0, -1, 1],[ 1, 0, 0]],
[[ 0, 1, 1],[ 1, 0, 0]],
[[ 1, 0, 1],[ 0, 1, 0]],
[[ 1, 1, 0],[ 0, 0, 1]],
],dtype=float),
'cI' : _np.array([
[[ 1, -1, 1],[ -1, 0, 1]],
[[ 1, 1, -1],[ 1, -1, 0]],
[[ -1, 1, 1],[ 0, 1, -1]],
[[ -1, 1, -1],[ 0, -1, -1]],
[[ -1, -1, 1],[ -1, 0, -1]],
[[ 1, -1, -1],[ -1, -1, 0]],
[[ 1, -1, -1],[ -1, 0, -1]],
[[ -1, 1, -1],[ -1, -1, 0]],
[[ -1, -1, 1],[ 0, -1, -1]],
[[ -1, 1, 1],[ 0, -1, 1]],
[[ 1, -1, 1],[ 1, 0, -1]],
[[ 1, 1, -1],[ -1, 1, 0]],
],dtype=float),
},
'Bain': {
'cF' : _np.array([
[[ 0, 1, 0],[ 1, 0, 0]],
[[ 0, 0, 1],[ 0, 1, 0]],
[[ 1, 0, 0],[ 0, 0, 1]],
],dtype=float),
'cI' : _np.array([
[[ 0, 1, 1],[ 1, 0, 0]],
[[ 1, 0, 1],[ 0, 1, 0]],
[[ 1, 1, 0],[ 0, 0, 1]],
],dtype=float),
},
'Burgers' : {
'cI' : _np.array([
[[ -1, 1, 1],[ 1, 1, 0]],
[[ -1, 1, -1],[ 1, 1, 0]],
[[ 1, 1, 1],[ 1, -1, 0]],
[[ 1, 1, -1],[ 1, -1, 0]],
[[ 1, 1, -1],[ 1, 0, 1]],
[[ -1, 1, 1],[ 1, 0, 1]],
[[ 1, 1, 1],[ -1, 0, 1]],
[[ 1, -1, 1],[ -1, 0, 1]],
[[ -1, 1, -1],[ 0, 1, 1]],
[[ 1, 1, -1],[ 0, 1, 1]],
[[ -1, 1, 1],[ 0, -1, 1]],
[[ 1, 1, 1],[ 0, -1, 1]],
],dtype=float),
'hP' : _np.array([
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
[[ -1, 2, -1, 0],[ 0, 0, 0, 1]],
[[ -1, -1, 2, 0],[ 0, 0, 0, 1]],
],dtype=float),
},
}

77
python/damask/util.py
View File

@ -27,6 +27,7 @@ __all__=[
'execution_stamp', 'execution_stamp',
'shapeshifter', 'shapeblender', 'shapeshifter', 'shapeblender',
'extend_docstring', 'extended_docstring', 'extend_docstring', 'extended_docstring',
'Bravais_to_Miller', 'Miller_to_Bravais',
'DREAM3D_base_group', 'DREAM3D_cell_data_group', 'DREAM3D_base_group', 'DREAM3D_cell_data_group',
'dict_prune', 'dict_flatten' 'dict_prune', 'dict_flatten'
] ]
@ -286,6 +287,8 @@ def project_stereographic(vector,direction='z',normalize=True,keepdims=False):
Examples Examples
-------- --------
>>> import damask
>>> import numpy as np
>>> project_stereographic(np.ones(3)) >>> project_stereographic(np.ones(3))
[0.3660254, 0.3660254] [0.3660254, 0.3660254]
>>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True) >>> project_stereographic(np.ones(3),direction='x',normalize=False,keepdims=True)
@ -338,7 +341,7 @@ def hybrid_IA(dist,N,rng_seed=None):
def shapeshifter(fro,to,mode='left',keep_ones=False): def shapeshifter(fro,to,mode='left',keep_ones=False):
""" """
Return a tuple that reshapes 'fro' to become broadcastable to 'to'. Return dimensions that reshape 'fro' to become broadcastable to 'to'.
Parameters Parameters
---------- ----------
@ -355,6 +358,22 @@ def shapeshifter(fro,to,mode='left',keep_ones=False):
Treat '1' in fro as literal value instead of dimensional placeholder. Treat '1' in fro as literal value instead of dimensional placeholder.
Defaults to False. Defaults to False.
Returns
-------
new_dims : tuple
Dimensions for reshape.
Example
-------
>>> import numpy as np
>>> from damask import util
>>> a = np.ones((3,4,2))
>>> b = np.ones(4)
>>> b_extended = b.reshape(util.shapeshifter(b.shape,a.shape))
>>> (a * np.broadcast_to(b_extended,a.shape)).shape
(3,4,2)
""" """
beg = dict(left ='(^.*\\b)', beg = dict(left ='(^.*\\b)',
right='(^.*?\\b)') right='(^.*?\\b)')
@ -499,6 +518,62 @@ def DREAM3D_cell_data_group(fname):
return cell_data_group return cell_data_group
def Bravais_to_Miller(*,uvtw=None,hkil=None):
"""
Transform 4 Miller–Bravais indices to 3 Miller indices of crystal direction [uvw] or plane normal (hkl).
Parameters
----------
uvtw|hkil : numpy.ndarray of shape (...,4)
Miller–Bravais indices of crystallographic direction [uvtw] or plane normal (hkil).
Returns
-------
uvw|hkl : numpy.ndarray of shape (...,3)
Miller indices of [uvw] direction or (hkl) plane normal.
"""
if (uvtw is not None) ^ (hkil is None):
raise KeyError('Specify either "uvtw" or "hkil"')
axis,basis = (np.array(uvtw),np.array([[1,0,-1,0],
[0,1,-1,0],
[0,0, 0,1]])) \
if hkil is None else \
(np.array(hkil),np.array([[1,0,0,0],
[0,1,0,0],
[0,0,0,1]]))
return np.einsum('il,...l',basis,axis)
def Miller_to_Bravais(*,uvw=None,hkl=None):
"""
Transform 3 Miller indices to 4 Miller–Bravais indices of crystal direction [uvtw] or plane normal (hkil).
Parameters
----------
uvw|hkl : numpy.ndarray of shape (...,3)
Miller indices of crystallographic direction [uvw] or plane normal (hkl).
Returns
-------
uvtw|hkil : numpy.ndarray of shape (...,4)
Miller–Bravais indices of [uvtw] direction or (hkil) plane normal.
"""
if (uvw is not None) ^ (hkl is None):
raise KeyError('Specify either "uvw" or "hkl"')
axis,basis = (np.array(uvw),np.array([[ 2,-1, 0],
[-1, 2, 0],
[-1,-1, 0],
[ 0, 0, 3]])/3) \
if hkl is None else \
(np.array(hkl),np.array([[ 1, 0, 0],
[ 0, 1, 0],
[-1,-1, 0],
[ 0, 0, 1]]))
return np.einsum('il,...l',basis,axis)
def dict_prune(d): def dict_prune(d):
""" """
Recursively remove empty dictionaries. Recursively remove empty dictionaries.

57
python/tests/test_Crystal.py Normal file
View File

@ -0,0 +1,57 @@
import pytest
import numpy as np
from damask import Crystal
class TestCrystal:
def test_double_to_lattice(self):
c = Crystal(lattice='cF')
with pytest.raises(KeyError):
c.to_lattice(direction=np.ones(3),plane=np.ones(3))
def test_double_to_frame(self):
c = Crystal(lattice='cF')
with pytest.raises(KeyError):
c.to_frame(uvw=np.ones(3),hkl=np.ones(3))
@pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma',
[
('aP',0.5,2.0,3.0,0.8,0.5,1.2),
('mP',1.0,2.0,3.0,np.pi/2,0.5,np.pi/2),
('oI',0.5,1.5,3.0,np.pi/2,np.pi/2,np.pi/2),
('tP',0.5,0.5,3.0,np.pi/2,np.pi/2,np.pi/2),
('hP',1.0,None,1.6,np.pi/2,np.pi/2,2*np.pi/3),
('cF',1.0,1.0,None,np.pi/2,np.pi/2,np.pi/2),
])
def test_bases_contraction(self,lattice,a,b,c,alpha,beta,gamma):
c = Crystal(lattice=lattice,
a=a,b=b,c=c,
alpha=alpha,beta=beta,gamma=gamma)
assert np.allclose(np.eye(3),np.einsum('ik,jk',c.basis_real,c.basis_reciprocal))
@pytest.mark.parametrize('keyFrame,keyLattice',[('uvw','direction'),('hkl','plane'),])
@pytest.mark.parametrize('vector',np.array([
[1.,1.,1.],
[-2.,3.,0.5],
[0.,0.,1.],
[1.,1.,1.],
[2.,2.,2.],
[0.,1.,1.],
]))
@pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma',
[
('aP',0.5,2.0,3.0,0.8,0.5,1.2),
('mP',1.0,2.0,3.0,np.pi/2,0.5,np.pi/2),
('oI',0.5,1.5,3.0,np.pi/2,np.pi/2,np.pi/2),
('tP',0.5,0.5,3.0,np.pi/2,np.pi/2,np.pi/2),
('hP',1.0,1.0,1.6,np.pi/2,np.pi/2,2*np.pi/3),
('cF',1.0,1.0,1.0,np.pi/2,np.pi/2,np.pi/2),
])
def test_to_frame_to_lattice(self,lattice,a,b,c,alpha,beta,gamma,vector,keyFrame,keyLattice):
c = Crystal(lattice=lattice,
a=a,b=b,c=c,
alpha=alpha,beta=beta,gamma=gamma)
assert np.allclose(vector,
c.to_frame(**{keyFrame:c.to_lattice(**{keyLattice:vector})}))

277
python/tests/test_Orientation.py
View File

@ -7,7 +7,6 @@ from damask import Orientation
from damask import Table from damask import Table
from damask import util from damask import util
from damask import grid_filters from damask import grid_filters
from damask import lattice
from damask import _orientation from damask import _orientation
crystal_families = set(_orientation.lattice_symmetries.values()) crystal_families = set(_orientation.lattice_symmetries.values())
@ -25,38 +24,42 @@ def set_of_rodrigues(set_of_quaternions):
class TestOrientation: class TestOrientation:
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('shape',[None,5,(4,6)]) @pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_equal(self,lattice,shape): def test_equal(self,family,shape):
R = Rotation.from_random(shape) R = Rotation.from_random(shape)
assert Orientation(R,lattice) == Orientation(R,lattice) if shape is None else \ assert Orientation(R,family=family) == Orientation(R,family=family) if shape is None else \
(Orientation(R,lattice) == Orientation(R,lattice)).all() (Orientation(R,family=family) == Orientation(R,family=family)).all()
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('shape',[None,5,(4,6)]) @pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_unequal(self,lattice,shape): def test_unequal(self,family,shape):
R = Rotation.from_random(shape) R = Rotation.from_random(shape)
assert not ( Orientation(R,lattice) != Orientation(R,lattice) if shape is None else \ assert not ( Orientation(R,family=family) != Orientation(R,family=family) if shape is None else \
(Orientation(R,lattice) != Orientation(R,lattice)).any()) (Orientation(R,family=family) != Orientation(R,family=family)).any())
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('shape',[None,5,(4,6)]) @pytest.mark.parametrize('shape',[None,5,(4,6)])
def test_close(self,lattice,shape): def test_close(self,family,shape):
R = Orientation.from_random(lattice=lattice,shape=shape) R = Orientation.from_random(family=family,shape=shape)
assert R.isclose(R.reduced).all() and R.allclose(R.reduced) assert R.isclose(R.reduced).all() and R.allclose(R.reduced)
@pytest.mark.parametrize('a,b',[ @pytest.mark.parametrize('a,b',[
(dict(rotation=[1,0,0,0],lattice='triclinic'), (dict(rotation=[1,0,0,0],family='triclinic'),
dict(rotation=[0.5,0.5,0.5,0.5],lattice='triclinic')), dict(rotation=[0.5,0.5,0.5,0.5],family='triclinic')),
(dict(rotation=[1,0,0,0],lattice='cubic'), (dict(rotation=[1,0,0,0],family='cubic'),
dict(rotation=[1,0,0,0],lattice='hexagonal')), dict(rotation=[1,0,0,0],family='hexagonal')),
])
def test_unequal_family(self,a,b):
assert Orientation(**a) != Orientation(**b)
@pytest.mark.parametrize('a,b',[
(dict(rotation=[1,0,0,0],lattice='cF',a=1), (dict(rotation=[1,0,0,0],lattice='cF',a=1),
dict(rotation=[1,0,0,0],lattice='cF',a=2)), dict(rotation=[1,0,0,0],lattice='cF',a=2)),
]) ])
def test_nonequal(self,a,b): def test_unequal_lattice(self,a,b):
assert Orientation(**a) != Orientation(**b) assert Orientation(**a) != Orientation(**b)
@pytest.mark.parametrize('kwargs',[ @pytest.mark.parametrize('kwargs',[
@ -72,7 +75,7 @@ class TestOrientation:
]) ])
def test_invalid_init(self,kwargs): def test_invalid_init(self,kwargs):
with pytest.raises(ValueError): with pytest.raises(ValueError):
Orientation(**kwargs).parameters # noqa Orientation(**kwargs)
@pytest.mark.parametrize('kwargs',[ @pytest.mark.parametrize('kwargs',[
dict(lattice='aP',a=1.0,b=1.1,c=1.2,alpha=np.pi/4,beta=np.pi/3,gamma=np.pi/2), dict(lattice='aP',a=1.0,b=1.1,c=1.2,alpha=np.pi/4,beta=np.pi/3,gamma=np.pi/2),
@ -100,47 +103,47 @@ class TestOrientation:
assert o != p assert o != p
def test_from_quaternion(self): def test_from_quaternion(self):
assert np.all(Orientation.from_quaternion(q=np.array([1,0,0,0]),lattice='triclinic').as_matrix() assert np.all(Orientation.from_quaternion(q=np.array([1,0,0,0]),family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_Euler_angles(self): def test_from_Euler_angles(self):
assert np.all(Orientation.from_Euler_angles(phi=np.zeros(3),lattice='triclinic').as_matrix() assert np.all(Orientation.from_Euler_angles(phi=np.zeros(3),family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_axis_angle(self): def test_from_axis_angle(self):
assert np.all(Orientation.from_axis_angle(axis_angle=[1,0,0,0],lattice='triclinic').as_matrix() assert np.all(Orientation.from_axis_angle(axis_angle=[1,0,0,0],family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_basis(self): def test_from_basis(self):
assert np.all(Orientation.from_basis(basis=np.eye(3),lattice='triclinic').as_matrix() assert np.all(Orientation.from_basis(basis=np.eye(3),family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_matrix(self): def test_from_matrix(self):
assert np.all(Orientation.from_matrix(R=np.eye(3),lattice='triclinic').as_matrix() assert np.all(Orientation.from_matrix(R=np.eye(3),family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_Rodrigues_vector(self): def test_from_Rodrigues_vector(self):
assert np.all(Orientation.from_Rodrigues_vector(rho=np.array([0,0,1,0]),lattice='triclinic').as_matrix() assert np.all(Orientation.from_Rodrigues_vector(rho=np.array([0,0,1,0]),family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_homochoric(self): def test_from_homochoric(self):
assert np.all(Orientation.from_homochoric(h=np.zeros(3),lattice='triclinic').as_matrix() assert np.all(Orientation.from_homochoric(h=np.zeros(3),family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_cubochoric(self): def test_from_cubochoric(self):
assert np.all(Orientation.from_cubochoric(x=np.zeros(3),lattice='triclinic').as_matrix() assert np.all(Orientation.from_cubochoric(x=np.zeros(3),family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_spherical_component(self): def test_from_spherical_component(self):
assert np.all(Orientation.from_spherical_component(center=Rotation(), assert np.all(Orientation.from_spherical_component(center=Rotation(),
sigma=0.0,N=1,lattice='triclinic').as_matrix() sigma=0.0,N=1,family='triclinic').as_matrix()
== np.eye(3)) == np.eye(3))
def test_from_fiber_component(self): def test_from_fiber_component(self):
r = Rotation.from_fiber_component(alpha=np.zeros(2),beta=np.zeros(2), r = Rotation.from_fiber_component(alpha=np.zeros(2),beta=np.zeros(2),
sigma=0.0,N=1,rng_seed=0) sigma=0.0,N=1,rng_seed=0)
assert np.all(Orientation.from_fiber_component(alpha=np.zeros(2),beta=np.zeros(2), assert np.all(Orientation.from_fiber_component(alpha=np.zeros(2),beta=np.zeros(2),
sigma=0.0,N=1,rng_seed=0,lattice='triclinic').quaternion sigma=0.0,N=1,rng_seed=0,family='triclinic').quaternion
== r.quaternion) == r.quaternion)
@pytest.mark.parametrize('kwargs',[ @pytest.mark.parametrize('kwargs',[
@ -185,26 +188,26 @@ class TestOrientation:
with pytest.raises(ValueError): with pytest.raises(ValueError):
Orientation(lattice='aP',a=1,b=2,c=3,alpha=45,beta=45,gamma=90.0001,degrees=True) # noqa Orientation(lattice='aP',a=1,b=2,c=3,alpha=45,beta=45,gamma=90.0001,degrees=True) # noqa
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('angle',[10,20,30,40]) @pytest.mark.parametrize('angle',[10,20,30,40])
def test_average(self,angle,lattice): def test_average(self,angle,family):
o = Orientation.from_axis_angle(lattice=lattice,axis_angle=[[0,0,1,10],[0,0,1,angle]],degrees=True) o = Orientation.from_axis_angle(family=family,axis_angle=[[0,0,1,10],[0,0,1,angle]],degrees=True)
avg_angle = o.average().as_axis_angle(degrees=True,pair=True)[1] avg_angle = o.average().as_axis_angle(degrees=True,pair=True)[1]
assert np.isclose(avg_angle,10+(angle-10)/2.) assert np.isclose(avg_angle,10+(angle-10)/2.)
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
def test_reduced_equivalent(self,lattice): def test_reduced_equivalent(self,family):
i = Orientation(lattice=lattice) i = Orientation(family=family)
o = Orientation.from_random(lattice=lattice) o = Orientation.from_random(family=family)
eq = o.equivalent eq = o.equivalent
FZ = np.argmin(abs(eq.misorientation(i.broadcast_to(len(eq))).as_axis_angle(pair=True)[1])) FZ = np.argmin(abs(eq.misorientation(i.broadcast_to(len(eq))).as_axis_angle(pair=True)[1]))
assert o.reduced == eq[FZ] assert o.reduced == eq[FZ]
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('N',[1,8,32]) @pytest.mark.parametrize('N',[1,8,32])
def test_disorientation(self,lattice,N): def test_disorientation(self,family,N):
o = Orientation.from_random(lattice=lattice,shape=N) o = Orientation.from_random(family=family,shape=N)
p = Orientation.from_random(lattice=lattice,shape=N) p = Orientation.from_random(family=family,shape=N)
d,ops = o.disorientation(p,return_operators=True) d,ops = o.disorientation(p,return_operators=True)
@ -218,72 +221,72 @@ class TestOrientation:
.misorientation(p[n].equivalent[ops[n][1]]) .misorientation(p[n].equivalent[ops[n][1]])
.as_quaternion()) .as_quaternion())
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('a,b',[ @pytest.mark.parametrize('a,b',[
((2,3,2),(2,3,2)), ((2,3,2),(2,3,2)),
((2,2),(4,4)), ((2,2),(4,4)),
((3,1),(1,3)), ((3,1),(1,3)),
(None,None), (None,None),
]) ])
def test_disorientation_blending(self,lattice,a,b): def test_disorientation_blending(self,family,a,b):
o = Orientation.from_random(lattice=lattice,shape=a) o = Orientation.from_random(family=family,shape=a)
p = Orientation.from_random(lattice=lattice,shape=b) p = Orientation.from_random(family=family,shape=b)
blend = util.shapeblender(o.shape,p.shape) blend = util.shapeblender(o.shape,p.shape)
for loc in np.random.randint(0,blend,(10,len(blend))): for loc in np.random.randint(0,blend,(10,len(blend))):
assert o[tuple(loc[:len(o.shape)])].disorientation(p[tuple(loc[-len(p.shape):])]) \ assert o[tuple(loc[:len(o.shape)])].disorientation(p[tuple(loc[-len(p.shape):])]) \
.isclose(o.disorientation(p)[tuple(loc)]) .isclose(o.disorientation(p)[tuple(loc)])
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
def test_disorientation360(self,lattice): def test_disorientation360(self,family):
o_1 = Orientation(Rotation(),lattice) o_1 = Orientation(Rotation(),family=family)
o_2 = Orientation.from_Euler_angles(lattice=lattice,phi=[360,0,0],degrees=True) o_2 = Orientation.from_Euler_angles(family=family,phi=[360,0,0],degrees=True)
assert np.allclose((o_1.disorientation(o_2)).as_matrix(),np.eye(3)) assert np.allclose((o_1.disorientation(o_2)).as_matrix(),np.eye(3))
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)]) @pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)])
def test_reduced_vectorization(self,lattice,shape): def test_reduced_vectorization(self,family,shape):
o = Orientation.from_random(lattice=lattice,shape=shape) o = Orientation.from_random(family=family,shape=shape)
for r, theO in zip(o.reduced.flatten(),o.flatten()): for r, theO in zip(o.reduced.flatten(),o.flatten()):
assert r == theO.reduced assert r == theO.reduced
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
def test_reduced_corner_cases(self,lattice): def test_reduced_corner_cases(self,family):
# test whether there is always a sym-eq rotation that falls into the FZ # test whether there is always a sym-eq rotation that falls into the FZ
N = np.random.randint(10,40) N = np.random.randint(10,40)
size = np.ones(3)*np.pi**(2./3.) size = np.ones(3)*np.pi**(2./3.)
grid = grid_filters.coordinates0_node([N+1,N+1,N+1],size,-size*.5) grid = grid_filters.coordinates0_node([N+1,N+1,N+1],size,-size*.5)
evenly_distributed = Orientation.from_cubochoric(x=grid[:-2,:-2,:-2],lattice=lattice) evenly_distributed = Orientation.from_cubochoric(x=grid[:-2,:-2,:-2],family=family)
assert evenly_distributed.shape == evenly_distributed.reduced.shape assert evenly_distributed.shape == evenly_distributed.reduced.shape
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)]) @pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)])
@pytest.mark.parametrize('vector',np.array([[1,0,0],[1,2,3],[-1,1,-1]])) @pytest.mark.parametrize('vector',np.array([[1,0,0],[1,2,3],[-1,1,-1]]))
@pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('proper',[True,False])
def test_to_SST_vectorization(self,lattice,shape,vector,proper): def test_to_SST_vectorization(self,family,shape,vector,proper):
o = Orientation.from_random(lattice=lattice,shape=shape) o = Orientation.from_random(family=family,shape=shape)
for r, theO in zip(o.to_SST(vector=vector,proper=proper).reshape((-1,3)),o.flatten()): for r, theO in zip(o.to_SST(vector=vector,proper=proper).reshape((-1,3)),o.flatten()):
assert np.allclose(r,theO.to_SST(vector=vector,proper=proper)) assert np.allclose(r,theO.to_SST(vector=vector,proper=proper))
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)]) @pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)])
@pytest.mark.parametrize('vector',np.array([[1,0,0],[1,2,3],[-1,1,-1]])) @pytest.mark.parametrize('vector',np.array([[1,0,0],[1,2,3],[-1,1,-1]]))
@pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('proper',[True,False])
@pytest.mark.parametrize('in_SST',[True,False]) @pytest.mark.parametrize('in_SST',[True,False])
def test_IPF_color_vectorization(self,lattice,shape,vector,proper,in_SST): def test_IPF_color_vectorization(self,family,shape,vector,proper,in_SST):
o = Orientation.from_random(lattice=lattice,shape=shape) o = Orientation.from_random(family=family,shape=shape)
for r, theO in zip(o.IPF_color(vector,in_SST=in_SST,proper=proper).reshape((-1,3)),o.flatten()): for r, theO in zip(o.IPF_color(vector,in_SST=in_SST,proper=proper).reshape((-1,3)),o.flatten()):
assert np.allclose(r,theO.IPF_color(vector,in_SST=in_SST,proper=proper)) assert np.allclose(r,theO.IPF_color(vector,in_SST=in_SST,proper=proper))
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('a,b',[ @pytest.mark.parametrize('a,b',[
((2,3,2),(2,3,2)), ((2,3,2),(2,3,2)),
((2,2),(4,4)), ((2,2),(4,4)),
((3,1),(1,3)), ((3,1),(1,3)),
(None,(3,)), (None,(3,)),
]) ])
def test_to_SST_blending(self,lattice,a,b): def test_to_SST_blending(self,family,a,b):
o = Orientation.from_random(lattice=lattice,shape=a) o = Orientation.from_random(family=family,shape=a)
v = np.random.random(b+(3,)) v = np.random.random(b+(3,))
blend = util.shapeblender(o.shape,b) blend = util.shapeblender(o.shape,b)
for loc in np.random.randint(0,blend,(10,len(blend))): for loc in np.random.randint(0,blend,(10,len(blend))):
@ -298,75 +301,45 @@ class TestOrientation:
{'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}]) {'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
@pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('proper',[True,False])
def test_IPF_cubic(self,color,proper): def test_IPF_cubic(self,color,proper):
cube = Orientation(lattice='cubic') cube = Orientation(family='cubic')
for direction in set(permutations(np.array(color['direction']))): for direction in set(permutations(np.array(color['direction']))):
assert np.allclose(np.array(color['RGB']), assert np.allclose(np.array(color['RGB']),
cube.IPF_color(vector=np.array(direction),proper=proper)) cube.IPF_color(vector=np.array(direction),proper=proper))
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('proper',[True,False])
def test_IPF_equivalent(self,set_of_quaternions,lattice,proper): def test_IPF_equivalent(self,set_of_quaternions,family,proper):
direction = np.random.random(3)*2.0-1.0 direction = np.random.random(3)*2.0-1.0
o = Orientation(rotation=set_of_quaternions,lattice=lattice).equivalent o = Orientation(rotation=set_of_quaternions,family=family).equivalent
color = o.IPF_color(vector=direction,proper=proper) color = o.IPF_color(vector=direction,proper=proper)
assert np.allclose(np.broadcast_to(color[0,...],color.shape),color) assert np.allclose(np.broadcast_to(color[0,...],color.shape),color)
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
def test_in_FZ_vectorization(self,set_of_rodrigues,lattice): def test_in_FZ_vectorization(self,set_of_rodrigues,family):
result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)),lattice=lattice).in_FZ.reshape(-1) result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)),family=family).in_FZ.reshape(-1)
for r,rho in zip(result,set_of_rodrigues[:len(result)]): for r,rho in zip(result,set_of_rodrigues[:len(result)]):
assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_FZ assert r == Orientation.from_Rodrigues_vector(rho=rho,family=family).in_FZ
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
def test_in_disorientation_FZ_vectorization(self,set_of_rodrigues,lattice): def test_in_disorientation_FZ_vectorization(self,set_of_rodrigues,family):
result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)), result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)),
lattice=lattice).in_disorientation_FZ.reshape(-1) family=family).in_disorientation_FZ.reshape(-1)
for r,rho in zip(result,set_of_rodrigues[:len(result)]): for r,rho in zip(result,set_of_rodrigues[:len(result)]):
assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_disorientation_FZ assert r == Orientation.from_Rodrigues_vector(rho=rho,family=family).in_disorientation_FZ
@pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('proper',[True,False])
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
def test_in_SST_vectorization(self,lattice,proper): def test_in_SST_vectorization(self,family,proper):
vecs = np.random.rand(20,4,3) vecs = np.random.rand(20,4,3)
result = Orientation(lattice=lattice).in_SST(vecs,proper).flatten() result = Orientation(family=family).in_SST(vecs,proper).flatten()
for r,v in zip(result,vecs.reshape((-1,3))): for r,v in zip(result,vecs.reshape((-1,3))):
assert np.all(r == Orientation(lattice=lattice).in_SST(v,proper)) assert np.all(r == Orientation(family=family).in_SST(v,proper))
@pytest.mark.parametrize('invalid_lattice',['fcc','bcc','hello'])
def test_invalid_lattice_init(self,invalid_lattice):
with pytest.raises(KeyError):
Orientation(lattice=invalid_lattice) # noqa
@pytest.mark.parametrize('invalid_family',[None,'fcc','bcc','hello']) @pytest.mark.parametrize('invalid_family',[None,'fcc','bcc','hello'])
def test_invalid_symmetry_family(self,invalid_family): def test_invalid_lattice_init(self,invalid_family):
with pytest.raises(KeyError): with pytest.raises(KeyError):
o = Orientation(lattice='cubic') Orientation(family=invalid_family)
o.family = invalid_family
o.symmetry_operations # noqa
def test_invalid_rot(self):
with pytest.raises(TypeError):
Orientation.from_random(lattice='cubic') * np.ones(3)
def test_missing_symmetry_immutable(self):
with pytest.raises(KeyError):
Orientation(lattice=None).immutable # noqa
def test_missing_symmetry_basis_real(self):
with pytest.raises(KeyError):
Orientation(lattice=None).basis_real # noqa
def test_missing_symmetry_basis_reciprocal(self):
with pytest.raises(KeyError):
Orientation(lattice=None).basis_reciprocal # noqa
def test_double_to_lattice(self):
with pytest.raises(KeyError):
Orientation().to_lattice(direction=np.ones(3),plane=np.ones(3)) # noqa
def test_double_to_frame(self):
with pytest.raises(KeyError):
Orientation().to_frame(uvw=np.ones(3),hkl=np.ones(3)) # noqa
@pytest.mark.parametrize('relation',[None,'Peter','Paul']) @pytest.mark.parametrize('relation',[None,'Peter','Paul'])
def test_unknown_relation(self,relation): def test_unknown_relation(self,relation):
@ -388,28 +361,22 @@ class TestOrientation:
a=a,b=b,c=c, a=a,b=b,c=c,
alpha=alpha,beta=beta,gamma=gamma).related(relation) # noqa alpha=alpha,beta=beta,gamma=gamma).related(relation) # noqa
@pytest.mark.parametrize('lattice',crystal_families) @pytest.mark.parametrize('family',crystal_families)
@pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('proper',[True,False])
def test_in_SST(self,lattice,proper): def test_in_SST(self,family,proper):
assert Orientation(lattice=lattice).in_SST(np.zeros(3),proper) assert Orientation(family=family).in_SST(np.zeros(3),proper)
@pytest.mark.parametrize('function',['in_SST','IPF_color']) @pytest.mark.parametrize('function',['in_SST','IPF_color'])
def test_invalid_argument(self,function): def test_invalid_argument(self,function):
o = Orientation(lattice='cubic') # noqa o = Orientation(family='cubic') # noqa
with pytest.raises(ValueError): with pytest.raises(ValueError):
eval(f'o.{function}(np.ones(4))') eval(f'o.{function}(np.ones(4))')
@pytest.mark.parametrize('model',lattice.relations)
def test_relationship_definition(self,model):
m,o = list(lattice.relations[model])
assert lattice.relations[model][m].shape[:-1] == lattice.relations[model][o].shape[:-1]
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['cF','cI']) @pytest.mark.parametrize('lattice',['cF','cI'])
def test_relationship_vectorize(self,set_of_quaternions,lattice,model): def test_relationship_vectorize(self,set_of_quaternions,lattice,model):
r = Orientation(rotation=set_of_quaternions[:200].reshape((50,4,4)),lattice=lattice).related(model) r = Orientation(rotation=set_of_quaternions[:200].reshape((50,4,4)),lattice=lattice).related(model)
for i in range(200): for i in range(200):
assert (r.reshape((-1,200))[:,i] == Orientation(set_of_quaternions[i],lattice).related(model)).all() assert (r.reshape((-1,200))[:,i] == Orientation(set_of_quaternions[i],lattice=lattice).related(model)).all()
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['cF','cI']) @pytest.mark.parametrize('lattice',['cF','cI'])
@ -444,45 +411,6 @@ class TestOrientation:
) )
assert np.allclose(o.to_frame(uvw=np.eye(3)),basis), 'Lattice basis disagrees with initialization' assert np.allclose(o.to_frame(uvw=np.eye(3)),basis), 'Lattice basis disagrees with initialization'
@pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma',
[
('aP',0.5,2.0,3.0,0.8,0.5,1.2),
('mP',1.0,2.0,3.0,np.pi/2,0.5,np.pi/2),
('oI',0.5,1.5,3.0,np.pi/2,np.pi/2,np.pi/2),
('tP',0.5,0.5,3.0,np.pi/2,np.pi/2,np.pi/2),
('hP',1.0,None,1.6,np.pi/2,np.pi/2,2*np.pi/3),
('cF',1.0,1.0,None,np.pi/2,np.pi/2,np.pi/2),
])
def test_bases_contraction(self,lattice,a,b,c,alpha,beta,gamma):
L = Orientation(lattice=lattice,
a=a,b=b,c=c,
alpha=alpha,beta=beta,gamma=gamma)
assert np.allclose(np.eye(3),np.einsum('ik,jk',L.basis_real,L.basis_reciprocal))
@pytest.mark.parametrize('keyFrame,keyLattice',[('uvw','direction'),('hkl','plane'),])
@pytest.mark.parametrize('vector',np.array([
[1.,1.,1.],
[-2.,3.,0.5],
[0.,0.,1.],
[1.,1.,1.],
[2.,2.,2.],
[0.,1.,1.],
]))
@pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma',
[
('aP',0.5,2.0,3.0,0.8,0.5,1.2),
('mP',1.0,2.0,3.0,np.pi/2,0.5,np.pi/2),
('oI',0.5,1.5,3.0,np.pi/2,np.pi/2,np.pi/2),
('tP',0.5,0.5,3.0,np.pi/2,np.pi/2,np.pi/2),
('hP',1.0,1.0,1.6,np.pi/2,np.pi/2,2*np.pi/3),
('cF',1.0,1.0,1.0,np.pi/2,np.pi/2,np.pi/2),
])
def test_to_frame_to_lattice(self,lattice,a,b,c,alpha,beta,gamma,vector,keyFrame,keyLattice):
L = Orientation(lattice=lattice,
a=a,b=b,c=c,
alpha=alpha,beta=beta,gamma=gamma)
assert np.allclose(vector,
L.to_frame(**{keyFrame:L.to_lattice(**{keyLattice:vector})}))
@pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma', @pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma',
[ [
@ -510,13 +438,22 @@ class TestOrientation:
assert o.to_pole(**{kw:vector,'with_symmetry':with_symmetry}).shape \ assert o.to_pole(**{kw:vector,'with_symmetry':with_symmetry}).shape \
== o.shape + (o.symmetry_operations.shape if with_symmetry else ()) + vector.shape == o.shape + (o.symmetry_operations.shape if with_symmetry else ()) + vector.shape
@pytest.mark.parametrize('lattice',['hP','cI','cF']) @pytest.mark.parametrize('lattice',['hP','cI','cF']) #tI not included yet
@pytest.mark.parametrize('mode',['slip','twin']) def test_Schmid(self,update,ref_path,lattice):
def test_Schmid(self,update,ref_path,lattice,mode): O = Orientation(lattice=lattice) # noqa
L = Orientation(lattice=lattice) for mode in ['slip','twin']:
reference = ref_path/f'{lattice}_{mode}.txt' reference = ref_path/f'{lattice}_{mode}.txt'
P = L.Schmid(mode) P = O.Schmid(N_slip='*') if mode == 'slip' else O.Schmid(N_twin='*')
if update: if update:
table = Table(P.reshape(-1,9),{'Schmid':(3,3,)}) table = Table(P.reshape(-1,9),{'Schmid':(3,3,)})
table.save(reference) table.save(reference)
assert np.allclose(P,Table.load(reference).get('Schmid')) assert np.allclose(P,Table.load(reference).get('Schmid'))
@pytest.mark.parametrize('lattice',['hP','cI','cF']) # tI not included yet
def test_Schmid_vectorize(self,lattice):
O = Orientation.from_random(shape=4,lattice=lattice) # noqa
for mode in ['slip','twin']:
Ps = O.Schmid(N_slip='*') if mode == 'slip' else O.Schmid(N_twin='*')
for i in range(4):
P = O[i].Schmid(N_slip='*') if mode == 'slip' else O[i].Schmid(N_twin='*')
assert np.allclose(P,Ps[:,i])

21
python/tests/test_Result.py
View File

@ -12,7 +12,6 @@ import vtk
import numpy as np import numpy as np
from damask import Result from damask import Result
from damask import Rotation
from damask import Orientation from damask import Orientation
from damask import tensor from damask import tensor
from damask import mechanics from damask import mechanics
@ -220,17 +219,15 @@ class TestResult:
in_file = default.place('S') in_file = default.place('S')
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@pytest.mark.skip(reason='requires rework of lattice.f90') @pytest.mark.parametrize('options',[{'uvw':[1,0,0]},{'hkl':[0,1,1]}])
@pytest.mark.parametrize('polar',[True,False]) def test_add_pole(self,default,options):
def test_add_pole(self,default,polar): default.add_pole(**options)
pole = np.array([1.,0.,0.]) rot = default.place('O')
default.add_pole('O',pole,polar) in_memory = Orientation(rot,lattice=rot.dtype.metadata['lattice']).to_pole(**options)
rot = Rotation(default.place('O')) brackets = ['[[]','[]]'] if 'uvw' in options.keys() else ['(',')'] # escape fnmatch
rotated_pole = rot * np.broadcast_to(pole,rot.shape+(3,)) label = '{}{} {} {}{}'.format(brackets[0],*(list(options.values())[0]),brackets[1])
xy = rotated_pole[:,0:2]/(1.+abs(pole[2])) in_file = default.place(f'p^{label}')
in_memory = xy if not polar else \ print(in_file - in_memory)
np.block([np.sqrt(xy[:,0:1]*xy[:,0:1]+xy[:,1:2]*xy[:,1:2]),np.arctan2(xy[:,1:2],xy[:,0:1])])
in_file = default.place('p^{}_[1 0 0)'.format(u'rφ' if polar else 'xy'))
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
def test_add_rotation(self,default): def test_add_rotation(self,default):

34
python/tests/test_lattice.py
View File

@ -1,34 +0,0 @@
import pytest
import numpy as np
from damask import lattice
class TestLattice:
def test_double_Bravais_to_Miller(self):
with pytest.raises(KeyError):
lattice.Bravais_to_Miller(uvtw=np.ones(4),hkil=np.ones(4)) # noqa
def test_double_Miller_to_Bravais(self):
with pytest.raises(KeyError):
lattice.Miller_to_Bravais(uvw=np.ones(4),hkl=np.ones(4)) # noqa
@pytest.mark.parametrize('vector',np.array([
[1,0,0],
[1,1,0],
[1,1,1],
[1,0,-2],
]))
@pytest.mark.parametrize('kw_Miller,kw_Bravais',[('uvw','uvtw'),('hkl','hkil')])
def test_Miller_Bravais_Miller(self,vector,kw_Miller,kw_Bravais):
assert np.all(vector == lattice.Bravais_to_Miller(**{kw_Bravais:lattice.Miller_to_Bravais(**{kw_Miller:vector})}))
@pytest.mark.parametrize('vector',np.array([
[1,0,-1,2],
[1,-1,0,3],
[1,1,-2,-3],
[0,0,0,1],
]))
@pytest.mark.parametrize('kw_Miller,kw_Bravais',[('uvw','uvtw'),('hkl','hkil')])
def test_Bravais_Miller_Bravais(self,vector,kw_Miller,kw_Bravais):
assert np.all(vector == lattice.Miller_to_Bravais(**{kw_Miller:lattice.Bravais_to_Miller(**{kw_Bravais:vector})}))

30
python/tests/test_util.py
View File

@ -158,3 +158,33 @@ class TestUtil:
({'A':{'B':{},'C':'D'}}, {'B':{},'C':'D'})]) ({'A':{'B':{},'C':'D'}}, {'B':{},'C':'D'})])
def test_flatten(self,full,reduced): def test_flatten(self,full,reduced):
assert util.dict_flatten(full) == reduced assert util.dict_flatten(full) == reduced
def test_double_Bravais_to_Miller(self):
with pytest.raises(KeyError):
util.Bravais_to_Miller(uvtw=np.ones(4),hkil=np.ones(4))
def test_double_Miller_to_Bravais(self):
with pytest.raises(KeyError):
util.Miller_to_Bravais(uvw=np.ones(4),hkl=np.ones(4))
@pytest.mark.parametrize('vector',np.array([
[1,0,0],
[1,1,0],
[1,1,1],
[1,0,-2],
]))
@pytest.mark.parametrize('kw_Miller,kw_Bravais',[('uvw','uvtw'),('hkl','hkil')])
def test_Miller_Bravais_Miller(self,vector,kw_Miller,kw_Bravais):
assert np.all(vector == util.Bravais_to_Miller(**{kw_Bravais:util.Miller_to_Bravais(**{kw_Miller:vector})}))
@pytest.mark.parametrize('vector',np.array([
[1,0,-1,2],
[1,-1,0,3],
[1,1,-2,-3],
[0,0,0,1],
]))
@pytest.mark.parametrize('kw_Miller,kw_Bravais',[('uvw','uvtw'),('hkl','hkil')])
def test_Bravais_Miller_Bravais(self,vector,kw_Miller,kw_Bravais):
assert np.all(vector == util.Miller_to_Bravais(**{kw_Miller:util.Bravais_to_Miller(**{kw_Bravais:vector})}))

7
src/IO.f90
View File

@ -16,7 +16,8 @@ module IO
private private
character(len=*), parameter, public :: & character(len=*), parameter, public :: &
IO_WHITESPACE = achar(44)//achar(32)//achar(9)//achar(10)//achar(13) !< whitespace characters IO_WHITESPACE = achar(44)//achar(32)//achar(9)//achar(10)//achar(13), & !< whitespace characters
IO_QUOTES = "'"//'"'
character, parameter, public :: & character, parameter, public :: &
IO_EOL = new_line('DAMASK'), & !< end of line character IO_EOL = new_line('DAMASK'), & !< end of line character
IO_COMMENT = '#' IO_COMMENT = '#'
@ -495,13 +496,13 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = '--- expected after YAML file header' msg = '--- expected after YAML file header'
case (709) case (709)
msg = 'Length mismatch' msg = 'Length mismatch'
case (710)
msg = 'Closing quotation mark missing in string'
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
! errors related to the grid solver ! errors related to the grid solver
case (831) case (831)
msg = 'mask consistency violated in grid load case' msg = 'mask consistency violated in grid load case'
case (832)
msg = 'ill-defined L (line partly defined) in grid load case'
case (833) case (833)
msg = 'non-positive ratio for geometric progression' msg = 'non-positive ratio for geometric progression'
case (834) case (834)

85
src/YAML_parse.f90
View File

@ -71,8 +71,8 @@ recursive function parse_flow(YAML_flow) result(node)
s = e s = e
d = s + scan(flow_string(s+1:),':') d = s + scan(flow_string(s+1:),':')
e = d + find_end(flow_string(d+1:),'}') e = d + find_end(flow_string(d+1:),'}')
key = trim(adjustl(flow_string(s+1:d-1))) key = trim(adjustl(flow_string(s+1:d-1)))
if(quotedString(key)) key = key(2:len(key)-1)
myVal => parse_flow(flow_string(d+1:e-1)) ! parse items (recursively) myVal => parse_flow(flow_string(d+1:e-1)) ! parse items (recursively)
select type (node) select type (node)
@ -97,7 +97,11 @@ recursive function parse_flow(YAML_flow) result(node)
allocate(tScalar::node) allocate(tScalar::node)
select type (node) select type (node)
class is (tScalar) class is (tScalar)
node = trim(adjustl(flow_string)) if(quotedString(flow_string)) then
node = trim(adjustl(flow_string(2:len(flow_string)-1)))
else
node = trim(adjustl(flow_string))
endif
end select end select
endif endif
@ -119,18 +123,38 @@ integer function find_end(str,e_char)
N_sq = 0 N_sq = 0
N_cu = 0 N_cu = 0
do i = 1, len_trim(str) i = 1
do while(i<=len_trim(str))
if (scan(str(i:i),IO_QUOTES) == 1) i = i + scan(str(i+1:),str(i:i))
if (N_sq==0 .and. N_cu==0 .and. scan(str(i:i),e_char//',') == 1) exit if (N_sq==0 .and. N_cu==0 .and. scan(str(i:i),e_char//',') == 1) exit
N_sq = N_sq + merge(1,0,str(i:i) == '[') N_sq = N_sq + merge(1,0,str(i:i) == '[')
N_cu = N_cu + merge(1,0,str(i:i) == '{') N_cu = N_cu + merge(1,0,str(i:i) == '{')
N_sq = N_sq - merge(1,0,str(i:i) == ']') N_sq = N_sq - merge(1,0,str(i:i) == ']')
N_cu = N_cu - merge(1,0,str(i:i) == '}') N_cu = N_cu - merge(1,0,str(i:i) == '}')
i = i + 1
enddo enddo
find_end = i find_end = i
end function find_end end function find_end
!--------------------------------------------------------------------------------------------------
! @brief check whether a string is enclosed with single or double quotes
!--------------------------------------------------------------------------------------------------
logical function quotedString(line)
character(len=*), intent(in) :: line
quotedString = .false.
if (scan(line(:1),IO_QUOTES) == 1) then
quotedString = .true.
if(line(len(line):len(line)) /= line(:1)) call IO_error(710,ext_msg=line)
endif
end function quotedString
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! @brief Returns Indentation. ! @brief Returns Indentation.
! @details It determines the indentation level for a given block/line. ! @details It determines the indentation level for a given block/line.
@ -333,6 +357,36 @@ subroutine remove_line_break(blck,s_blck,e_char,flow_line)
end subroutine remove_line_break end subroutine remove_line_break
!--------------------------------------------------------------------------------------------------
!> @brief return the scalar list item without line break
!--------------------------------------------------------------------------------------------------
subroutine list_item_inline(blck,s_blck,inline) !ToDo: SR: merge with remove_line_break eventually
character(len=*), intent(in) :: blck !< YAML in mixed style
integer, intent(inout) :: s_blck
character(len=:), allocatable, intent(out) :: inline
character(len=:), allocatable :: line
integer :: indent,indent_next
line = IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))
indent = indentDepth(blck(s_blck:))
inline = line(indent+3:)
s_blck = s_blck + index(blck(s_blck:),IO_EOL)
indent_next = indentDepth(blck(s_blck:))
do while(indent_next > indent)
inline = inline//' '//trim(adjustl(IO_rmComment(blck(s_blck:s_blck + index(blck(s_blck:),IO_EOL) - 2))))
s_blck = s_blck + index(blck(s_blck:),IO_EOL)
indent_next = indentDepth(blck(s_blck:))
enddo
if(scan(inline,",") > 0) inline = '"'//inline//'"'
end subroutine list_item_inline
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! @brief reads a line of YAML block which is already in flow style ! @brief reads a line of YAML block which is already in flow style
! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser ! @details Dicts should be enlcosed within '{}' for it to be consistent with DAMASK YAML parser
@ -463,7 +517,7 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
integer, intent(inout) :: s_blck, & !< start position in blck integer, intent(inout) :: s_blck, & !< start position in blck
s_flow, & !< start position in flow s_flow, & !< start position in flow
offset !< stores leading '- ' in nested lists offset !< stores leading '- ' in nested lists
character(len=:), allocatable :: line,flow_line character(len=:), allocatable :: line,flow_line,inline
integer :: e_blck,indent integer :: e_blck,indent
indent = indentDepth(blck(s_blck:),offset) indent = indentDepth(blck(s_blck:),offset)
@ -509,8 +563,8 @@ recursive subroutine lst(blck,flow,s_blck,s_flow,offset)
else ! list item in the same line else ! list item in the same line
line = line(indentDepth(line)+3:) line = line(indentDepth(line)+3:)
if(isScalar(line)) then if(isScalar(line)) then
call line_toFlow(flow,s_flow,line) call list_item_inline(blck,s_blck,inline)
s_blck = e_blck +2 call line_toFlow(flow,s_flow,inline)
offset = 0 offset = 0
elseif(isFlow(line)) then elseif(isFlow(line)) then
s_blck = s_blck + index(blck(s_blck:),'-') s_blck = s_blck + index(blck(s_blck:),'-')
@ -723,6 +777,8 @@ subroutine selfTest
if (indentDepth('a') /= 0) error stop 'indentDepth' if (indentDepth('a') /= 0) error stop 'indentDepth'
if (indentDepth('x ') /= 0) error stop 'indentDepth' if (indentDepth('x ') /= 0) error stop 'indentDepth'
if (.not. quotedString("'a'")) error stop 'quotedString'
if ( isFlow(' a')) error stop 'isFLow' if ( isFlow(' a')) error stop 'isFLow'
if (.not. isFlow('{')) error stop 'isFlow' if (.not. isFlow('{')) error stop 'isFlow'
if (.not. isFlow(' [')) error stop 'isFlow' if (.not. isFlow(' [')) error stop 'isFlow'
@ -809,14 +865,14 @@ subroutine selfTest
multi_line_flow1: block multi_line_flow1: block
character(len=*), parameter :: flow_multi = & character(len=*), parameter :: flow_multi = &
"%YAML 1.1"//IO_EOL//& '%YAML 1.1'//IO_EOL//&
"---"//IO_EOL//& '---'//IO_EOL//&
"a: [b,"//IO_EOL//& 'a: ["b",'//IO_EOL//&
"c: "//IO_EOL//& 'c: '//IO_EOL//&
"d, e]"//IO_EOL '"d", "e"]'//IO_EOL
character(len=*), parameter :: flow = & character(len=*), parameter :: flow = &
"{a: [b, {c: d}, e]}" '{a: ["b", {c: "d"}, "e"]}'
if( .not. to_flow(flow_multi) == flow) error stop 'to_flow' if( .not. to_flow(flow_multi) == flow) error stop 'to_flow'
end block multi_line_flow1 end block multi_line_flow1
@ -848,14 +904,15 @@ subroutine selfTest
" "//IO_EOL//& " "//IO_EOL//&
" "//IO_EOL//& " "//IO_EOL//&
" param_1: [a: b, c, {d: {e: [f: g, h]}}]"//IO_EOL//& " param_1: [a: b, c, {d: {e: [f: g, h]}}]"//IO_EOL//&
" - c: d"//IO_EOL//& " - c:d"//IO_EOL//&
" e.f,"//IO_EOL//&
" bb:"//IO_EOL//& " bb:"//IO_EOL//&
" "//IO_EOL//& " "//IO_EOL//&
" - "//IO_EOL//& " - "//IO_EOL//&
" {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL//& " {param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}"//IO_EOL//&
"..."//IO_EOL "..."//IO_EOL
character(len=*), parameter :: mixed_flow = & character(len=*), parameter :: mixed_flow = &
"{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, {c: d}], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}" '{aa: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}, "c:d e.f,"], bb: [{param_1: [{a: b}, c, {d: {e: [{f: g}, h]}}]}]}'
if(.not. to_flow(block_flow) == mixed_flow) error stop 'to_flow' if(.not. to_flow(block_flow) == mixed_flow) error stop 'to_flow'
end block basic_mixed end block basic_mixed

73
src/grid/DAMASK_grid.f90
View File

@ -53,12 +53,11 @@ program DAMASK_grid
integer, parameter :: & integer, parameter :: &
subStepFactor = 2 !< for each substep, divide the last time increment by 2.0 subStepFactor = 2 !< for each substep, divide the last time increment by 2.0
real(pReal) :: & real(pReal) :: &
T_0 = 300.0_pReal, & t = 0.0_pReal, & !< elapsed time
time = 0.0_pReal, & !< elapsed time t_0 = 0.0_pReal, & !< begin of interval
time0 = 0.0_pReal, & !< begin of interval Delta_t = 1.0_pReal, & !< current time interval
timeinc = 1.0_pReal, & !< current time interval Delta_t_prev = 0.0_pReal, & !< previous time interval
timeIncOld = 0.0_pReal, & !< previous time interval t_remaining = 0.0_pReal !< remaining time of current load case
remainingLoadCaseTime = 0.0_pReal !< remaining time of current load case
logical :: & logical :: &
guess, & !< guess along former trajectory guess, & !< guess along former trajectory
stagIterate, & stagIterate, &
@ -230,12 +229,6 @@ program DAMASK_grid
reportAndCheck: if (worldrank == 0) then reportAndCheck: if (worldrank == 0) then
print'(/,a,i0)', ' load case: ', l print'(/,a,i0)', ' load case: ', l
print*, ' estimate_rate:', loadCases(l)%estimate_rate print*, ' estimate_rate:', loadCases(l)%estimate_rate
if (loadCases(l)%deformation%myType == 'L') then
do j = 1, 3
if (any(loadCases(l)%deformation%mask(j,1:3) .eqv. .true.) .and. &
any(loadCases(l)%deformation%mask(j,1:3) .eqv. .false.)) errorID = 832 ! each row should be either fully or not at all defined
enddo
endif
if (loadCases(l)%deformation%myType == 'F') then if (loadCases(l)%deformation%myType == 'F') then
print*, ' F:' print*, ' F:'
else else
@ -243,14 +236,14 @@ program DAMASK_grid
endif endif
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if (loadCases(l)%deformation%mask(i,j)) then if (loadCases(l)%deformation%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j)
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' write(IO_STDOUT,'(2x,12a)',advance='no') ' x '
else
write(IO_STDOUT,'(2x,f12.7)',advance='no') loadCases(l)%deformation%values(i,j)
endif endif
enddo; write(IO_STDOUT,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831 if (any(loadCases(l)%stress%mask .eqv. loadCases(l)%deformation%mask)) errorID = 831
if (any(loadCases(l)%stress%mask .and. transpose(loadCases(l)%stress%mask) .and. (math_I3<1))) & if (any(.not.(loadCases(l)%stress%mask .or. transpose(loadCases(l)%stress%mask)) .and. (math_I3<1))) &
errorID = 838 ! no rotation is allowed by stress BC errorID = 838 ! no rotation is allowed by stress BC
if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:' if (loadCases(l)%stress%myType == 'P') print*, ' P / MPa:'
@ -259,9 +252,9 @@ program DAMASK_grid
if (loadCases(l)%stress%myType /= '') then if (loadCases(l)%stress%myType /= '') then
do i = 1, 3; do j = 1, 3 do i = 1, 3; do j = 1, 3
if (loadCases(l)%stress%mask(i,j)) then if (loadCases(l)%stress%mask(i,j)) then
write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal
else
write(IO_STDOUT,'(2x,12a)',advance='no') ' x ' write(IO_STDOUT,'(2x,12a)',advance='no') ' x '
else
write(IO_STDOUT,'(2x,f12.4)',advance='no') loadCases(l)%stress%values(i,j)*1e-6_pReal
endif endif
enddo; write(IO_STDOUT,'(/)',advance='no') enddo; write(IO_STDOUT,'(/)',advance='no')
enddo enddo
@ -303,7 +296,7 @@ program DAMASK_grid
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict) initial_conditions => config_load%get('initial_conditions',defaultVal=emptyDict)
thermal => initial_conditions%get('thermal',defaultVal=emptyDict) thermal => initial_conditions%get('thermal',defaultVal=emptyDict)
call grid_thermal_spectral_init(thermal%get_asFloat('T',defaultVal = T_0)) call grid_thermal_spectral_init(thermal%get_asFloat('T'))
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
call grid_damage_spectral_init call grid_damage_spectral_init
@ -330,7 +323,7 @@ program DAMASK_grid
endif writeUndeformed endif writeUndeformed
loadCaseLooping: do l = 1, size(loadCases) loadCaseLooping: do l = 1, size(loadCases)
time0 = time ! load case start time t_0 = t ! load case start time
guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc guess = loadCases(l)%estimate_rate ! change of load case? homogeneous guess for the first inc
incLooping: do inc = 1, loadCases(l)%N incLooping: do inc = 1, loadCases(l)%N
@ -338,31 +331,31 @@ program DAMASK_grid
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! forwarding time ! forwarding time
timeIncOld = timeinc ! last timeinc that brought former inc to an end Delta_t_prev = Delta_t ! last time intervall that brought former inc to an end
if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then ! linear scale if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then ! linear scale
timeinc = loadCases(l)%t/real(loadCases(l)%N,pReal) Delta_t = loadCases(l)%t/real(loadCases(l)%N,pReal)
else else
timeinc = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) & Delta_t = loadCases(l)%t * (loadCases(l)%r**(inc-1)-loadCases(l)%r**inc) &
/ (1.0_pReal-loadCases(l)%r**loadCases(l)%N) / (1.0_pReal-loadCases(l)%r**loadCases(l)%N)
endif endif
timeinc = timeinc * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step Delta_t = Delta_t * real(subStepFactor,pReal)**real(-cutBackLevel,pReal) ! depending on cut back level, decrease time step
skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc? skipping: if (totalIncsCounter <= interface_restartInc) then ! not yet at restart inc?
time = time + timeinc ! just advance time, skip already performed calculation t = t + Delta_t ! just advance time, skip already performed calculation
guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference guess = .true. ! QUESTION:why forced guessing instead of inheriting loadcase preference
else skipping else skipping
stepFraction = 0 ! fraction scaled by stepFactor**cutLevel stepFraction = 0 ! fraction scaled by stepFactor**cutLevel
subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel) subStepLooping: do while (stepFraction < subStepFactor**cutBackLevel)
remainingLoadCaseTime = loadCases(l)%t+time0 - time t_remaining = loadCases(l)%t + t_0 - t
time = time + timeinc ! forward target time t = t + Delta_t ! forward target time
stepFraction = stepFraction + 1 ! count step stepFraction = stepFraction + 1 ! count step
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! report begin of new step ! report begin of new step
print'(/,a)', ' ###########################################################################' print'(/,a)', ' ###########################################################################'
print'(1x,a,es12.5,6(a,i0))', & print'(1x,a,es12.5,6(a,i0))', &
'Time', time, & 'Time', t, &
's: Increment ', inc,'/',loadCases(l)%N,& 's: Increment ', inc,'/',loadCases(l)%N,&
'-', stepFraction,'/',subStepFactor**cutBackLevel,& '-', stepFraction,'/',subStepFactor**cutBackLevel,&
' of load case ', l,'/',size(loadCases) ' of load case ', l,'/',size(loadCases)
@ -377,10 +370,10 @@ program DAMASK_grid
select case(ID(field)) select case(ID(field))
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
call mechanical_forward (& call mechanical_forward (&
cutBack,guess,timeinc,timeIncOld,remainingLoadCaseTime, & cutBack,guess,Delta_t,Delta_t_prev,t_remaining, &
deformation_BC = loadCases(l)%deformation, & deformation_BC = loadCases(l)%deformation, &
stress_BC = loadCases(l)%stress, & stress_BC = loadCases(l)%stress, &
rotation_BC = loadCases(l)%rot) rotation_BC = loadCases(l)%rot)
case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack) case(FIELD_THERMAL_ID); call grid_thermal_spectral_forward(cutBack)
case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack) case(FIELD_DAMAGE_ID); call grid_damage_spectral_forward(cutBack)
@ -398,9 +391,9 @@ program DAMASK_grid
case(FIELD_MECH_ID) case(FIELD_MECH_ID)
solres(field) = mechanical_solution(incInfo) solres(field) = mechanical_solution(incInfo)
case(FIELD_THERMAL_ID) case(FIELD_THERMAL_ID)
solres(field) = grid_thermal_spectral_solution(timeinc) solres(field) = grid_thermal_spectral_solution(Delta_t)
case(FIELD_DAMAGE_ID) case(FIELD_DAMAGE_ID)
solres(field) = grid_damage_spectral_solution(timeinc) solres(field) = grid_damage_spectral_solution(Delta_t)
end select end select
if (.not. solres(field)%converged) exit ! no solution found if (.not. solres(field)%converged) exit ! no solution found
@ -418,11 +411,11 @@ program DAMASK_grid
if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged if ( (all(solres(:)%converged .and. solres(:)%stagConverged)) & ! converged
.and. .not. solres(1)%termIll) then ! and acceptable solution found .and. .not. solres(1)%termIll) then ! and acceptable solution found
call mechanical_updateCoords call mechanical_updateCoords
timeIncOld = timeinc Delta_t_prev = Delta_t
cutBack = .false. cutBack = .false.
guess = .true. ! start guessing after first converged (sub)inc guess = .true. ! start guessing after first converged (sub)inc
if (worldrank == 0) then if (worldrank == 0) then
write(statUnit,*) totalIncsCounter, time, cutBackLevel, & write(statUnit,*) totalIncsCounter, t, cutBackLevel, &
solres(1)%converged, solres(1)%iterationsNeeded solres(1)%converged, solres(1)%iterationsNeeded
flush(statUnit) flush(statUnit)
endif endif
@ -430,8 +423,8 @@ program DAMASK_grid
cutBack = .true. cutBack = .true.
stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator stepFraction = (stepFraction - 1) * subStepFactor ! adjust to new denominator
cutBackLevel = cutBackLevel + 1 cutBackLevel = cutBackLevel + 1
time = time - timeinc ! rewind time t = t - Delta_t
timeinc = timeinc/real(subStepFactor,pReal) ! cut timestep Delta_t = Delta_t/real(subStepFactor,pReal) ! cut timestep
print'(/,a)', ' cutting back ' print'(/,a)', ' cutting back '
else ! no more options to continue else ! no more options to continue
if (worldrank == 0) close(statUnit) if (worldrank == 0) close(statUnit)
@ -453,7 +446,7 @@ program DAMASK_grid
if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then if (mod(inc,loadCases(l)%f_out) == 0 .or. signal) then
print'(1/,a)', ' ... writing results to file ......................................' print'(1/,a)', ' ... writing results to file ......................................'
flush(IO_STDOUT) flush(IO_STDOUT)
call CPFEM_results(totalIncsCounter,time) call CPFEM_results(totalIncsCounter,t)
endif endif
if(signal) call interface_setSIGUSR1(.false.) if(signal) call interface_setSIGUSR1(.false.)
call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(interface_SIGUSR2,signal,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
@ -497,8 +490,8 @@ subroutine getMaskedTensor(values,mask,tensor)
do i = 1,3 do i = 1,3
row => tensor%get(i) row => tensor%get(i)
do j = 1,3 do j = 1,3
mask(i,j) = row%get_asString(j) /= 'x' mask(i,j) = row%get_asString(j) == 'x'
if (mask(i,j)) values(i,j) = row%get_asFloat(j) if (.not. mask(i,j)) values(i,j) = row%get_asFloat(j)
enddo enddo
enddo enddo

10
src/grid/grid_damage_spectral.f90
View File

@ -160,10 +160,10 @@ end subroutine grid_damage_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral damage scheme with internal iterations !> @brief solution for the spectral damage scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_damage_spectral_solution(timeinc) result(solution) function grid_damage_spectral_solution(Delta_t) result(solution)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc !< increment in time for current solution Delta_t !< increment in time for current solution
integer :: i, j, k, ce integer :: i, j, k, ce
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
@ -176,7 +176,7 @@ function grid_damage_spectral_solution(timeinc) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%Delta_t = Delta_t
call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESSolve(damage_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(damage_snes,reason,ierr); CHKERRQ(ierr)
@ -284,7 +284,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) & scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_phi(phi_current(i,j,k),ce)) &
+ homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) & + homogenization_mu_phi(ce)*(phi_lastInc(i,j,k) - phi_current(i,j,k)) &
+ mu_ref*phi_current(i,j,k) + mu_ref*phi_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -292,7 +292,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of damage field with green operator ! convolution of damage field with green operator
call utilities_FFTscalarForward call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc) call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t)
call utilities_FFTscalarBackward call utilities_FFTscalarBackward
where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) & where(scalarField_real(1:grid(1),1:grid(2),1:grid3) > phi_lastInc) &

20
src/grid/grid_mech_FEM.f90
View File

@ -324,7 +324,7 @@ function grid_mechanical_FEM_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_FEM_solution end function grid_mechanical_FEM_solution
@ -363,20 +363,20 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
else else
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + matmul(merge(.0_pReal,deformation_BC%values,deformation_BC%mask),F_aim_lastInc)
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + merge(.0_pReal,deformation_BC%values,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
endif endif
if (guess) then if (guess) then
@ -397,9 +397,9 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if (stress_BC%myType=='P') P_aim = P_aim & if (stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t
if (stress_BC%myType=='dot_P') P_aim = P_aim & if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr) call VecAXPY(solution_current,Delta_t,solution_rate,ierr); CHKERRQ(ierr)
@ -412,7 +412,7 @@ subroutine grid_mechanical_FEM_forward(cutBack,guess,Delta_t,Delta_t_old,t_remai
! set module wide available data ! set module wide available data
params%stress_mask = stress_BC%mask params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%Delta_t = Delta_t
end subroutine grid_mechanical_FEM_forward end subroutine grid_mechanical_FEM_forward
@ -568,13 +568,13 @@ subroutine formResidual(da_local,x_local, &
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(P_current,& call utilities_constitutiveResponse(P_current,&
P_av,C_volAvg,devNull, & P_av,C_volAvg,devNull, &
F,params%timeinc,params%rotation_BC) F,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask)))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! constructing residual ! constructing residual

20
src/grid/grid_mech_spectral_basic.f90
View File

@ -277,7 +277,7 @@ function grid_mechanical_spectral_basic_solution(incInfoIn) result(solution)
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_basic_solution end function grid_mechanical_spectral_basic_solution
@ -314,20 +314,20 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + matmul(merge(.0_pReal,deformation_BC%values,deformation_BC%mask),F_aim_lastInc)
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + merge(.0_pReal,deformation_BC%values,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
endif endif
Fdot = utilities_calculateRate(guess, & Fdot = utilities_calculateRate(guess, &
@ -342,9 +342,9 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if (stress_BC%myType=='P') P_aim = P_aim & if (stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t
if (stress_BC%myType=='dot_P') P_aim = P_aim & if (stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3]) rotation_BC%rotate(F_aim,active=.true.)),[9,grid(1),grid(2),grid3])
@ -354,7 +354,7 @@ subroutine grid_mechanical_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_
! set module wide available data ! set module wide available data
params%stress_mask = stress_BC%mask params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%Delta_t = Delta_t
end subroutine grid_mechanical_spectral_basic_forward end subroutine grid_mechanical_spectral_basic_forward
@ -492,14 +492,14 @@ subroutine formResidual(in, F, &
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory) call utilities_constitutiveResponse(residuum, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, & P_av,C_volAvg,C_minMaxAvg, &
F,params%timeinc,params%rotation_BC) F,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc deltaF_aim = math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
F_aim = F_aim - deltaF_aim F_aim = F_aim - deltaF_aim
err_BC = maxval(abs(merge(P_av - P_aim,.0_pReal,params%stress_mask))) err_BC = maxval(abs(merge(.0_pReal,P_av - P_aim,params%stress_mask)))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! updated deformation gradient using fix point algorithm of basic scheme ! updated deformation gradient using fix point algorithm of basic scheme

22
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -309,7 +309,7 @@ function grid_mechanical_spectral_polarisation_solution(incInfoIn) result(soluti
solution%iterationsNeeded = totalIter solution%iterationsNeeded = totalIter
solution%termIll = terminallyIll solution%termIll = terminallyIll
terminallyIll = .false. terminallyIll = .false.
P_aim = merge(P_aim,P_av,params%stress_mask) P_aim = merge(P_av,P_aim,params%stress_mask)
end function grid_mechanical_spectral_polarisation_solution end function grid_mechanical_spectral_polarisation_solution
@ -350,20 +350,20 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
C_volAvgLastInc = C_volAvg C_volAvgLastInc = C_volAvg
C_minMaxAvgLastInc = C_minMaxAvg C_minMaxAvgLastInc = C_minMaxAvg
F_aimDot = merge(merge((F_aim-F_aim_lastInc)/Delta_t_old,0.0_pReal,stress_BC%mask), 0.0_pReal, guess) ! estimate deformation rate for prescribed stress components F_aimDot = merge(merge(.0_pReal,(F_aim-F_aim_lastInc)/Delta_t_old,stress_BC%mask),.0_pReal,guess) ! estimate deformation rate for prescribed stress components
F_aim_lastInc = F_aim F_aim_lastInc = F_aim
!----------------------------------------------------------------------------------------------- !-----------------------------------------------------------------------------------------------
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='L') then ! calculate F_aimDot from given L and current F
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(matmul(deformation_BC%values, F_aim_lastInc),.0_pReal,deformation_BC%mask) + matmul(merge(.0_pReal,deformation_BC%values,deformation_BC%mask),F_aim_lastInc)
elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed elseif (deformation_BC%myType=='dot_F') then ! F_aimDot is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge(deformation_BC%values,.0_pReal,deformation_BC%mask) + merge(.0_pReal,deformation_BC%values,deformation_BC%mask)
elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed elseif (deformation_BC%myType=='F') then ! aim at end of load case is prescribed
F_aimDot = F_aimDot & F_aimDot = F_aimDot &
+ merge((deformation_BC%values - F_aim_lastInc)/t_remaining,.0_pReal,deformation_BC%mask) + merge(.0_pReal,(deformation_BC%values - F_aim_lastInc)/t_remaining,deformation_BC%mask)
endif endif
Fdot = utilities_calculateRate(guess, & Fdot = utilities_calculateRate(guess, &
@ -382,9 +382,9 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
! update average and local deformation gradients ! update average and local deformation gradients
F_aim = F_aim_lastInc + F_aimDot * Delta_t F_aim = F_aim_lastInc + F_aimDot * Delta_t
if(stress_BC%myType=='P') P_aim = P_aim & if(stress_BC%myType=='P') P_aim = P_aim &
+ merge((stress_BC%values - P_aim)/t_remaining,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,(stress_BC%values - P_aim)/t_remaining,stress_BC%mask)*Delta_t
if(stress_BC%myType=='dot_P') P_aim = P_aim & if(stress_BC%myType=='dot_P') P_aim = P_aim &
+ merge(stress_BC%values,0.0_pReal,stress_BC%mask)*Delta_t + merge(.0_pReal,stress_BC%values,stress_BC%mask)*Delta_t
F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average F = reshape(utilities_forwardField(Delta_t,F_lastInc,Fdot, & ! estimate of F at end of time+Delta_t that matches rotated F_aim on average
rotation_BC%rotate(F_aim,active=.true.)),& rotation_BC%rotate(F_aim,active=.true.)),&
@ -408,7 +408,7 @@ subroutine grid_mechanical_spectral_polarisation_forward(cutBack,guess,Delta_t,D
! set module wide available data ! set module wide available data
params%stress_mask = stress_BC%mask params%stress_mask = stress_BC%mask
params%rotation_BC = rotation_BC params%rotation_BC = rotation_BC
params%timeinc = Delta_t params%Delta_t = Delta_t
end subroutine grid_mechanical_spectral_polarisation_forward end subroutine grid_mechanical_spectral_polarisation_forward
@ -592,14 +592,14 @@ subroutine formResidual(in, FandF_tau, &
! evaluate constitutive response ! evaluate constitutive response
call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory) call utilities_constitutiveResponse(residual_F, & ! "residuum" gets field of first PK stress (to save memory)
P_av,C_volAvg,C_minMaxAvg, & P_av,C_volAvg,C_minMaxAvg, &
F - residual_F_tau/num%beta,params%timeinc,params%rotation_BC) F - residual_F_tau/num%beta,params%Delta_t,params%rotation_BC)
call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,terminallyIll,1,MPI_LOGICAL,MPI_LOR,MPI_COMM_WORLD,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress BC handling ! stress BC handling
F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc F_aim = F_aim - math_mul3333xx33(S, P_av - P_aim) ! S = 0.0 for no bc
err_BC = maxval(abs(merge(P_av-P_aim, & err_BC = maxval(abs(merge(math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)), &
math_mul3333xx33(C_scale,F_aim-params%rotation_BC%rotate(F_av)),& P_av-P_aim, &
params%stress_mask))) params%stress_mask)))
! calculate divergence ! calculate divergence
tensorField_real = 0.0_pReal tensorField_real = 0.0_pReal

14
src/grid/grid_thermal_spectral.f90
View File

@ -155,10 +155,10 @@ end subroutine grid_thermal_spectral_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief solution for the spectral thermal scheme with internal iterations !> @brief solution for the spectral thermal scheme with internal iterations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function grid_thermal_spectral_solution(timeinc) result(solution) function grid_thermal_spectral_solution(Delta_t) result(solution)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc !< increment in time for current solution Delta_t !< increment in time for current solution
integer :: i, j, k, ce integer :: i, j, k, ce
type(tSolutionState) :: solution type(tSolutionState) :: solution
PetscInt :: devNull PetscInt :: devNull
@ -171,7 +171,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! set module wide availabe data ! set module wide availabe data
params%timeinc = timeinc params%Delta_t = Delta_t
call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr) call SNESSolve(thermal_snes,PETSC_NULL_VEC,solution_vec,ierr); CHKERRQ(ierr)
call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr) call SNESGetConvergedReason(thermal_snes,reason,ierr); CHKERRQ(ierr)
@ -195,7 +195,7 @@ function grid_thermal_spectral_solution(timeinc) result(solution)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce) call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
enddo; enddo; enddo enddo; enddo; enddo
call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr) call VecMin(solution_vec,devNull,T_min,ierr); CHKERRQ(ierr)
@ -233,7 +233,7 @@ subroutine grid_thermal_spectral_forward(cutBack)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%timeinc,ce) call homogenization_thermal_setField(T_current(i,j,k),(T_current(i,j,k)-T_lastInc(i,j,k))/params%Delta_t,ce)
enddo; enddo; enddo enddo; enddo; enddo
else else
T_lastInc = T_current T_lastInc = T_current
@ -279,7 +279,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
ce = 0 ce = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
ce = ce + 1 ce = ce + 1
scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + homogenization_f_T(ce)) & scalarField_real(i,j,k) = params%Delta_t*(scalarField_real(i,j,k) + homogenization_f_T(ce)) &
+ homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) & + homogenization_mu_T(ce) * (T_lastInc(i,j,k) - T_current(i,j,k)) &
+ mu_ref*T_current(i,j,k) + mu_ref*T_current(i,j,k)
enddo; enddo; enddo enddo; enddo; enddo
@ -287,7 +287,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! convolution of temperature field with green operator ! convolution of temperature field with green operator
call utilities_FFTscalarForward call utilities_FFTscalarForward
call utilities_fourierGreenConvolution(K_ref, mu_ref, params%timeinc) call utilities_fourierGreenConvolution(K_ref, mu_ref, params%Delta_t)
call utilities_FFTscalarBackward call utilities_FFTscalarBackward
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

28
src/grid/spectral_utilities.f90
View File

@ -88,7 +88,7 @@ module spectral_utilities
type, public :: tBoundaryCondition !< set of parameters defining a boundary condition type, public :: tBoundaryCondition !< set of parameters defining a boundary condition
real(pReal), dimension(3,3) :: values = 0.0_pReal real(pReal), dimension(3,3) :: values = 0.0_pReal
logical, dimension(3,3) :: mask = .false. logical, dimension(3,3) :: mask = .true.
character(len=:), allocatable :: myType character(len=:), allocatable :: myType
end type tBoundaryCondition end type tBoundaryCondition
@ -96,7 +96,7 @@ module spectral_utilities
real(pReal), dimension(3,3) :: stress_BC real(pReal), dimension(3,3) :: stress_BC
logical, dimension(3,3) :: stress_mask logical, dimension(3,3) :: stress_mask
type(rotation) :: rotation_BC type(rotation) :: rotation_BC
real(pReal) :: timeinc real(pReal) :: Delta_t
end type tSolutionParams end type tSolutionParams
type :: tNumerics type :: tNumerics
@ -684,7 +684,7 @@ function utilities_maskedCompliance(rot_BC,mask_stress,C)
logical :: errmatinv logical :: errmatinv
character(len=pStringLen):: formatString character(len=pStringLen):: formatString
mask_stressVector = reshape(transpose(mask_stress), [9]) mask_stressVector = .not. reshape(transpose(mask_stress), [9])
size_reduced = count(mask_stressVector) size_reduced = count(mask_stressVector)
if(size_reduced > 0) then if(size_reduced > 0) then
temp99_real = math_3333to99(rot_BC%rotate(C)) temp99_real = math_3333to99(rot_BC%rotate(C))
@ -791,16 +791,16 @@ end subroutine utilities_fourierTensorDivergence
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate constitutive response from homogenization_F0 to F during timeinc !> @brief calculate constitutive response from homogenization_F0 to F during Delta_t
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
F,timeinc,rotation_BC) F,Delta_t,rotation_BC)
real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness real(pReal), intent(out), dimension(3,3,3,3) :: C_volAvg, C_minmaxAvg !< average stiffness
real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress real(pReal), intent(out), dimension(3,3) :: P_av !< average PK stress
real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress real(pReal), intent(out), dimension(3,3,grid(1),grid(2),grid3) :: P !< PK stress
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: F !< deformation gradient target
real(pReal), intent(in) :: timeinc !< loading time real(pReal), intent(in) :: Delta_t !< loading time
type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame type(rotation), intent(in), optional :: rotation_BC !< rotation of load frame
@ -815,11 +815,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call homogenization_mechanical_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field call homogenization_mechanical_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3]) ! calculate P field
if (.not. terminallyIll) & if (.not. terminallyIll) &
call homogenization_thermal_response(timeinc,[1,1],[1,product(grid(1:2))*grid3]) call homogenization_thermal_response(Delta_t,[1,1],[1,product(grid(1:2))*grid3])
if (.not. terminallyIll) & if (.not. terminallyIll) &
call homogenization_mechanical_response2(timeinc,[1,1],[1,product(grid(1:2))*grid3]) call homogenization_mechanical_response2(Delta_t,[1,1],[1,product(grid(1:2))*grid3])
P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3]) P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
@ -870,14 +870,14 @@ end subroutine utilities_constitutiveResponse
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates forward rate, either guessing or just add delta/timeinc !> @brief calculates forward rate, either guessing or just add delta/Delta_t
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate) pure function utilities_calculateRate(heterogeneous,field0,field,dt,avRate)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
avRate !< homogeneous addon avRate !< homogeneous addon
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
dt !< timeinc between field0 and field dt !< Delta_t between field0 and field
logical, intent(in) :: & logical, intent(in) :: &
heterogeneous !< calculate field of rates heterogeneous !< calculate field of rates
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
@ -899,10 +899,10 @@ end function utilities_calculateRate
!> @brief forwards a field with a pointwise given rate, if aim is given, !> @brief forwards a field with a pointwise given rate, if aim is given,
!> ensures that the average matches the aim !> ensures that the average matches the aim
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function utilities_forwardField(timeinc,field_lastInc,rate,aim) function utilities_forwardField(Delta_t,field_lastInc,rate,aim)
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
timeinc !< timeinc of current step Delta_t !< Delta_t of current step
real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: & real(pReal), intent(in), dimension(3,3,grid(1),grid(2),grid3) :: &
field_lastInc, & !< initial field field_lastInc, & !< initial field
rate !< rate by which to forward rate !< rate by which to forward
@ -913,7 +913,7 @@ function utilities_forwardField(timeinc,field_lastInc,rate,aim)
real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim real(pReal), dimension(3,3) :: fieldDiff !< <a + adot*t> - aim
PetscErrorCode :: ierr PetscErrorCode :: ierr
utilities_forwardField = field_lastInc + rate*timeinc utilities_forwardField = field_lastInc + rate*Delta_t
if (present(aim)) then !< correct to match average if (present(aim)) then !< correct to match average
fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt fieldDiff = sum(sum(sum(utilities_forwardField,dim=5),dim=4),dim=3)*wgt
call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,fieldDiff,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)