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 python-vtk-improvements

This commit is contained in:
Philip Eisenlohr 2022-02-18 14:46:40 -05:00
parent 75f535b60c d8fab8cb07
commit 0adc827278
23 changed files with 732 additions and 638 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
.gitlab-ci.yml
View File

@ -96,14 +96,14 @@ mesh_GNU:
grid_GNU-64bit: grid_GNU-64bit:
stage: compile stage: compile
script: script:
- module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit - module load Compiler/GNU/10 Libraries/PETSc/3.16.4/64bit
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit - pytest -k 'compile and grid' --basetemp ${TESTROOT}/compile_grid_GNU-64bit
mesh_GNU-64bit: mesh_GNU-64bit:
stage: compile stage: compile
script: script:
- module load Compiler/GNU/10 Libraries/PETSc/3.16.2/64bit - module load Compiler/GNU/10 Libraries/PETSc/3.16.4/64bit
- cd PRIVATE/testing/pytest - cd PRIVATE/testing/pytest
- pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU-64bit - pytest -k 'compile and mesh' --basetemp ${TESTROOT}/compile_mesh_GNU-64bit

2
PRIVATE

@ -1 +1 @@
Subproject commit 0d639a9ba41db279b0d2825c8e8eddf0ccd91326 Subproject commit e663a548b10ee3f9ced35c5a5105bd6824d3eebf

73
examples/config/Phase_Dislotwin_TWIP-Steel-FeMnC.yaml
View File

@ -1,31 +1,62 @@
type: dislotwin type: dislotwin
output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, tau_hat_tw, f_tr] output: [rho_mob, rho_dip, gamma_sl, Lambda_sl, tau_pass, f_tw, Lambda_tw, f_tr]
D: 2.0e-5
# Slip
N_sl: [12] N_sl: [12]
b_sl: [2.56e-10] b_sl: [2.56e-10]
Q_sl: [3.5e-19]
p_sl: [0.325]
q_sl: [1.55]
B: [0.001]
i_sl: [30.0]
rho_mob_0: [1.0e+12] rho_mob_0: [1.0e+12]
rho_dip_0: [1.0] rho_dip_0: [1.0]
v_0: [1.0e+4] v_0: [1.0e+4]
Q_sl: [3.7e-19] tau_0: [3.5e+8]
p_sl: [1.0]
q_sl: [1.0]
tau_0: [1.5e+8]
i_sl: [10.0] # Adj. parameter controlling dislocation mean free path
D_0: 4.0e-5 # Vacancy diffusion prefactor / m^2/s
D_a: 1.0 # minimum dipole distance / b D_a: 1.0 # minimum dipole distance / b
Q_cl: 4.5e-19 # Activation energy for climb / J Q_cl: 4.5e-19 # Activation energy for climb / J
h_sl-sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.137, 0.122] # Interaction coefficients (Kubin et al. 2008) h_sl-sl: [0.122, 0.122, 0.625, 0.07, 0.137, 0.137, 0.122] # Interaction coefficients (Kubin et al. 2008)
# twinning parameters h_tw-sl: [0.0] # ToDo: values not known
# Twin
N_tw: [12] N_tw: [12]
b_tw: [1.47e-10] # Burgers vector length of twin system / b b_tw: [1.40e-10]
t_tw: [5.0e-8] # Twin stack mean thickness / m L_tw: 1.8e-7
L_tw: 442.0 # Length of twin nuclei / b i_tw: 10.0
x_c_tw: 1.0e-9 # critical distance for formation of twin nucleus / m t_tw: [5.0e-8]
V_cs: 1.67e-29 # cross slip volume / m^3 p_tw: [3.0] # A
p_tw: [10.0] # r-exponent in twin formation probability
i_tw: 1.0 # Adj. parameter controlling twin mean free path h_tw-tw: [0.0, 1.0] # ToDo: values not known
h_sl-tw: [0.0, 1.0, 1.0] # dislocation-twin interaction coefficients h_sl-tw: [0.0, 1.0, 1.0] # ToDo: values not known
h_tw-tw: [0.0, 1.0] # twin-twin interaction coefficients
T_ref: 0.0 # Transformation
Gamma_sf: -0.0396 # stacking fault energy / J/m^2 at zero K; TWIP steel: -0.0526; Cu: -0.0396 N_tr: [12]
Gamma_sf,T: 0.0002 # temperature dependence / J/(m^2 K) of stacking fault energy b_tr: [1.40e-10]
L_tr: 1.8e-7
i_tr: 3.0
t_tr: [1.0e-7]
p_tr: [3.0] # B
V_mol: 7.09e-6
Delta_G: 1.33343932e+02
Delta_G,T: 2.56640412
Delta_G,T^2: 1.49524179e-03
h_tr-tr: [1.0, 1.0] # guessing
h_sl-tr: [0.0, 1.0, 1.0] # guessing
# Twin & Transformation
T_ref: 298.15
Gamma_sf: 2.89394017e-02
Gamma_sf,T: 1.22823814e-04
Gamma_sf,T^2: 1.47322968e-07
x_c: 1.0e-9
V_cs: 1.67e-29
a_cF: 3.62e-10
c/a_hP: 1.633
# Slip & Twin & Transformation
D: 2.0e-5

2
python/damask/VERSION
View File

@ -1 +1 @@
v3.0.0-alpha5-651-gd4f711416 v3.0.0-alpha6

13
python/damask/_config.py
View File

@ -144,10 +144,9 @@ class Config(dict):
Configuration from file. Configuration from file.
""" """
if isinstance(fname, (str, Path)): fhandle = open(fname) if isinstance(fname, (str, Path)) else \
fhandle = open(fname) fname
else:
fhandle = fname
return cls(yaml.safe_load(fhandle)) return cls(yaml.safe_load(fhandle))
def save(self, def save(self,
@ -164,10 +163,8 @@ class Config(dict):
Keyword arguments parsed to yaml.dump. Keyword arguments parsed to yaml.dump.
""" """
if isinstance(fname, (str, Path)): fhandle = open(fname,'w',newline='\n') if isinstance(fname, (str, Path)) else \
fhandle = open(fname,'w',newline='\n') fname
else:
fhandle = fname
if 'width' not in kwargs: if 'width' not in kwargs:
kwargs['width'] = 256 kwargs['width'] = 256

54
python/damask/_crystal.py
View File

@ -2,10 +2,11 @@ from typing import Union, Dict, List, Tuple
import numpy as np import numpy as np
from ._typehints import FloatSequence, CrystalFamily, CrystalLattice, CrystalKinematics
from . import util from . import util
from . import Rotation from . import Rotation
lattice_symmetries = { lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
'aP': 'triclinic', 'aP': 'triclinic',
'mP': 'monoclinic', 'mP': 'monoclinic',
@ -30,9 +31,9 @@ lattice_symmetries = {
class Crystal(): class Crystal():
"""Crystal lattice.""" """Crystal lattice."""
def __init__(self,*, def __init__(self, *,
family = None, family: CrystalFamily = None,
lattice = None, lattice: CrystalLattice = None,
a: float = None, b: float = None, c: float = None, a: float = None, b: float = None, c: float = None,
alpha: float = None, beta: float = None, gamma: float = None, alpha: float = None, beta: float = None, gamma: float = None,
degrees: bool = False): degrees: bool = False):
@ -130,9 +131,8 @@ class Crystal():
Crystal to check for equality. Crystal to check for equality.
""" """
if not isinstance(other, Crystal): return NotImplemented if not isinstance(other, Crystal) else \
return NotImplemented self.lattice == other.lattice and \
return self.lattice == other.lattice and \
self.parameters == other.parameters and \ self.parameters == other.parameters and \
self.family == other.family self.family == other.family
@ -208,7 +208,7 @@ class Crystal():
... } ... }
""" """
_basis = { _basis: Dict[CrystalFamily, Dict[str, np.ndarray]] = {
'cubic': {'improper':np.array([ [-1. , 0. , 1. ], 'cubic': {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ], [ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]), [ 0. , np.sqrt(3.) , 0. ] ]),
@ -315,58 +315,58 @@ class Crystal():
self.lattice[-1],None),dtype=float) self.lattice[-1],None),dtype=float)
def to_lattice(self, *, def to_lattice(self, *,
direction: np.ndarray = None, direction: FloatSequence = None,
plane: np.ndarray = None) -> np.ndarray: plane: FloatSequence = None) -> np.ndarray:
""" """
Calculate lattice vector corresponding to crystal frame direction or plane normal. Calculate lattice vector corresponding to crystal frame direction or plane normal.
Parameters Parameters
---------- ----------
direction|plane : numpy.ndarray of shape (...,3) direction|plane : numpy.ndarray, shape (...,3)
Vector along direction or plane normal. Vector along direction or plane normal.
Returns Returns
------- -------
Miller : numpy.ndarray of shape (...,3) Miller : numpy.ndarray, shape (...,3)
Lattice vector of direction or plane. Lattice vector of direction or plane.
Use util.scale_to_coprime to convert to (integer) Miller indices. Use util.scale_to_coprime to convert to (integer) Miller indices.
""" """
if (direction is not None) ^ (plane is None): if (direction is not None) ^ (plane is None):
raise KeyError('specify either "direction" or "plane"') raise KeyError('specify either "direction" or "plane"')
axis,basis = (np.array(direction),self.basis_reciprocal.T) \ basis,axis = (self.basis_reciprocal,np.array(direction)) \
if plane is None else \ if plane is None else \
(np.array(plane),self.basis_real.T) (self.basis_real,np.array(plane))
return np.einsum('il,...l',basis,axis) return np.einsum('li,...l',basis,axis)
def to_frame(self, *, def to_frame(self, *,
uvw: np.ndarray = None, uvw: FloatSequence = None,
hkl: np.ndarray = None) -> np.ndarray: hkl: FloatSequence = None) -> np.ndarray:
""" """
Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl). Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
Parameters Parameters
---------- ----------
uvw|hkl : numpy.ndarray of shape (...,3) uvw|hkl : numpy.ndarray, shape (...,3)
Miller indices of crystallographic direction or plane normal. Miller indices of crystallographic direction or plane normal.
Returns Returns
------- -------
vector : numpy.ndarray of shape (...,3) vector : numpy.ndarray, shape (...,3)
Crystal frame vector along [uvw] direction or (hkl) plane normal. Crystal frame vector along [uvw] direction or (hkl) plane normal.
""" """
if (uvw is not None) ^ (hkl is None): if (uvw is not None) ^ (hkl is None):
raise KeyError('specify either "uvw" or "hkl"') raise KeyError('specify either "uvw" or "hkl"')
axis,basis = (np.array(uvw),self.basis_real) \ basis,axis = (self.basis_real,np.array(uvw)) \
if hkl is None else \ if hkl is None else \
(np.array(hkl),self.basis_reciprocal) (self.basis_reciprocal,np.array(hkl))
return np.einsum('il,...l',basis,axis) return np.einsum('il,...l',basis,axis)
def kinematics(self, def kinematics(self,
mode: str) -> Dict[str, List[np.ndarray]]: mode: CrystalKinematics) -> Dict[str, List[np.ndarray]]:
""" """
Return crystal kinematics systems. Return crystal kinematics systems.
@ -381,7 +381,7 @@ class Crystal():
Directions and planes of deformation mode families. Directions and planes of deformation mode families.
""" """
_kinematics = { _kinematics: Dict[CrystalLattice, Dict[CrystalKinematics, List[np.ndarray]]] = {
'cF': { 'cF': {
'slip': [np.array([ 'slip': [np.array([
[+0,+1,-1, +1,+1,+1], [+0,+1,-1, +1,+1,+1],
@ -626,7 +626,7 @@ class Crystal():
def relation_operations(self, def relation_operations(self,
model: str) -> Tuple[str, Rotation]: model: str) -> Tuple[CrystalLattice, Rotation]:
""" """
Crystallographic orientation relationships for phase transformations. Crystallographic orientation relationships for phase transformations.
@ -658,7 +658,7 @@ class Crystal():
https://doi.org/10.1016/j.actamat.2004.11.021 https://doi.org/10.1016/j.actamat.2004.11.021
""" """
_orientation_relationships = { _orientation_relationships: Dict[str, Dict[CrystalLattice,np.ndarray]] = {
'KS': { 'KS': {
'cF' : np.array([ 'cF' : np.array([
[[-1, 0, 1],[ 1, 1, 1]], [[-1, 0, 1],[ 1, 1, 1]],

34
python/damask/_grid.py
View File

@ -63,8 +63,8 @@ class Grid:
mat_N = self.N_materials mat_N = self.N_materials
return util.srepr([ return util.srepr([
f'cells: {util.srepr(self.cells, " × ")}', f'cells: {util.srepr(self.cells, " × ")}',
f'size: {util.srepr(self.size, " × ")} / m³', f'size: {util.srepr(self.size, " × ")} m³',
f'origin: {util.srepr(self.origin," ")} / m', f'origin: {util.srepr(self.origin," ")} m',
f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else f'# materials: {mat_N}' + ('' if mat_min == 0 and mat_max+1 == mat_N else
f' (min: {mat_min}, max: {mat_max})') f' (min: {mat_min}, max: {mat_max})')
]) ])
@ -182,7 +182,7 @@ class Grid:
Grid-based geometry from file. Grid-based geometry from file.
""" """
v = VTK.load(fname if str(fname).endswith(('.vti','.vtr')) else str(fname)+'.vti') # compatibility hack v = VTK.load(fname if str(fname).endswith('.vti') else str(fname)+'.vti')
cells = np.array(v.vtk_data.GetDimensions())-1 cells = np.array(v.vtk_data.GetDimensions())-1
bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T bbox = np.array(v.vtk_data.GetBounds()).reshape(3,2).T
comments = v.comments comments = v.comments
@ -603,8 +603,8 @@ class Grid:
>>> import damask >>> import damask
>>> damask.Grid.from_minimal_surface([64]*3,np.ones(3)*1.e-4,'Gyroid') >>> damask.Grid.from_minimal_surface([64]*3,np.ones(3)*1.e-4,'Gyroid')
cells : 64 x 64 x 64 cells : 64 x 64 x 64
size : 0.0001 x 0.0001 x 0.0001 / m³ size : 0.0001 x 0.0001 x 0.0001 m³
origin: 0.0 0.0 0.0 / m origin: 0.0 0.0 0.0 m
# materials: 2 # materials: 2
Minimal surface of 'Neovius' type. non-default material IDs. Minimal surface of 'Neovius' type. non-default material IDs.
@ -614,8 +614,8 @@ class Grid:
>>> damask.Grid.from_minimal_surface([80]*3,np.ones(3)*5.e-4, >>> damask.Grid.from_minimal_surface([80]*3,np.ones(3)*5.e-4,
... 'Neovius',materials=(1,5)) ... 'Neovius',materials=(1,5))
cells : 80 x 80 x 80 cells : 80 x 80 x 80
size : 0.0005 x 0.0005 x 0.0005 / m³ size : 0.0005 x 0.0005 x 0.0005 m³
origin: 0.0 0.0 0.0 / m origin: 0.0 0.0 0.0 m
# materials: 2 (min: 1, max: 5) # materials: 2 (min: 1, max: 5)
""" """
@ -735,8 +735,8 @@ class Grid:
>>> g = damask.Grid(np.zeros([64]*3,int), np.ones(3)*1e-4) >>> g = damask.Grid(np.zeros([64]*3,int), np.ones(3)*1e-4)
>>> g.add_primitive(np.ones(3)*5e-5,np.ones(3)*5e-5,1) >>> g.add_primitive(np.ones(3)*5e-5,np.ones(3)*5e-5,1)
cells : 64 x 64 x 64 cells : 64 x 64 x 64
size : 0.0001 x 0.0001 x 0.0001 / m³ size : 0.0001 x 0.0001 x 0.0001 m³
origin: 0.0 0.0 0.0 / m origin: 0.0 0.0 0.0 m
# materials: 2 # materials: 2
Add a cube at the origin. Add a cube at the origin.
@ -746,8 +746,8 @@ class Grid:
>>> g = damask.Grid(np.zeros([64]*3,int), np.ones(3)*1e-4) >>> g = damask.Grid(np.zeros([64]*3,int), np.ones(3)*1e-4)
>>> g.add_primitive(np.ones(3,int)*32,np.zeros(3),np.inf) >>> g.add_primitive(np.ones(3,int)*32,np.zeros(3),np.inf)
cells : 64 x 64 x 64 cells : 64 x 64 x 64
size : 0.0001 x 0.0001 x 0.0001 / m³ size : 0.0001 x 0.0001 x 0.0001 m³
origin: 0.0 0.0 0.0 / m origin: 0.0 0.0 0.0 m
# materials: 2 # materials: 2
""" """
@ -805,8 +805,8 @@ class Grid:
>>> g = damask.Grid(np.zeros([32]*3,int), np.ones(3)*1e-4) >>> g = damask.Grid(np.zeros([32]*3,int), np.ones(3)*1e-4)
>>> g.mirror('xy',True) >>> g.mirror('xy',True)
cells : 64 x 64 x 32 cells : 64 x 64 x 32
size : 0.0002 x 0.0002 x 0.0001 / m³ size : 0.0002 x 0.0002 x 0.0001 m³
origin: 0.0 0.0 0.0 / m origin: 0.0 0.0 0.0 m
# materials: 1 # materials: 1
""" """
@ -886,8 +886,8 @@ class Grid:
>>> g = damask.Grid(np.zeros([32]*3,int),np.ones(3)*1e-4) >>> g = damask.Grid(np.zeros([32]*3,int),np.ones(3)*1e-4)
>>> g.scale(g.cells*2) >>> g.scale(g.cells*2)
cells : 64 x 64 x 64 cells : 64 x 64 x 64
size : 0.0001 x 0.0001 x 0.0001 / m³ size : 0.0001 x 0.0001 x 0.0001 m³
origin: 0.0 0.0 0.0 / m origin: 0.0 0.0 0.0 m
# materials: 1 # materials: 1
""" """
@ -1036,8 +1036,8 @@ class Grid:
>>> g = damask.Grid(np.zeros([32]*3,int),np.ones(3)*1e-4) >>> g = damask.Grid(np.zeros([32]*3,int),np.ones(3)*1e-4)
>>> g.canvas([32,32,16]) >>> g.canvas([32,32,16])
cells : 33 x 32 x 16 cells : 33 x 32 x 16
size : 0.0001 x 0.0001 x 5e-05 / m³ size : 0.0001 x 0.0001 x 5e-05 m³
origin: 0.0 0.0 0.0 / m origin: 0.0 0.0 0.0 m
# materials: 1 # materials: 1
""" """

238
python/damask/_orientation.py
View File

@ -1,8 +1,10 @@
import inspect import inspect
import copy import copy
from typing import Union, Callable, List, Dict, Any, Tuple, TypeVar
import numpy as np import numpy as np
from ._typehints import FloatSequence, IntSequence, CrystalFamily, CrystalLattice
from . import Rotation from . import Rotation
from . import Crystal from . import Crystal
from . import util from . import util
@ -33,6 +35,7 @@ _parameter_doc = \
""" """
MyType = TypeVar('MyType', bound='Orientation')
class Orientation(Rotation,Crystal): class Orientation(Rotation,Crystal):
""" """
@ -93,12 +96,13 @@ class Orientation(Rotation,Crystal):
@util.extend_docstring(_parameter_doc) @util.extend_docstring(_parameter_doc)
def __init__(self, def __init__(self,
rotation = np.array([1.0,0.0,0.0,0.0]), *, rotation: Union[FloatSequence, Rotation] = np.array([1.,0.,0.,0.]),
family = None, *,
lattice = None, family: CrystalFamily = None,
a = None,b = None,c = None, lattice: CrystalLattice = None,
alpha = None,beta = None,gamma = None, a: float = None, b: float = None, c: float = None,
degrees = False): alpha: float = None, beta: float = None, gamma: float = None,
degrees: bool = False):
""" """
New orientation. New orientation.
@ -115,13 +119,13 @@ class Orientation(Rotation,Crystal):
a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma, degrees=degrees) a=a,b=b,c=c, alpha=alpha,beta=beta,gamma=gamma, degrees=degrees)
def __repr__(self): def __repr__(self) -> str:
"""Represent.""" """Represent."""
return '\n'.join([Crystal.__repr__(self), return '\n'.join([Crystal.__repr__(self),
Rotation.__repr__(self)]) Rotation.__repr__(self)])
def __copy__(self: MyType,
def __copy__(self,rotation=None): rotation: Union[FloatSequence, Rotation] = None) -> MyType:
"""Create deep copy.""" """Create deep copy."""
dup = copy.deepcopy(self) dup = copy.deepcopy(self)
if rotation is not None: if rotation is not None:
@ -131,7 +135,9 @@ class Orientation(Rotation,Crystal):
copy = __copy__ copy = __copy__
def __eq__(self,other):
def __eq__(self,
other: object) -> bool:
""" """
Equal to other. Equal to other.
@ -141,12 +147,15 @@ class Orientation(Rotation,Crystal):
Orientation to check for equality. Orientation to check for equality.
""" """
if not isinstance(other, Orientation):
return NotImplemented
matching_type = self.family == other.family and \ matching_type = self.family == other.family and \
self.lattice == other.lattice and \ self.lattice == other.lattice and \
self.parameters == other.parameters 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: object) -> bool:
""" """
Not equal to other. Not equal to other.
@ -156,10 +165,14 @@ class Orientation(Rotation,Crystal):
Orientation to check for equality. Orientation to check for equality.
""" """
return np.logical_not(self==other) return np.logical_not(self==other) if isinstance(other, Orientation) else NotImplemented
def isclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): def isclose(self: MyType,
other: MyType,
rtol: float = 1e-5,
atol: float = 1e-8,
equal_nan: bool = True) -> bool:
""" """
Report where values are approximately equal to corresponding ones of other Orientation. Report where values are approximately equal to corresponding ones of other Orientation.
@ -176,7 +189,7 @@ class Orientation(Rotation,Crystal):
Returns Returns
------- -------
mask : numpy.ndarray bool mask : numpy.ndarray of bool, shape (self.shape)
Mask indicating where corresponding orientations are close. Mask indicating where corresponding orientations are close.
""" """
@ -187,7 +200,11 @@ class Orientation(Rotation,Crystal):
def allclose(self,other,rtol=1e-5,atol=1e-8,equal_nan=True): def allclose(self: MyType,
other: MyType,
rtol: float = 1e-5,
atol: float = 1e-8,
equal_nan: bool = True) -> bool:
""" """
Test whether all values are approximately equal to corresponding ones of other Orientation. Test whether all values are approximately equal to corresponding ones of other Orientation.
@ -208,10 +225,11 @@ class Orientation(Rotation,Crystal):
Whether all values are close between both orientations. Whether all values are close between both orientations.
""" """
return np.all(self.isclose(other,rtol,atol,equal_nan)) return bool(np.all(self.isclose(other,rtol,atol,equal_nan)))
def __mul__(self,other): def __mul__(self: MyType,
other: Union[Rotation, 'Orientation']) -> MyType:
""" """
Compose this orientation with other. Compose this orientation with other.
@ -226,14 +244,15 @@ class Orientation(Rotation,Crystal):
Compound rotation self*other, i.e. first other then self rotation. Compound rotation self*other, i.e. first other then self rotation.
""" """
if isinstance(other,Orientation) or isinstance(other,Rotation): if isinstance(other, (Orientation,Rotation)):
return self.copy(rotation=Rotation.__mul__(self,Rotation(other.quaternion))) return self.copy(Rotation(self.quaternion)*Rotation(other.quaternion))
else: else:
raise TypeError('use "O@b", i.e. matmul, to apply Orientation "O" to object "b"') raise TypeError('use "O@b", i.e. matmul, to apply Orientation "O" to object "b"')
@staticmethod @staticmethod
def _split_kwargs(kwargs,target): def _split_kwargs(kwargs: Dict[str, Any],
target: Callable) -> Tuple[Dict[str, Any], ...]:
""" """
Separate keyword arguments in 'kwargs' targeted at 'target' from general keyword arguments of Orientation objects. Separate keyword arguments in 'kwargs' targeted at 'target' from general keyword arguments of Orientation objects.
@ -252,7 +271,7 @@ class Orientation(Rotation,Crystal):
Valid keyword arguments of Orientation object. Valid keyword arguments of Orientation object.
""" """
kws = () kws: Tuple[Dict[str, Any], ...] = ()
for t in (target,Orientation.__init__): for t in (target,Orientation.__init__):
kws += ({key: kwargs[key] for key in set(inspect.signature(t).parameters) & set(kwargs)},) kws += ({key: kwargs[key] for key in set(inspect.signature(t).parameters) & set(kwargs)},)
@ -264,105 +283,108 @@ class Orientation(Rotation,Crystal):
@classmethod @classmethod
@util.extended_docstring(Rotation.from_random,_parameter_doc) @util.extended_docstring(Rotation.from_random, _parameter_doc)
def from_random(cls,**kwargs): def from_random(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random)
return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_quaternion,_parameter_doc) @util.extended_docstring(Rotation.from_quaternion,_parameter_doc)
def from_quaternion(cls,**kwargs): def from_quaternion(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion)
return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc) @util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc)
def from_Euler_angles(cls,**kwargs): def from_Euler_angles(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles)
return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_axis_angle,_parameter_doc) @util.extended_docstring(Rotation.from_axis_angle,_parameter_doc)
def from_axis_angle(cls,**kwargs): def from_axis_angle(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle)
return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_basis,_parameter_doc) @util.extended_docstring(Rotation.from_basis,_parameter_doc)
def from_basis(cls,**kwargs): def from_basis(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis)
return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_matrix,_parameter_doc) @util.extended_docstring(Rotation.from_matrix,_parameter_doc)
def from_matrix(cls,**kwargs): def from_matrix(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix)
return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc) @util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc)
def from_Rodrigues_vector(cls,**kwargs): def from_Rodrigues_vector(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector)
return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_homochoric,_parameter_doc) @util.extended_docstring(Rotation.from_homochoric,_parameter_doc)
def from_homochoric(cls,**kwargs): def from_homochoric(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric)
return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_cubochoric,_parameter_doc) @util.extended_docstring(Rotation.from_cubochoric,_parameter_doc)
def from_cubochoric(cls,**kwargs): def from_cubochoric(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric)
return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_spherical_component,_parameter_doc) @util.extended_docstring(Rotation.from_spherical_component,_parameter_doc)
def from_spherical_component(cls,**kwargs): def from_spherical_component(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component)
return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extended_docstring(Rotation.from_fiber_component,_parameter_doc) @util.extended_docstring(Rotation.from_fiber_component,_parameter_doc)
def from_fiber_component(cls,**kwargs): def from_fiber_component(cls, **kwargs) -> 'Orientation':
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component) kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component)
return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori) return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori)
@classmethod @classmethod
@util.extend_docstring(_parameter_doc) @util.extend_docstring(_parameter_doc)
def from_directions(cls,uvw,hkl,**kwargs): def from_directions(cls,
uvw: FloatSequence,
hkl: FloatSequence,
**kwargs) -> 'Orientation':
""" """
Initialize orientation object from two crystallographic directions. Initialize orientation object from two crystallographic directions.
Parameters Parameters
---------- ----------
uvw : list, numpy.ndarray of shape (...,3) uvw : numpy.ndarray, shape (...,3)
lattice direction aligned with lab frame x-direction. Lattice direction aligned with lab frame x-direction.
hkl : list, numpy.ndarray of shape (...,3) hkl : numpy.ndarray, shape (...,3)
lattice plane normal aligned with lab frame z-direction. Lattice plane normal aligned with lab frame z-direction.
""" """
o = cls(**kwargs) o = cls(**kwargs)
x = o.to_frame(uvw=uvw) x = o.to_frame(uvw=uvw)
z = o.to_frame(hkl=hkl) z = o.to_frame(hkl=hkl)
om = np.stack([x,np.cross(z,x),z],axis=-2) om = np.stack([x,np.cross(z,x),z],axis=-2)
return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True)))) return o.copy(Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True))))
@property @property
def equivalent(self): def equivalent(self: MyType) -> MyType:
""" """
Orientations that are symmetrically equivalent. Orientations that are symmetrically equivalent.
@ -372,11 +394,11 @@ class Orientation(Rotation,Crystal):
""" """
sym_ops = self.symmetry_operations sym_ops = self.symmetry_operations
o = sym_ops.broadcast_to(sym_ops.shape+self.shape,mode='right') 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(o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
@property @property
def reduced(self): def reduced(self: MyType) -> MyType:
"""Select symmetrically equivalent orientation that falls into fundamental zone according to symmetry.""" """Select symmetrically equivalent orientation that falls into fundamental zone according to symmetry."""
eq = self.equivalent eq = self.equivalent
ok = eq.in_FZ ok = eq.in_FZ
@ -387,13 +409,13 @@ class Orientation(Rotation,Crystal):
@property @property
def in_FZ(self): def in_FZ(self) -> Union[np.bool_, np.ndarray]:
""" """
Check whether orientation falls into fundamental zone of own symmetry. Check whether orientation falls into fundamental zone of own symmetry.
Returns Returns
------- -------
in : numpy.ndarray of bool, quaternion.shape in : numpy.ndarray of bool, shape (self.shape)
Whether Rodrigues-Frank vector falls into fundamental zone. Whether Rodrigues-Frank vector falls into fundamental zone.
Notes Notes
@ -431,13 +453,13 @@ class Orientation(Rotation,Crystal):
@property @property
def in_disorientation_FZ(self): def in_disorientation_FZ(self) -> np.ndarray:
""" """
Check whether orientation falls into fundamental zone of disorientations. Check whether orientation falls into fundamental zone of disorientations.
Returns Returns
------- -------
in : numpy.ndarray of bool, quaternion.shape in : numpy.ndarray of bool, shape (self.shape)
Whether Rodrigues-Frank vector falls into disorientation FZ. Whether Rodrigues-Frank vector falls into disorientation FZ.
References References
@ -471,8 +493,9 @@ class Orientation(Rotation,Crystal):
else: else:
return np.ones_like(rho[...,0],dtype=bool) return np.ones_like(rho[...,0],dtype=bool)
def disorientation(self,
def disorientation(self,other,return_operators=False): other: 'Orientation',
return_operators: bool = False) -> object:
""" """
Calculate disorientation between myself and given other orientation. Calculate disorientation between myself and given other orientation.
@ -490,7 +513,7 @@ class Orientation(Rotation,Crystal):
------- -------
disorientation : Orientation disorientation : Orientation
Disorientation between self and other. Disorientation between self and other.
operators : numpy.ndarray int of shape (...,2), conditional operators : numpy.ndarray of int, shape (...,2), conditional
Index of symmetrically equivalent orientation that rotated vector to the SST. Index of symmetrically equivalent orientation that rotated vector to the SST.
Notes Notes
@ -530,7 +553,7 @@ class Orientation(Rotation,Crystal):
raise NotImplementedError('disorientation between different crystal families') raise NotImplementedError('disorientation between different crystal families')
blend = util.shapeblender(self.shape,other.shape) blend = util.shapeblender(self.shape,other.shape)
s = self.equivalent s = self.equivalent
o = other.equivalent o = other.equivalent
s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right') s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right')
@ -557,13 +580,15 @@ class Orientation(Rotation,Crystal):
) )
def average(self,weights=None,return_cloud=False): def average(self,
weights: FloatSequence = None,
return_cloud: bool = False):
""" """
Return orientation average over last dimension. Return orientation average over last dimension.
Parameters Parameters
---------- ----------
weights : numpy.ndarray, optional weights : numpy.ndarray, shape (self.shape), optional
Relative weights of orientations. Relative weights of orientations.
return_cloud : bool, optional return_cloud : bool, optional
Return the set of symmetrically equivalent orientations that was used in averaging. Return the set of symmetrically equivalent orientations that was used in averaging.
@ -583,31 +608,30 @@ class Orientation(Rotation,Crystal):
""" """
eq = self.equivalent eq = self.equivalent
m = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,)) m = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,)) # type: ignore
.broadcast_to(eq.shape))\ .broadcast_to(eq.shape)).as_axis_angle()[...,3] # type: ignore
.as_axis_angle()[...,3]
r = Rotation(np.squeeze(np.take_along_axis(eq.quaternion, r = Rotation(np.squeeze(np.take_along_axis(eq.quaternion,
np.argmin(m,axis=0)[np.newaxis,...,np.newaxis], np.argmin(m,axis=0)[np.newaxis,...,np.newaxis],
axis=0), axis=0),
axis=0)) axis=0))
return ( return ((self.copy(Rotation(r).average(weights)),self.copy(Rotation(r))) if return_cloud else
(self.copy(rotation=Rotation(r).average(weights)), self.copy(Rotation(r).average(weights))
self.copy(rotation=Rotation(r)))
if return_cloud else
self.copy(rotation=Rotation(r).average(weights))
) )
def to_SST(self,vector,proper=False,return_operators=False): def to_SST(self,
vector: FloatSequence,
proper: bool = False,
return_operators: bool = False) -> np.ndarray:
""" """
Rotate vector to ensure it falls into (improper or proper) standard stereographic triangle of crystal symmetry. Rotate vector to ensure it falls into (improper or proper) standard stereographic triangle of crystal symmetry.
Parameters Parameters
---------- ----------
vector : numpy.ndarray of shape (...,3) vector : numpy.ndarray, shape (...,3)
Lab frame vector to align with crystal frame direction. Lab frame vector to align with crystal frame direction.
Shape of vector blends with shape of own rotation array. Shape of vector blends with shape of own rotation array.
For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs. For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
proper : bool, optional proper : bool, optional
Consider only vectors with z >= 0, hence combine two neighboring SSTs. Consider only vectors with z >= 0, hence combine two neighboring SSTs.
Defaults to False. Defaults to False.
@ -617,15 +641,18 @@ class Orientation(Rotation,Crystal):
Returns Returns
------- -------
vector_SST : numpy.ndarray of shape (...,3) vector_SST : numpy.ndarray, shape (...,3)
Rotated vector falling into SST. Rotated vector falling into SST.
operators : numpy.ndarray int of shape (...), conditional operators : numpy.ndarray of int, shape (...), conditional
Index of symmetrically equivalent orientation that rotated vector to SST. Index of symmetrically equivalent orientation that rotated vector to SST.
""" """
vector_ = np.array(vector,float)
if vector_.shape[-1] != 3:
raise ValueError('input is not a field of three-dimensional vectors')
eq = self.equivalent eq = self.equivalent
blend = util.shapeblender(eq.shape,np.array(vector).shape[:-1]) blend = util.shapeblender(eq.shape,vector_.shape[:-1])
poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(np.array(vector),blend+(3,)) poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(vector_,blend+(3,))
ok = self.in_SST(poles,proper=proper) ok = self.in_SST(poles,proper=proper)
ok &= np.cumsum(ok,axis=0) == 1 ok &= np.cumsum(ok,axis=0) == 1
loc = np.where(ok) loc = np.where(ok)
@ -637,13 +664,15 @@ class Orientation(Rotation,Crystal):
) )
def in_SST(self,vector,proper=False): def in_SST(self,
vector: FloatSequence,
proper: bool = False) -> Union[np.bool_, np.ndarray]:
""" """
Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry. Check whether given crystal frame vector falls into standard stereographic triangle of own symmetry.
Parameters Parameters
---------- ----------
vector : numpy.ndarray of shape (...,3) vector : numpy.ndarray, shape (...,3)
Vector to check. Vector to check.
proper : bool, optional proper : bool, optional
Consider only vectors with z >= 0, hence combine two neighboring SSTs. Consider only vectors with z >= 0, hence combine two neighboring SSTs.
@ -655,39 +684,43 @@ class Orientation(Rotation,Crystal):
Whether vector falls into SST. Whether vector falls into SST.
""" """
if not isinstance(vector,np.ndarray) or vector.shape[-1] != 3: vector_ = np.array(vector,float)
if vector_.shape[-1] != 3:
raise ValueError('input is not a field of three-dimensional vectors') raise ValueError('input is not a field of three-dimensional vectors')
if self.standard_triangle is None: # direct exit for no symmetry if self.standard_triangle is None: # direct exit for no symmetry
return np.ones_like(vector[...,0],bool) return np.ones_like(vector_[...,0],bool)
if proper: if proper:
components_proper = np.around(np.einsum('...ji,...i', components_proper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['proper'], vector.shape+(3,)), np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)),
vector), 12) vector_), 12)
components_improper = np.around(np.einsum('...ji,...i', components_improper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)), np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
vector), 12) vector_), 12)
return np.all(components_proper >= 0.0,axis=-1) \ return np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 0.0,axis=-1) | np.all(components_improper >= 0.0,axis=-1)
else: else:
components = np.around(np.einsum('...ji,...i', components = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector.shape+(3,)), np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
np.block([vector[...,:2],np.abs(vector[...,2:3])])), 12) np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
return np.all(components >= 0.0,axis=-1) return np.all(components >= 0.0,axis=-1)
def IPF_color(self,vector,in_SST=True,proper=False): def IPF_color(self,
vector: FloatSequence,
in_SST: bool = True,
proper: bool = False) -> np.ndarray:
""" """
Map vector to RGB color within standard stereographic triangle of own symmetry. Map vector to RGB color within standard stereographic triangle of own symmetry.
Parameters Parameters
---------- ----------
vector : numpy.ndarray of shape (...,3) vector : numpy.ndarray, shape (...,3)
Vector to colorize. Vector to colorize.
Shape of vector blends with shape of own rotation array. Shape of vector blends with shape of own rotation array.
For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs. For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
in_SST : bool, optional in_SST : bool, optional
Consider symmetrically equivalent orientations such that poles are located in SST. Consider symmetrically equivalent orientations such that poles are located in SST.
Defaults to True. Defaults to True.
@ -697,7 +730,7 @@ class Orientation(Rotation,Crystal):
Returns Returns
------- -------
rgb : numpy.ndarray of shape (...,3) rgb : numpy.ndarray, shape (...,3)
RGB array of IPF colors. RGB array of IPF colors.
Examples Examples
@ -721,35 +754,35 @@ class Orientation(Rotation,Crystal):
if proper: if proper:
components_proper = np.around(np.einsum('...ji,...i', components_proper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)), np.broadcast_to(self.standard_triangle['proper'], vector_.shape+(3,)),
vector_), 12) vector_), 12)
components_improper = np.around(np.einsum('...ji,...i', components_improper = np.around(np.einsum('...ji,...i',
np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)), np.broadcast_to(self.standard_triangle['improper'], vector_.shape+(3,)),
vector_), 12) vector_), 12)
in_SST = np.all(components_proper >= 0.0,axis=-1) \ in_SST_ = np.all(components_proper >= 0.0,axis=-1) \
| np.all(components_improper >= 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 = np.where((in_SST_ & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis],
components_proper,components_improper) components_proper,components_improper)
else: else:
components = np.around(np.einsum('...ji,...i', components = np.around(np.einsum('...ji,...i',
np.broadcast_to(self .standard_triangle['improper'], vector_.shape+(3,)), np.broadcast_to(self .standard_triangle['improper'], vector_.shape+(3,)),
np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12) np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12)
in_SST = np.all(components >= 0.0,axis=-1) in_SST_ = np.all(components >= 0.0,axis=-1)
with np.errstate(invalid='ignore',divide='ignore'): with np.errstate(invalid='ignore',divide='ignore'):
rgb = (components/np.linalg.norm(components,axis=-1,keepdims=True))**0.5 # smoothen color ramps 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.clip(rgb,0.,1.) # clip intensity
rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1 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 rgb[np.broadcast_to(~in_SST_[...,np.newaxis],rgb.shape)] = 0.0
return rgb return rgb
@property @property
def symmetry_operations(self): def symmetry_operations(self) -> Rotation:
"""Symmetry operations as Rotations.""" """Symmetry operations as Rotations."""
_symmetry_operations = { _symmetry_operations: Dict[CrystalFamily, List] = {
'cubic': [ 'cubic': [
[ 1.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, 1.0, 0.0, 0.0 ],
@ -819,22 +852,25 @@ class Orientation(Rotation,Crystal):
#################################################################################################### ####################################################################################################
# functions that require lattice, not just family # functions that require lattice, not just family
def to_pole(self,*,uvw=None,hkl=None,with_symmetry=False): def to_pole(self, *,
uvw: FloatSequence = None,
hkl: FloatSequence = None,
with_symmetry: bool = False) -> np.ndarray:
""" """
Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl). Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl).
Parameters Parameters
---------- ----------
uvw|hkl : numpy.ndarray of shape (...,3) uvw|hkl : numpy.ndarray, shape (...,3)
Miller indices of crystallographic direction or plane normal. Miller indices of crystallographic direction or plane normal.
Shape of vector blends with shape of own rotation array. Shape of vector blends with shape of own rotation array.
For example, a rotation array of shape (3,2) and a (2,4) vector array result in (3,2,4) outputs. For example, a rotation array, shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
with_symmetry : bool, optional with_symmetry : bool, optional
Calculate all N symmetrically equivalent vectors. Calculate all N symmetrically equivalent vectors.
Returns Returns
------- -------
vector : numpy.ndarray of shape (...,3) or (...,N,3) vector : numpy.ndarray, shape (...,3) or (...,N,3)
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.
""" """
@ -846,23 +882,24 @@ class Orientation(Rotation,Crystal):
blend += sym_ops.shape blend += sym_ops.shape
v = sym_ops.broadcast_to(shape) \ v = sym_ops.broadcast_to(shape) \
@ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,)) @ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,))
return ~(self.broadcast_to(blend)) \ return ~(self.broadcast_to(blend))@ np.broadcast_to(v,blend+(3,))
@ np.broadcast_to(v,blend+(3,))
def Schmid(self,*,N_slip=None,N_twin=None): def Schmid(self, *,
N_slip: IntSequence = None,
N_twin: IntSequence = None) -> np.ndarray:
u""" u"""
Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems. Calculate Schmid matrix P = d ⨂ n in the lab frame for selected deformation systems.
Parameters Parameters
---------- ----------
N_slip|N_twin : iterable of int N_slip|N_twin : '*' or iterable of int
Number of deformation systems per family of the deformation system. Number of deformation systems per family of the deformation system.
Use '*' to select all. Use '*' to select all.
Returns Returns
------- -------
P : numpy.ndarray of shape (N,...,3,3) P : numpy.ndarray, shape (N,...,3,3)
Schmid matrix for each of the N deformation systems. Schmid matrix for each of the N deformation systems.
Examples Examples
@ -887,6 +924,8 @@ class Orientation(Rotation,Crystal):
(self.kinematics('twin'),N_twin) (self.kinematics('twin'),N_twin)
if active == '*': active = [len(a) for a in kinematics['direction']] if active == '*': active = [len(a) for a in kinematics['direction']]
if not active:
raise ValueError('Schmid matrix not defined')
d = self.to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)])) 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 = 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.einsum('...i,...j',d/np.linalg.norm(d,axis=1,keepdims=True),
@ -897,7 +936,8 @@ class Orientation(Rotation,Crystal):
@ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape) @ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape)
def related(self,model): def related(self: MyType,
model: str) -> MyType:
""" """
Orientations derived from the given relationship. Orientations derived from the given relationship.

38
python/damask/_result.py
View File

@ -168,17 +168,21 @@ class Result:
def __repr__(self): def __repr__(self):
"""Show summary of file content.""" """Show summary of file content."""
with h5py.File(self.fname,'r') as f:
header = [f'Created by {f.attrs["creator"]}',
f' on {f.attrs["created"]}',
f' executing "{f.attrs["call"]}"']
visible_increments = self.visible['increments'] visible_increments = self.visible['increments']
first = self.view(increments=visible_increments[0:1]).list_data() first = self.view(increments=visible_increments[0:1]).list_data()
last = '' if len(visible_increments) < 2 else \ last = [] if len(visible_increments) < 2 else \
self.view(increments=visible_increments[-1:]).list_data() self.view(increments=visible_increments[-1:]).list_data()
in_between = '' if len(visible_increments) < 3 else \ in_between = [] if len(visible_increments) < 3 else \
''.join([f'\n{inc}\n ...\n' for inc in visible_increments[1:-1]]) [f'\n{inc}\n ...' for inc in visible_increments[1:-1]]
return util.srepr(first + in_between + last) return util.srepr([util.deemph(header)] + first + in_between + last)
def _manage_view(self,action,what,datasets): def _manage_view(self,action,what,datasets):
@ -537,24 +541,32 @@ class Result:
def list_data(self): def list_data(self):
"""Return information on all active datasets in the file.""" """
msg = '' Collect information on all active datasets in the file.
Returns
-------
data : list of str
Line-formatted information about active datasets.
"""
msg = []
with h5py.File(self.fname,'r') as f: with h5py.File(self.fname,'r') as f:
for inc in self.visible['increments']: for inc in self.visible['increments']:
msg = ''.join([msg,f'\n{inc} ({self.times[self.increments.index(inc)]}s)\n']) msg += [f'\n{inc} ({self.times[self.increments.index(inc)]} s)']
for ty in ['phase','homogenization']: for ty in ['phase','homogenization']:
msg = ' '.join([msg,f'{ty}\n']) msg += [f' {ty}']
for label in self.visible[ty+'s']: for label in self.visible[ty+'s']:
msg = ' '.join([msg,f'{label}\n']) msg += [f' {label}']
for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()): for field in _match(self.visible['fields'],f['/'.join([inc,ty,label])].keys()):
msg = ' '.join([msg,f'{field}\n']) msg += [f' {field}']
for d in f['/'.join([inc,ty,label,field])].keys(): for d in f['/'.join([inc,ty,label,field])].keys():
dataset = f['/'.join([inc,ty,label,field,d])] dataset = f['/'.join([inc,ty,label,field,d])]
unit = f' / {dataset.attrs["unit"]}' if h5py3 else \ unit = dataset.attrs["unit"] if h5py3 else \
f' / {dataset.attrs["unit"].decode()}' dataset.attrs["unit"].decode()
description = dataset.attrs['description'] if h5py3 else \ description = dataset.attrs['description'] if h5py3 else \
dataset.attrs['description'].decode() dataset.attrs['description'].decode()
msg = ' '.join([msg,f'{d}{unit}: {description}\n']) msg += [f' {d} / {unit}: {description}']
return msg return msg

424
python/damask/_rotation.py
View File

File diff suppressed because it is too large Load Diff

6
python/damask/_table.py
View File

@ -194,7 +194,7 @@ class Table:
Returns Returns
------- -------
mask : numpy.ndarray bool mask : numpy.ndarray of bool
Mask indicating where corresponding table values are close. Mask indicating where corresponding table values are close.
""" """
@ -263,10 +263,8 @@ class Table:
f.seek(0) f.seek(0)
comments = [] comments = []
line = f.readline().strip() while (line := f.readline().strip()).startswith('#'):
while line.startswith('#'):
comments.append(line.lstrip('#').strip()) comments.append(line.lstrip('#').strip())
line = f.readline().strip()
labels = line.split() labels = line.split()
shapes = {} shapes = {}

5
python/damask/_typehints.py
View File

@ -1,6 +1,6 @@
"""Functionality for typehints.""" """Functionality for typehints."""
from typing import Sequence, Union, TextIO from typing import Sequence, Union, Literal, TextIO
from pathlib import Path from pathlib import Path
import numpy as np import numpy as np
@ -9,6 +9,9 @@ import numpy as np
FloatSequence = Union[np.ndarray,Sequence[float]] FloatSequence = Union[np.ndarray,Sequence[float]]
IntSequence = Union[np.ndarray,Sequence[int]] IntSequence = Union[np.ndarray,Sequence[int]]
FileHandle = Union[TextIO, str, Path] FileHandle = Union[TextIO, str, Path]
CrystalFamily = Union[None,Literal['triclinic', 'monoclinic', 'orthorhombic', 'tetragonal', 'hexagonal', 'cubic']]
CrystalLattice = Union[None,Literal['aP', 'mP', 'mS', 'oP', 'oS', 'oI', 'oF', 'tP', 'tI', 'hP', 'cP', 'cI', 'cF']]
CrystalKinematics = Literal['slip', 'twin']
NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator] NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.Generator]
# BitGenerator does not exists in older numpy versions # BitGenerator does not exists in older numpy versions
#NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.BitGenerator, np.random.Generator] #NumpyRngSeed = Union[int, IntSequence, np.random.SeedSequence, np.random.BitGenerator, np.random.Generator]

8
python/damask/util.py
View File

@ -9,7 +9,7 @@ import re
import fractions import fractions
from collections import abc from collections import abc
from functools import reduce from functools import reduce
from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal, SupportsIndex, Sequence from typing import Union, Tuple, Iterable, Callable, Dict, List, Any, Literal
from pathlib import Path from pathlib import Path
import numpy as np import numpy as np
@ -428,7 +428,7 @@ def hybrid_IA(dist: np.ndarray,
def shapeshifter(fro: Tuple[int, ...], def shapeshifter(fro: Tuple[int, ...],
to: Tuple[int, ...], to: Tuple[int, ...],
mode: Literal['left','right'] = 'left', mode: Literal['left','right'] = 'left',
keep_ones: bool = False) -> Sequence[SupportsIndex]: keep_ones: bool = False) -> Tuple[int, ...]:
""" """
Return dimensions that reshape 'fro' to become broadcastable to 'to'. Return dimensions that reshape 'fro' to become broadcastable to 'to'.
@ -491,7 +491,7 @@ def shapeshifter(fro: Tuple[int, ...],
def shapeblender(a: Tuple[int, ...], def shapeblender(a: Tuple[int, ...],
b: Tuple[int, ...]) -> Sequence[SupportsIndex]: b: Tuple[int, ...]) -> Tuple[int, ...]:
""" """
Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'. Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
@ -794,7 +794,7 @@ class ProgressBar:
fraction = (iteration+1) / self.total fraction = (iteration+1) / self.total
if filled_length := int(self.bar_length * fraction) > int(self.bar_length * self.fraction_last) or \ if (filled_length := int(self.bar_length * fraction)) > int(self.bar_length * self.fraction_last) or \
datetime.datetime.now() - self.time_last_update > datetime.timedelta(seconds=10): datetime.datetime.now() - self.time_last_update > datetime.timedelta(seconds=10):
self.time_last_update = datetime.datetime.now() self.time_last_update = datetime.datetime.now()
bar = '█' * filled_length + '░' * (self.bar_length - filled_length) bar = '█' * filled_length + '░' * (self.bar_length - filled_length)

8
python/tests/test_Rotation.py
View File

@ -313,13 +313,13 @@ def ho2ax(ho):
if iszero(hmag_squared): if iszero(hmag_squared):
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
else: else:
hm = hmag_squared hm = np.ones_like(hmag_squared)
# convert the magnitude to the rotation angle # convert the magnitude to the rotation angle
s = tfit[0] + tfit[1] * hmag_squared s = 0.0
for i in range(2,16): for t in tfit:
s += t * hm
hm *= hmag_squared hm *= hmag_squared
s += tfit[i] * hm
ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0)))
return ax return ax

4
src/discretization.f90
View File

@ -12,7 +12,8 @@ module discretization
integer, public, protected :: & integer, public, protected :: &
discretization_nIPs, & discretization_nIPs, &
discretization_Nelems discretization_Nelems, &
discretization_Ncells
integer, public, protected, dimension(:), allocatable :: & integer, public, protected, dimension(:), allocatable :: &
discretization_materialAt !ToDo: discretization_ID_material discretization_materialAt !ToDo: discretization_ID_material
@ -53,6 +54,7 @@ subroutine discretization_init(materialAt,&
discretization_Nelems = size(materialAt,1) discretization_Nelems = size(materialAt,1)
discretization_nIPs = size(IPcoords0,2)/discretization_Nelems discretization_nIPs = size(IPcoords0,2)/discretization_Nelems
discretization_Ncells = discretization_Nelems*discretization_nIPs
discretization_materialAt = materialAt discretization_materialAt = materialAt

4
src/grid/discretization_grid.f90
View File

@ -97,8 +97,8 @@ subroutine discretization_grid_init(restart)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
print'(/,1x,a,i0,a,i0,a,i0)', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3) print'(/,1x,a,i0,a,i0,a,i0)', 'cells: ', cells(1), ' × ', cells(2), ' × ', cells(3)
print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' / m³' print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'size: ', geomSize(1), ' × ', geomSize(2), ' × ', geomSize(3), ' m³'
print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'origin: ', origin(1), ' ', origin(2), ' ', origin(3), ' / m' print '(1x,a,es8.2,a,es8.2,a,es8.2,a)', 'origin: ', origin(1), ' ', origin(2), ' ', origin(3), ' m'
if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)') if (worldsize>cells(3)) call IO_error(894, ext_msg='number of processes exceeds cells(3)')

13
src/homogenization_mechanical.f90
View File

@ -131,20 +131,15 @@ module subroutine mechanical_homogenize(Delta_t,ce)
integer :: co integer :: co
homogenization_P(1:3,1:3,ce) = phase_P(1,ce) homogenization_P(1:3,1:3,ce) = phase_P(1,ce)*material_v(1,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce)*material_v(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) & homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) &
+ phase_P(co,ce) + phase_P(co,ce)*material_v(co,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) & homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
+ phase_mechanical_dPdF(Delta_t,co,ce) + phase_mechanical_dPdF(Delta_t,co,ce)*material_v(co,ce)
end do end do
homogenization_P(1:3,1:3,ce) = homogenization_P(1:3,1:3,ce) &
/ real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
/ real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end subroutine mechanical_homogenize end subroutine mechanical_homogenize

18
src/homogenization_thermal.f90
View File

@ -105,13 +105,11 @@ module function homogenization_mu_T(ce) result(mu)
integer :: co integer :: co
mu = phase_mu_T(1,ce) mu = phase_mu_T(1,ce)*material_v(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
mu = mu + phase_mu_T(co,ce) mu = mu + phase_mu_T(co,ce)*material_v(co,ce)
end do end do
mu = mu / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function homogenization_mu_T end function homogenization_mu_T
@ -126,13 +124,11 @@ module function homogenization_K_T(ce) result(K)
integer :: co integer :: co
K = phase_K_T(1,ce) K = phase_K_T(1,ce)*material_v(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
K = K + phase_K_T(co,ce) K = K + phase_K_T(co,ce)*material_v(co,ce)
end do end do
K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function homogenization_K_T end function homogenization_K_T
@ -147,13 +143,11 @@ module function homogenization_f_T(ce) result(f)
integer :: co integer :: co
f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce)) f = phase_f_T(material_phaseID(1,ce),material_phaseEntry(1,ce))*material_v(1,ce)
do co = 2, homogenization_Nconstituents(material_homogenizationID(ce)) do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce)) f = f + phase_f_T(material_phaseID(co,ce),material_phaseEntry(co,ce))*material_v(co,ce)
end do end do
f = f/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function homogenization_f_T end function homogenization_f_T

39
src/lattice.f90
View File

@ -528,22 +528,22 @@ end function lattice_C66_twin
!> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation !> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66,lattice_target, & function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
cOverA_trans,a_bcc,a_fcc) cOverA_trans,a_fcc,a_bcc)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol)
real(pReal), dimension(6,6), intent(in) :: C_parent66 real(pReal), dimension(6,6), intent(in) :: C_parent66
real(pReal), optional, intent(in) :: cOverA_trans, a_fcc, a_bcc
real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans
real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66 real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66
real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S
type(tRotation) :: R type(tRotation) :: R
real(pReal) :: a_bcc, a_fcc, cOverA_trans
integer :: i integer :: i
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! elasticity matrix of the target phase in cube orientation ! elasticity matrix of the target phase in cube orientation
if (lattice_target == 'hP') then if (lattice_target == 'hP' .and. present(cOverA_trans)) then
! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19) ! https://doi.org/10.1063/1.1663858 eq. (16), eq. (18), eq. (19)
! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48) ! https://doi.org/10.1016/j.actamat.2016.07.032 eq. (47), eq. (48)
if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) & if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) &
@ -562,7 +562,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(3,3) = C_bar66(3,3)
C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2)))
C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP') C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP')
elseif (lattice_target == 'cI') then elseif (lattice_target == 'cI' .and. present(a_fcc) .and. present(a_bcc)) then
if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) & if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) &
call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target)) call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target))
C_target_unrotated66 = C_parent66 C_target_unrotated66 = C_parent66
@ -1484,26 +1484,27 @@ end function lattice_SchmidMatrix_twin
!> @brief Schmid matrix for transformation !> @brief Schmid matrix for transformation
!> details only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix) function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_fcc,a_bcc) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol) character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol)
real(pReal), intent(in) :: cOverA !< c/a ratio real(pReal), optional, intent(in) :: cOverA, a_bcc, a_fcc
real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix
real(pReal), dimension(3,3,sum(Ntrans)) :: devNull real(pReal), dimension(3,3,sum(Ntrans)) :: devNull
real(pReal) :: a_bcc, a_fcc
if (lattice_target /= 'cI' .and. lattice_target /= 'hP') &
call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
if (lattice_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & if (lattice_target == 'hP' .and. present(cOverA)) then
if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) &
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA=cOverA)
if (lattice_target == 'cI' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) & else if (lattice_target == 'cI' .and. present(a_fcc) .and. present(a_bcc)) then
if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) &
call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target)) call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,a_fcc=a_fcc,a_bcc=a_bcc)
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc) else
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
end if
end function lattice_SchmidMatrix_trans end function lattice_SchmidMatrix_trans
@ -1971,7 +1972,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
real(pReal), dimension(3,3,sum(Ntrans)), intent(out) :: & real(pReal), dimension(3,3,sum(Ntrans)), intent(out) :: &
Q, & !< Total rotation: Q = R*B Q, & !< Total rotation: Q = R*B
S !< Eigendeformation tensor for phase transformation S !< Eigendeformation tensor for phase transformation
real(pReal), intent(in) :: & real(pReal), optional, intent(in) :: &
cOverA, & !< c/a for target hex lattice cOverA, & !< c/a for target hex lattice
a_bcc, & !< lattice parameter a for bcc target lattice a_bcc, & !< lattice parameter a for bcc target lattice
a_fcc !< lattice parameter a for fcc parent lattice a_fcc !< lattice parameter a for fcc parent lattice
@ -2050,7 +2051,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
0.0, 0.0, 1.0, 45.0 & 0.0, 0.0, 1.0, 45.0 &
],shape(FCCTOBCC_BAINROT)) ],shape(FCCTOBCC_BAINROT))
if (a_bcc > 0.0_pReal .and. a_fcc > 0.0_pReal .and. dEq0(cOverA)) then ! fcc -> bcc if (present(a_bcc) .and. present(a_fcc)) then
do i = 1,sum(Ntrans) do i = 1,sum(Ntrans)
call R%fromAxisAngle(FCCTOBCC_SYSTEMTRANS(:,i),degrees=.true.,P=1) call R%fromAxisAngle(FCCTOBCC_SYSTEMTRANS(:,i),degrees=.true.,P=1)
call B%fromAxisAngle(FCCTOBCC_BAINROT(:,i), degrees=.true.,P=1) call B%fromAxisAngle(FCCTOBCC_BAINROT(:,i), degrees=.true.,P=1)
@ -2064,7 +2065,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
Q(1:3,1:3,i) = matmul(R%asMatrix(),B%asMatrix()) Q(1:3,1:3,i) = matmul(R%asMatrix(),B%asMatrix())
S(1:3,1:3,i) = matmul(R%asMatrix(),U) - MATH_I3 S(1:3,1:3,i) = matmul(R%asMatrix(),U) - MATH_I3
enddo enddo
elseif (cOverA > 0.0_pReal .and. dEq0(a_bcc)) then ! fcc -> hex else if (present(cOverA)) then
ss = MATH_I3 ss = MATH_I3
sd = MATH_I3 sd = MATH_I3
ss(1,3) = sqrt(2.0_pReal)/4.0_pReal ss(1,3) = sqrt(2.0_pReal)/4.0_pReal
@ -2078,10 +2079,10 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
Q(1:3,2,i) = y Q(1:3,2,i) = y
Q(1:3,3,i) = z Q(1:3,3,i) = z
S(1:3,1:3,i) = matmul(Q(1:3,1:3,i), matmul(matmul(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only S(1:3,1:3,i) = matmul(Q(1:3,1:3,i), matmul(matmul(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only
enddo end do
else else
call IO_error(132,ext_msg='buildTransformationSystem') call IO_error(132,ext_msg='buildTransformationSystem')
endif end if
end subroutine buildTransformationSystem end subroutine buildTransformationSystem

21
src/material.f90
View File

@ -45,6 +45,9 @@ module material
material_phaseID, & ! TODO: rename to material_ID_phase material_phaseID, & ! TODO: rename to material_ID_phase
material_phaseEntry ! TODO: rename to material_entry_phase material_phaseEntry ! TODO: rename to material_entry_phase
real(pReal), dimension(:,:), allocatable, public, protected :: &
material_v ! fraction
public :: & public :: &
material_init material_init
@ -123,11 +126,13 @@ subroutine parse()
allocate(counterPhase(phases%length),source=0) allocate(counterPhase(phases%length),source=0)
allocate(counterHomogenization(homogenizations%length),source=0) allocate(counterHomogenization(homogenizations%length),source=0)
allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0) allocate(material_homogenizationID(discretization_Ncells),source=0)
allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0) allocate(material_homogenizationEntry(discretization_Ncells),source=0)
allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseID(homogenization_maxNconstituents,discretization_Ncells),source=0)
allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_Ncells),source=0)
allocate(material_v(homogenization_maxNconstituents,discretization_Ncells),source=0.0_pReal)
do el = 1, discretization_Nelems do el = 1, discretization_Nelems
material => materials%get(discretization_materialAt(el)) material => materials%get(discretization_materialAt(el))
@ -140,11 +145,11 @@ subroutine parse()
material_homogenizationEntry(ce) = counterHomogenization(ho) material_homogenizationEntry(ce) = counterHomogenization(ho)
end do end do
v = 0.0_pReal
constituents => material%get('constituents') constituents => material%get('constituents')
do co = 1, constituents%length do co = 1, constituents%length
constituent => constituents%get(co) constituent => constituents%get(co)
v = v + constituent%get_asFloat('v')
v = constituent%get_asFloat('v')
ph = phases%getIndex(constituent%get_asString('phase')) ph = phases%getIndex(constituent%get_asString('phase'))
do ip = 1, discretization_nIPs do ip = 1, discretization_nIPs
@ -152,10 +157,12 @@ subroutine parse()
material_phaseID(co,ce) = ph material_phaseID(co,ce) = ph
counterPhase(ph) = counterPhase(ph) + 1 counterPhase(ph) = counterPhase(ph) + 1
material_phaseEntry(co,ce) = counterPhase(ph) material_phaseEntry(co,ce) = counterPhase(ph)
material_v(co,ce) = v
end do end do
end do end do
if (dNeq(v,1.0_pReal,1.e-12_pReal)) call IO_error(153,ext_msg='constituent') if (dNeq(sum(material_v(1:constituents%length,ce)),1.0_pReal,1.e-9_pReal)) &
call IO_error(153,ext_msg='constituent')
end do end do

9
src/mesh/discretization_mesh.f90
View File

@ -49,9 +49,6 @@ module discretization_mesh
mesh_ipVolume, & !< volume associated with IP (initially!) mesh_ipVolume, & !< volume associated with IP (initially!)
mesh_node0 !< node x,y,z coordinates (initially!) mesh_node0 !< node x,y,z coordinates (initially!)
real(pReal), pointer, dimension(:) :: &
mesh_node0_temp
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!) mesh_ipCoordinates !< IP x,y,z coordinates (after deformation!)
@ -88,6 +85,8 @@ subroutine discretization_mesh_init(restart)
num_mesh num_mesh
integer :: p_i, dim !< integration order (quadrature rule) integer :: p_i, dim !< integration order (quadrature rule)
type(tvec) :: coords_node0 type(tvec) :: coords_node0
real(pReal), pointer, dimension(:) :: &
mesh_node0_temp
print'(/,1x,a)', '<<<+- discretization_mesh init -+>>>' print'(/,1x,a)', '<<<+- discretization_mesh init -+>>>'
@ -124,7 +123,7 @@ subroutine discretization_mesh_init(restart)
dimPlex = int(dim,pPETSCINT) dimPlex = int(dim,pPETSCINT)
if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error' if (err_MPI /= 0_MPI_INTEGER_KIND) error stop 'MPI error'
if (worldrank == 0) then if (worldsize == 1) then
call DMClone(globalMesh,geomMesh,err_PETSc) call DMClone(globalMesh,geomMesh,err_PETSc)
else else
call DMPlexDistribute(globalMesh,0_pPETSCINT,sf,geomMesh,err_PETSc) call DMPlexDistribute(globalMesh,0_pPETSCINT,sf,geomMesh,err_PETSc)
@ -155,7 +154,7 @@ subroutine discretization_mesh_init(restart)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
! Get initial nodal coordinates ! Get initial nodal coordinates
call DMGetCoordinates(geomMesh,coords_node0,err_PETSc) call DMGetCoordinatesLocal(geomMesh,coords_node0,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)
call VecGetArrayF90(coords_node0, mesh_node0_temp,err_PETSc) call VecGetArrayF90(coords_node0, mesh_node0_temp,err_PETSc)
CHKERRQ(err_PETSc) CHKERRQ(err_PETSc)

17
src/phase_mechanical_plastic.f90
View File

@ -140,15 +140,15 @@ submodule(phase:mechanical) plastic
dotState dotState
end function plastic_kinehardening_dotState end function plastic_kinehardening_dotState
module subroutine dislotwin_dotState(Mp,T,ph,en) module function dislotwin_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
end subroutine dislotwin_dotState real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
end function dislotwin_dotState
module function dislotungsten_dotState(Mp,ph,en) result(dotState) module function dislotungsten_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
@ -170,12 +170,10 @@ submodule(phase:mechanical) plastic
en en
end subroutine nonlocal_dotState end subroutine nonlocal_dotState
module subroutine dislotwin_dependentState(T,ph,en) module subroutine dislotwin_dependentState(ph,en)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
real(pReal), intent(in) :: &
T
end subroutine dislotwin_dependentState end subroutine dislotwin_dependentState
module subroutine dislotungsten_dependentState(ph,en) module subroutine dislotungsten_dependentState(ph,en)
@ -326,8 +324,7 @@ module function plastic_dotState(subdt,ph,en) result(dotState)
dotState = plastic_kinehardening_dotState(Mp,ph,en) dotState = plastic_kinehardening_dotState(Mp,ph,en)
case (PLASTIC_DISLOTWIN_ID) plasticType case (PLASTIC_DISLOTWIN_ID) plasticType
call dislotwin_dotState(Mp,thermal_T(ph,en),ph,en) dotState = dislotwin_dotState(Mp,ph,en)
dotState = plasticState(ph)%dotState(:,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
dotState = dislotungsten_dotState(Mp,ph,en) dotState = dislotungsten_dotState(Mp,ph,en)
@ -355,7 +352,7 @@ module subroutine plastic_dependentState(ph,en)
plasticType: select case (phase_plasticity(ph)) plasticType: select case (phase_plasticity(ph))
case (PLASTIC_DISLOTWIN_ID) plasticType case (PLASTIC_DISLOTWIN_ID) plasticType
call dislotwin_dependentState(thermal_T(ph,en),ph,en) call dislotwin_dependentState(ph,en)
case (PLASTIC_DISLOTUNGSTEN_ID) plasticType case (PLASTIC_DISLOTUNGSTEN_ID) plasticType
call dislotungsten_dependentState(ph,en) call dislotungsten_dependentState(ph,en)

336
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -17,22 +17,23 @@ submodule(phase:plastic) dislotwin
p_sb = 1.0_pReal, & !< p-exponent in shear band velocity p_sb = 1.0_pReal, & !< p-exponent in shear band velocity
q_sb = 1.0_pReal, & !< q-exponent in shear band velocity q_sb = 1.0_pReal, & !< q-exponent in shear band velocity
i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning i_tw = 1.0_pReal, & !< adjustment parameter to calculate MFP for twinning
L_tw = 1.0_pReal, & !< length of twin nuclei in Burgers vectors: TODO unit should be meters i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation
L_tr = 1.0_pReal, & !< length of trans nuclei in Burgers vectors: TODO unit should be meters L_tw = 1.0_pReal, & !< length of twin nuclei
x_c_tw = 1.0_pReal, & !< critical distance for formation of twin nucleus L_tr = 1.0_pReal, & !< length of trans nuclei
x_c_tr = 1.0_pReal, & !< critical distance for formation of trans nucleus x_c = 1.0_pReal, & !< critical distance for formation of twin/trans nucleus
V_cs = 1.0_pReal, & !< cross slip volume V_cs = 1.0_pReal, & !< cross slip volume
xi_sb = 1.0_pReal, & !< value for shearband resistance xi_sb = 1.0_pReal, & !< value for shearband resistance
v_sb = 1.0_pReal, & !< value for shearband velocity_0 v_sb = 1.0_pReal, & !< value for shearband velocity_0
E_sb = 1.0_pReal, & !< activation energy for shear bands E_sb = 1.0_pReal, & !< activation energy for shear bands
delta_G = 1.0_pReal, & !< free energy difference between austensite and martensite
i_tr = 1.0_pReal, & !< adjustment parameter to calculate MFP for transformation
h = 1.0_pReal, & !< stack height of hex nucleus h = 1.0_pReal, & !< stack height of hex nucleus
T_ref = T_ROOM, & gamma_char_tr = sqrt(0.125_pReal), & !< Characteristic shear for transformation
a_cI = 1.0_pReal, & a_cF = 1.0_pReal, &
a_cF = 1.0_pReal cOverA_hP = 1.0_pReal, &
real(pReal), dimension(2) :: & V_mol = 1.0_pReal, &
Gamma_sf = 0.0_pReal rho = 1.0_pReal
type(tPolynomial) :: &
Gamma_sf, & !< stacking fault energy
Delta_G !< free energy difference between austensite and martensite
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
b_sl, & !< absolute length of Burgers vector [m] for each slip system b_sl, & !< absolute length of Burgers vector [m] for each slip system
b_tw, & !< absolute length of Burgers vector [m] for each twin system b_tw, & !< absolute length of Burgers vector [m] for each twin system
@ -48,7 +49,7 @@ submodule(phase:plastic) dislotwin
r, & !< exponent in twin nucleation rate r, & !< exponent in twin nucleation rate
s, & !< exponent in trans nucleation rate s, & !< exponent in trans nucleation rate
tau_0, & !< strength due to elements in solid solution tau_0, & !< strength due to elements in solid solution
gamma_char, & !< characteristic shear for twins gamma_char_tw, & !< characteristic shear for twins
B, & !< drag coefficient B, & !< drag coefficient
d_caron !< distance of spontaneous annhihilation d_caron !< distance of spontaneous annhihilation
real(pReal), allocatable, dimension(:,:) :: & real(pReal), allocatable, dimension(:,:) :: &
@ -77,13 +78,22 @@ submodule(phase:plastic) dislotwin
character(len=pStringLen), allocatable, dimension(:) :: & character(len=pStringLen), allocatable, dimension(:) :: &
output output
logical :: & logical :: &
ExtendedDislocations, & !< consider split into partials for climb calculation extendedDislocations, & !< consider split into partials for climb calculation
fccTwinTransNucleation, & !< twinning and transformation models are for fcc fccTwinTransNucleation, & !< twinning and transformation models are for fcc
omitDipoles !< flag controlling consideration of dipole formation omitDipoles !< flag controlling consideration of dipole formation
character(len=:), allocatable, dimension(:) :: & character(len=:), allocatable, dimension(:) :: &
systems_sl, & systems_sl, &
systems_tw systems_tw
end type !< container type for internal constitutive parameters end type tParameters !< container type for internal constitutive parameters
type :: tIndexDotState
integer, dimension(2) :: &
rho_mob, &
rho_dip, &
gamma_sl, &
f_tw, &
f_tr
end type tIndexDotState
type :: tDislotwinState type :: tDislotwinState
real(pReal), dimension(:,:), pointer :: & real(pReal), dimension(:,:), pointer :: &
@ -99,17 +109,14 @@ submodule(phase:plastic) dislotwin
Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation
Lambda_tw, & !< mean free path between 2 obstacles seen by a growing twin Lambda_tw, & !< mean free path between 2 obstacles seen by a growing twin
Lambda_tr, & !< mean free path between 2 obstacles seen by a growing martensite Lambda_tr, & !< mean free path between 2 obstacles seen by a growing martensite
tau_pass, & !< threshold stress for slip tau_pass !< threshold stress for slip
tau_hat_tw, & !< threshold stress for twinning
tau_hat_tr !< threshold stress for transformation
end type tDislotwinDependentState end type tDislotwinDependentState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! containers for parameters and state ! containers for parameters and state
type(tParameters), allocatable, dimension(:) :: param type(tParameters), allocatable, dimension(:) :: param
type(tDislotwinState), allocatable, dimension(:) :: & type(tIndexDotState), allocatable, dimension(:) :: indexDotState
dotState, & type(tDislotwinState), allocatable, dimension(:) :: state
state
type(tDislotwinDependentState), allocatable, dimension(:) :: dependentState type(tDislotwinDependentState), allocatable, dimension(:) :: dependentState
contains contains
@ -129,6 +136,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
N_sl N_sl
real(pReal) :: a_cF
real(pReal), allocatable, dimension(:) :: & real(pReal), allocatable, dimension(:) :: &
rho_mob_0, & !< initial unipolar dislocation density per slip system rho_mob_0, & !< initial unipolar dislocation density per slip system
rho_dip_0 !< initial dipole dislocation density per slip system rho_dip_0 !< initial dipole dislocation density per slip system
@ -159,15 +167,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(phases%length)) allocate(param(phases%length))
allocate(indexDotState(phases%length))
allocate(state(phases%length)) allocate(state(phases%length))
allocate(dotState(phases%length))
allocate(dependentState(phases%length)) allocate(dependentState(phases%length))
do ph = 1, phases%length do ph = 1, phases%length
if (.not. myPlasticity(ph)) cycle if (.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
idx_dot => indexDotState(ph))
phase => phases%get(ph) phase => phases%get(ph)
mech => phase%get('mechanical') mech => phase%get('mechanical')
@ -179,7 +187,6 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray)
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! slip related parameters ! slip related parameters
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@ -209,7 +216,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%Q_cl = pl%get_asFloat('Q_cl') prm%Q_cl = pl%get_asFloat('Q_cl')
prm%ExtendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.) prm%extendedDislocations = pl%get_asBool('extend_dislocations',defaultVal = .false.)
prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.) prm%omitDipoles = pl%get_asBool('omit_dipoles',defaultVal = .false.)
! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
@ -250,7 +257,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! twin related parameters ! twin related parameters
prm%N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray) prm%N_tw = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
prm%sum_N_tw = sum(abs(prm%N_tw)) prm%sum_N_tw = sum(abs(prm%N_tw))
twinActive: if (prm%sum_N_tw > 0) then twinActive: if (prm%sum_N_tw > 0) then
prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph)) prm%systems_tw = lattice_labels_twin(prm%N_tw,phase_lattice(ph))
@ -262,16 +269,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%t_tw = pl%get_as1dFloat('t_tw', requiredSize=size(prm%N_tw)) prm%t_tw = pl%get_as1dFloat('t_tw', requiredSize=size(prm%N_tw))
prm%r = pl%get_as1dFloat('p_tw', requiredSize=size(prm%N_tw)) prm%r = pl%get_as1dFloat('p_tw', requiredSize=size(prm%N_tw))
prm%x_c_tw = pl%get_asFloat('x_c_tw')
prm%L_tw = pl%get_asFloat('L_tw') prm%L_tw = pl%get_asFloat('L_tw')
prm%i_tw = pl%get_asFloat('i_tw') prm%i_tw = pl%get_asFloat('i_tw')
prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph)) prm%gamma_char_tw = lattice_characteristicShear_Twin(prm%N_tw,phase_lattice(ph),phase_cOverA(ph))
if (.not. prm%fccTwinTransNucleation) then
prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw')
prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,prm%N_tw)
endif
! expand: family => system ! expand: family => system
prm%b_tw = math_expand(prm%b_tw,prm%N_tw) prm%b_tw = math_expand(prm%b_tw,prm%N_tw)
@ -279,54 +280,47 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%r = math_expand(prm%r,prm%N_tw) prm%r = math_expand(prm%r,prm%N_tw)
! sanity checks ! sanity checks
if ( prm%x_c_tw < 0.0_pReal) extmsg = trim(extmsg)//' x_c_tw' if (.not. prm%fccTwinTransNucleation) extmsg = trim(extmsg)//' TWIP for non-fcc'
if ( prm%L_tw < 0.0_pReal) extmsg = trim(extmsg)//' L_tw' if ( prm%L_tw < 0.0_pReal) extmsg = trim(extmsg)//' L_tw'
if ( prm%i_tw < 0.0_pReal) extmsg = trim(extmsg)//' i_tw' if ( prm%i_tw < 0.0_pReal) extmsg = trim(extmsg)//' i_tw'
if (any(prm%b_tw < 0.0_pReal)) extmsg = trim(extmsg)//' b_tw' if (any(prm%b_tw < 0.0_pReal)) extmsg = trim(extmsg)//' b_tw'
if (any(prm%t_tw < 0.0_pReal)) extmsg = trim(extmsg)//' t_tw' if (any(prm%t_tw < 0.0_pReal)) extmsg = trim(extmsg)//' t_tw'
if (any(prm%r < 0.0_pReal)) extmsg = trim(extmsg)//' p_tw' if (any(prm%r < 0.0_pReal)) extmsg = trim(extmsg)//' p_tw'
if (.not. prm%fccTwinTransNucleation) then
if (any(prm%dot_N_0_tw < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tw'
end if
else twinActive else twinActive
allocate(prm%gamma_char,prm%b_tw,prm%dot_N_0_tw,prm%t_tw,prm%r,source=emptyRealArray) allocate(prm%gamma_char_tw,prm%b_tw,prm%t_tw,prm%r,source=emptyRealArray)
allocate(prm%h_tw_tw(0,0)) allocate(prm%h_tw_tw(0,0))
end if twinActive end if twinActive
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! transformation related parameters ! transformation related parameters
prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray) prm%N_tr = pl%get_as1dInt('N_tr', defaultVal=emptyIntArray)
prm%sum_N_tr = sum(abs(prm%N_tr)) prm%sum_N_tr = sum(abs(prm%N_tr))
transActive: if (prm%sum_N_tr > 0) then transActive: if (prm%sum_N_tr > 0) then
prm%b_tr = pl%get_as1dFloat('b_tr') prm%b_tr = pl%get_as1dFloat('b_tr')
prm%b_tr = math_expand(prm%b_tr,prm%N_tr) prm%b_tr = math_expand(prm%b_tr,prm%N_tr)
prm%h = pl%get_asFloat('h', defaultVal=0.0_pReal) ! ToDo: This is not optional! prm%i_tr = pl%get_asFloat('i_tr')
prm%i_tr = pl%get_asFloat('i_tr', defaultVal=0.0_pReal) ! ToDo: This is not optional! prm%Delta_G = polynomial(pl%asDict(),'Delta_G','T')
prm%delta_G = pl%get_asFloat('delta_G') prm%L_tr = pl%get_asFloat('L_tr')
prm%x_c_tr = pl%get_asFloat('x_c_tr', defaultVal=0.0_pReal) ! ToDo: This is not optional! a_cF = pl%get_asFloat('a_cF')
prm%L_tr = pl%get_asFloat('L_tr') prm%h = 5.0_pReal * a_cF/sqrt(3.0_pReal)
prm%a_cI = pl%get_asFloat('a_cI', defaultVal=0.0_pReal) prm%cOverA_hP = pl%get_asFloat('c/a_hP')
prm%a_cF = pl%get_asFloat('a_cF', defaultVal=0.0_pReal) prm%rho = 4.0_pReal/(sqrt(3.0_pReal)*a_cF**2)/N_A
prm%V_mol = pl%get_asFloat('V_mol')
prm%lattice_tr = pl%get_asString('lattice_tr')
prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dFloat('h_tr-tr'),& prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,pl%get_as1dFloat('h_tr-tr'),&
phase_lattice(ph)) phase_lattice(ph))
prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,prm%lattice_tr, & prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr,'hP',prm%cOverA_hP)
0.0_pReal, &
prm%a_cI, &
prm%a_cF)
prm%t_tr = pl%get_as1dFloat('t_tr') prm%t_tr = pl%get_as1dFloat('t_tr')
prm%t_tr = math_expand(prm%t_tr,prm%N_tr) prm%t_tr = math_expand(prm%t_tr,prm%N_tr)
prm%s = pl%get_as1dFloat('p_tr',defaultVal=[0.0_pReal]) prm%s = pl%get_as1dFloat('p_tr')
prm%s = math_expand(prm%s,prm%N_tr) prm%s = math_expand(prm%s,prm%N_tr)
! sanity checks ! sanity checks
if ( prm%x_c_tr < 0.0_pReal) extmsg = trim(extmsg)//' x_c_tr' if (.not. prm%fccTwinTransNucleation) extmsg = trim(extmsg)//' TRIP for non-fcc'
if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr' if ( prm%L_tr < 0.0_pReal) extmsg = trim(extmsg)//' L_tr'
if ( prm%V_mol < 0.0_pReal) extmsg = trim(extmsg)//' V_mol'
if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr' if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr'
if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr' if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr'
if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr' if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr'
@ -356,15 +350,16 @@ module function plastic_dislotwin_init() result(myPlasticity)
if (prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) & if (prm%sum_N_sl + prm%sum_N_tw + prm%sum_N_tw > 0) &
prm%D = pl%get_asFloat('D') prm%D = pl%get_asFloat('D')
if (prm%sum_N_tw + prm%sum_N_tr > 0) & if (prm%sum_N_tw + prm%sum_N_tr > 0) then
prm%x_c = pl%get_asFloat('x_c')
prm%V_cs = pl%get_asFloat('V_cs') prm%V_cs = pl%get_asFloat('V_cs')
if (prm%x_c < 0.0_pReal) extmsg = trim(extmsg)//' x_c'
if (prm%sum_N_tw + prm%sum_N_tr > 0 .or. prm%ExtendedDislocations) then if (prm%V_cs < 0.0_pReal) extmsg = trim(extmsg)//' V_cs'
prm%T_ref = pl%get_asFloat('T_ref')
prm%Gamma_sf(1) = pl%get_asFloat('Gamma_sf')
prm%Gamma_sf(2) = pl%get_asFloat('Gamma_sf,T',defaultVal=0.0_pReal)
end if end if
if (prm%sum_N_tw + prm%sum_N_tr > 0 .or. prm%extendedDislocations) &
prm%Gamma_sf = polynomial(pl%asDict(),'Gamma_sf','T')
slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), & prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,prm%N_tw,pl%get_as1dFloat('h_sl-tw'), &
phase_lattice(ph)) phase_lattice(ph))
@ -385,55 +380,51 @@ module function plastic_dislotwin_init() result(myPlasticity)
+ size(['f_tr']) * prm%sum_N_tr + size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState sizeState = sizeDotState
call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0) call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,0)
deallocate(plasticState(ph)%dotState) ! ToDo: remove dotState completely
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
idx_dot%rho_mob = [startIndex,endIndex]
stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:) stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob= spread(rho_mob_0,2,Nmembers) stt%rho_mob= spread(rho_mob_0,2,Nmembers)
dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
idx_dot%rho_dip = [startIndex,endIndex]
stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:) stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip= spread(rho_dip_0,2,Nmembers) stt%rho_dip= spread(rho_dip_0,2,Nmembers)
dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
idx_dot%gamma_sl = [startIndex,endIndex]
stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:) stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
idx_dot%f_tw = [startIndex,endIndex]
stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:) stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tw',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tw'
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tr endIndex = endIndex + prm%sum_N_tr
idx_dot%f_tr = [startIndex,endIndex]
stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:) stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_f_tr',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr' if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_f_tr'
allocate(dst%Lambda_sl (prm%sum_N_sl,Nmembers),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl,Nmembers),source=0.0_pReal) allocate(dst%Lambda_sl(prm%sum_N_sl,Nmembers),source=0.0_pReal)
allocate(dst%Lambda_tw(prm%sum_N_tw,Nmembers),source=0.0_pReal)
allocate(dst%Lambda_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal) allocate(dst%Lambda_tr(prm%sum_N_tr,Nmembers),source=0.0_pReal)
allocate(dst%tau_hat_tw (prm%sum_N_tw,Nmembers),source=0.0_pReal)
allocate(dst%Lambda_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal)
allocate(dst%tau_hat_tr (prm%sum_N_tr,Nmembers),source=0.0_pReal)
end associate end associate
@ -482,7 +473,7 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
end if twinActive end if twinActive
transActive: if (prm%sum_N_tr > 0) then transActive: if (prm%sum_N_tr > 0) then
C66_tr = lattice_C66_trans(prm%N_tr,C,prm%lattice_tr,0.0_pReal,prm%a_cI,prm%a_cF) C66_tr = lattice_C66_trans(prm%N_tr,C,'hP',prm%cOverA_hP)
do i=1,prm%sum_N_tr do i=1,prm%sum_N_tr
homogenizedC = homogenizedC & homogenizedC = homogenizedC &
+ stt%f_tr(i,en)*C66_tr(1:6,1:6,i) + stt%f_tr(i,en)*C66_tr(1:6,1:6,i)
@ -559,7 +550,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
+ ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
end do slipContribution end do slipContribution
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw) if (prm%sum_N_tw > 0) call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
twinContibution: do i = 1, prm%sum_N_tw twinContibution: do i = 1, prm%sum_N_tw
Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i) Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -567,7 +558,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
+ ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
end do twinContibution end do twinContibution
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr) if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
transContibution: do i = 1, prm%sum_N_tr transContibution: do i = 1, prm%sum_N_tr
Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@ -581,7 +572,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
shearBandingContribution: if (dNeq0(prm%v_sb)) then shearBandingContribution: if (dNeq0(prm%v_sb)) then
E_kB_T = prm%E_sb/(K_B*T) E_kB_T = prm%E_sb/(K_B*T)
call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design? call math_eigh33(eigValues,eigVectors,Mp) ! is Mp symmetric by design?
do i = 1,6 do i = 1,6
P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
@ -612,15 +603,15 @@ end subroutine dislotwin_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate the rate of change of microstructure. !> @brief Calculate the rate of change of microstructure.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine dislotwin_dotState(Mp,T,ph,en) module function dislotwin_dotState(Mp,ph,en) result(dotState)
real(pReal), dimension(3,3), intent(in):: & real(pReal), dimension(3,3), intent(in):: &
Mp !< Mandel stress Mp !< Mandel stress
real(pReal), intent(in) :: &
T !< temperature at integration point
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
real(pReal), dimension(plasticState(ph)%sizeDotState) :: &
dotState
integer :: i integer :: i
real(pReal) :: & real(pReal) :: &
@ -639,22 +630,27 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
real(pReal), dimension(param(ph)%sum_N_tr) :: & real(pReal), dimension(param(ph)%sum_N_tr) :: &
dot_gamma_tr dot_gamma_tr
real(pReal) :: & real(pReal) :: &
mu, & mu, nu, &
nu, & T
Gamma
associate(prm => param(ph), stt => state(ph), dot => dotState(ph), dst => dependentState(ph))
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph), &
dot_rho_mob => dotState(indexDotState(ph)%rho_mob(1):indexDotState(ph)%rho_mob(2)), &
dot_rho_dip => dotState(indexDotState(ph)%rho_dip(1):indexDotState(ph)%rho_dip(2)), &
abs_dot_gamma_sl => dotState(indexDotState(ph)%gamma_sl(1):indexDotState(ph)%gamma_sl(2)), &
dot_f_tw => dotState(indexDotState(ph)%f_tw(1):indexDotState(ph)%f_tw(2)), &
dot_f_tr => dotState(indexDotState(ph)%f_tr(1):indexDotState(ph)%f_tr(2)))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en) nu = elastic_nu(ph,en)
Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref) T = thermal_T(ph,en)
f_matrix = 1.0_pReal & f_matrix = 1.0_pReal &
- sum(stt%f_tw(1:prm%sum_N_tw,en)) & - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
- sum(stt%f_tr(1:prm%sum_N_tr,en)) - sum(stt%f_tr(1:prm%sum_N_tr,en))
call kinetics_sl(Mp,T,ph,en,dot_gamma_sl) call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
dot%gamma_sl(:,en) = abs(dot_gamma_sl) abs_dot_gamma_sl = abs(dot_gamma_sl)
slipState: do i = 1, prm%sum_N_sl slipState: do i = 1, prm%sum_N_sl
tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i)) tau = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
@ -668,16 +664,18 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
d_hat = math_clip(d_hat, left = prm%d_caron(i)) d_hat = math_clip(d_hat, left = prm%d_caron(i))
dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) & dot_rho_dip_formation(i) = 2.0_pReal*(d_hat-prm%d_caron(i))/prm%b_sl(i) &
* stt%rho_mob(i,en)*abs(dot_gamma_sl(i)) * stt%rho_mob(i,en)*abs_dot_gamma_sl(i)
if (dEq(d_hat,prm%d_caron(i))) then if (dEq(d_hat,prm%d_caron(i))) then
dot_rho_dip_climb(i) = 0.0_pReal dot_rho_dip_climb(i) = 0.0_pReal
else else
! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 ! Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i)))
b_d = merge(24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) * Gamma / (mu*prm%b_sl(i)), & if (prm%extendedDislocations) then
1.0_pReal, & b_d = 24.0_pReal*PI*(1.0_pReal - nu)/(2.0_pReal + nu) * prm%Gamma_sf%at(T) / (mu*prm%b_sl(i))
prm%ExtendedDislocations) else
b_d = 1.0_pReal
end if
v_cl = 2.0_pReal*prm%omega*b_d**2*exp(-prm%Q_cl/(K_B*T)) & v_cl = 2.0_pReal*prm%omega*b_d**2*exp(-prm%Q_cl/(K_B*T)) &
* (exp(abs(sigma_cl)*prm%b_sl(i)**3/(K_B*T)) - 1.0_pReal) * (exp(abs(sigma_cl)*prm%b_sl(i)**3/(K_B*T)) - 1.0_pReal)
@ -687,38 +685,36 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
end if significantSlipStress end if significantSlipStress
end do slipState end do slipState
dot%rho_mob(:,en) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,en)) & dot_rho_mob = abs_dot_gamma_sl/(prm%b_sl*dst%Lambda_sl(:,en)) &
- dot_rho_dip_formation & - dot_rho_dip_formation &
- 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs(dot_gamma_sl) - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_mob(:,en)*abs_dot_gamma_sl
dot%rho_dip(:,en) = dot_rho_dip_formation & dot_rho_dip = dot_rho_dip_formation &
- 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) & - 2.0_pReal*prm%d_caron/prm%b_sl * stt%rho_dip(:,en)*abs_dot_gamma_sl &
- dot_rho_dip_climb - dot_rho_dip_climb
call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw) if (prm%sum_N_tw > 0) call kinetics_tw(Mp,T,abs_dot_gamma_sl,ph,en,dot_gamma_tw)
dot%f_tw(:,en) = f_matrix*dot_gamma_tw/prm%gamma_char dot_f_tw = f_matrix*dot_gamma_tw/prm%gamma_char_tw
call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr) if (prm%sum_N_tr > 0) call kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,dot_gamma_tr)
dot%f_tr(:,en) = f_matrix*dot_gamma_tr dot_f_tr = f_matrix*dot_gamma_tr/prm%gamma_char_tr
end associate end associate
end subroutine dislotwin_dotState end function dislotwin_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate derived quantities from state. !> @brief Calculate derived quantities from state.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine dislotwin_dependentState(T,ph,en) module subroutine dislotwin_dependentState(ph,en)
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
real(pReal), intent(in) :: &
T
real(pReal) :: & real(pReal) :: &
sumf_tw,Gamma,sumf_tr sumf_tw, sumf_tr
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
inv_lambda_sl inv_lambda_sl
real(pReal), dimension(param(ph)%sum_N_tw) :: & real(pReal), dimension(param(ph)%sum_N_tw) :: &
@ -728,19 +724,15 @@ module subroutine dislotwin_dependentState(T,ph,en)
inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite
f_over_t_tr f_over_t_tr
real(pReal) :: & real(pReal) :: &
mu, & mu
nu
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en)
sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en)) sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en)) sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
!* rescaled volume fraction for topology !* rescaled volume fraction for topology
f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ... f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw ! this is per system ...
f_over_t_tr = sumf_tr/prm%t_tr ! but this not f_over_t_tr = sumf_tr/prm%t_tr ! but this not
@ -762,13 +754,6 @@ module subroutine dislotwin_dependentState(T,ph,en)
!* threshold stress for dislocation motion !* threshold stress for dislocation motion
dst%tau_pass(:,en) = mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en))) dst%tau_pass(:,en) = mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
!* threshold stress for growing twin/martensite
dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
+ 3.0_pReal*prm%b_tw*mu/(prm%L_tw*prm%b_tw)
dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
+ 3.0_pReal*prm%b_tr*mu/(prm%L_tr*prm%b_tr) &
+ prm%h*prm%delta_G/(3.0_pReal*prm%b_tr)
end associate end associate
end subroutine dislotwin_dependentState end subroutine dislotwin_dependentState
@ -813,9 +798,6 @@ module subroutine plastic_dislotwin_results(ph,group)
case('Lambda_tw') case('Lambda_tw')
call results_writeDataset(dst%Lambda_tw,group,trim(prm%output(ou)), & call results_writeDataset(dst%Lambda_tw,group,trim(prm%output(ou)), &
'mean free path for twinning','m',prm%systems_tw) 'mean free path for twinning','m',prm%systems_tw)
case('tau_hat_tw')
call results_writeDataset(dst%tau_hat_tw,group,trim(prm%output(ou)), &
'threshold stress for twinning','Pa',prm%systems_tw)
case('f_tr') case('f_tr')
if (prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), & if (prm%sum_N_tr>0) call results_writeDataset(stt%f_tr,group,trim(prm%output(ou)), &
@ -847,15 +829,14 @@ pure subroutine kinetics_sl(Mp,T,ph,en, &
integer, intent(in) :: & integer, intent(in) :: &
ph, & ph, &
en en
real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: &
dot_gamma_sl dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: &
ddot_gamma_dtau_sl, & ddot_gamma_dtau_sl, &
tau_sl tau_sl
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
ddot_gamma_dtau ddot_gamma_dtau
real(pReal), dimension(param(ph)%sum_N_sl) :: & real(pReal), dimension(param(ph)%sum_N_sl) :: &
tau, & tau, &
stressRatio, & stressRatio, &
@ -914,7 +895,7 @@ end subroutine kinetics_sl
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,& pure subroutine kinetics_tw(Mp,T,abs_dot_gamma_sl,ph,en,&
dot_gamma_tw,ddot_gamma_dtau_tw) dot_gamma_tw,ddot_gamma_dtau_tw)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
@ -925,18 +906,17 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
ph, & ph, &
en en
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
dot_gamma_sl abs_dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
dot_gamma_tw dot_gamma_tw
real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tw), optional, intent(out) :: &
ddot_gamma_dtau_tw ddot_gamma_dtau_tw
real :: & real :: &
tau, tau_r, & tau, tau_r, tau_hat, &
dot_N_0, & dot_N_0, &
x0, V, & x0, V, &
Gamma, & Gamma_sf, &
mu, nu, & mu, nu, &
P_ncs, dP_ncs_dtau, & P_ncs, dP_ncs_dtau, &
P, dP_dtau P, dP_dtau
@ -947,50 +927,37 @@ pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
associate(prm => param(ph), stt => state(ph), dst => dependentState(ph)) associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
isFCC: if (prm%fccTwinTransNucleation) then mu = elastic_mu(ph,en)
mu = elastic_mu(ph,en) nu = elastic_nu(ph,en)
nu = elastic_nu(ph,en) Gamma_sf = prm%Gamma_sf%at(T)
Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref)
do i = 1, prm%sum_N_tw tau_hat = 3.0_pReal*prm%b_tw(1)*mu/prm%L_tw &
tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) + Gamma_sf/(3.0_pReal*prm%b_tw(1))
x0 = mu*prm%b_tw(i)**2*(2.0_pReal+nu)/(Gamma*8.0_pReal*PI*(1.0_pReal-nu)) ! ToDo: In the paper, the Burgers vector for slip is used x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu))
tau_r = mu*prm%b_tw(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(PI/3.0_pReal)/x0) ! ToDo: In the paper, the Burgers vector for slip is used tau_r = mu*prm%b_sl(1)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c)+cos(PI/3.0_pReal)/x0)
if (tau > tol_math_check .and. tau < tau_r) then do i = 1, prm%sum_N_tw
P = exp(-(dst%tau_hat_tw(i,en)/tau)**prm%r(i)) tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
dP_dTau = prm%r(i) * (dst%tau_hat_tw(i,en)/tau)**prm%r(i)/tau * P
s = prm%fcc_twinNucleationSlipPair(1:2,i) if (tau > tol_math_check .and. tau < tau_r) then
dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) & P = exp(-(tau_hat/tau)**prm%r(i))
/ (prm%L_tw*prm%b_sl(i)) dP_dTau = prm%r(i) * (tau_hat/tau)**prm%r(i)/tau * P
P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) s = prm%fcc_twinNucleationSlipPair(1:2,i)
dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) dot_N_0 = sum(abs_dot_gamma_sl(s(2:1:-1))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/(prm%L_tw*3.0_pReal)
V = PI/4.0_pReal*dst%Lambda_tw(i,en)**2*prm%t_tw(i) P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau))
dot_gamma_tw(i) = V*dot_N_0*P_ncs*P dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
if (present(ddot_gamma_dtau_tw)) &
ddot_gamma_dtau_tw(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)
else
dot_gamma_tw(i) = 0.0_pReal
if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal
end if
end do
else isFCC V = PI/4.0_pReal*dst%Lambda_tw(i,en)**2*prm%t_tw(i)
do i = 1, prm%sum_N_tw dot_gamma_tw(i) = V*dot_N_0*P_ncs*P*prm%gamma_char_tw(i)
error stop 'not implemented' if (present(ddot_gamma_dtau_tw)) &
tau = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) ddot_gamma_dtau_tw(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*prm%gamma_char_tw(i)
if (tau > tol_math_check) then else
dot_gamma_tw(i) = 0.0_pReal dot_gamma_tw(i) = 0.0_pReal
if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal
else end if
dot_gamma_tw(i) = 0.0_pReal end do
if (present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw(i) = 0.0_pReal
end if
end do
end if isFCC
end associate end associate
@ -1004,7 +971,7 @@ end subroutine kinetics_tw
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end. ! have the optional arguments at the end.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,& pure subroutine kinetics_tr(Mp,T,abs_dot_gamma_sl,ph,en,&
dot_gamma_tr,ddot_gamma_dtau_tr) dot_gamma_tr,ddot_gamma_dtau_tr)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
@ -1015,18 +982,17 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
ph, & ph, &
en en
real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: & real(pReal), dimension(param(ph)%sum_N_sl), intent(in) :: &
dot_gamma_sl abs_dot_gamma_sl
real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tr), intent(out) :: &
dot_gamma_tr dot_gamma_tr
real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: & real(pReal), dimension(param(ph)%sum_N_tr), optional, intent(out) :: &
ddot_gamma_dtau_tr ddot_gamma_dtau_tr
real :: & real :: &
tau, tau_r, & tau, tau_r, tau_hat, &
dot_N_0, & dot_N_0, &
x0, V, & x0, V, &
Gamma, & Gamma_sf, &
mu, nu, & mu, nu, &
P_ncs, dP_ncs_dtau, & P_ncs, dP_ncs_dtau, &
P, dP_dtau P, dP_dtau
@ -1039,28 +1005,30 @@ pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
mu = elastic_mu(ph,en) mu = elastic_mu(ph,en)
nu = elastic_nu(ph,en) nu = elastic_nu(ph,en)
Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * (T-prm%T_ref) Gamma_sf = prm%Gamma_sf%at(T)
tau_hat = 3.0_pReal*prm%b_tr(1)*mu/prm%L_tr &
+ (Gamma_sf + (prm%h/prm%V_mol - 2.0_pReal*prm%rho)*prm%Delta_G%at(T))/(3.0_pReal*prm%b_tr(1))
x0 = mu*prm%b_sl(1)**2*(2.0_pReal+nu)/(Gamma_sf*8.0_pReal*PI*(1.0_pReal-nu))
tau_r = mu*prm%b_sl(1)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c)+cos(PI/3.0_pReal)/x0)
do i = 1, prm%sum_N_tr do i = 1, prm%sum_N_tr
tau = math_tensordot(Mp,prm%P_tr(1:3,1:3,i)) tau = math_tensordot(Mp,prm%P_tr(1:3,1:3,i))
x0 = mu*prm%b_tr(i)**2*(2.0_pReal+nu)/(Gamma*8.0_pReal*PI*(1.0_pReal-nu)) ! ToDo: In the paper, the Burgers vector for slip is used
tau_r = mu*prm%b_tr(i)/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(PI/3.0_pReal)/x0) ! ToDo: In the paper, the Burgers vector for slip is used
if (tau > tol_math_check .and. tau < tau_r) then if (tau > tol_math_check .and. tau < tau_r) then
P = exp(-(dst%tau_hat_tr(i,en)/tau)**prm%s(i)) P = exp(-(tau_hat/tau)**prm%s(i))
dP_dTau = prm%s(i) * (dst%tau_hat_tr(i,en)/tau)**prm%s(i)/tau * P dP_dTau = prm%s(i) * (tau_hat/tau)**prm%s(i)/tau * P
s = prm%fcc_twinNucleationSlipPair(1:2,i) s = prm%fcc_twinNucleationSlipPair(1:2,i)
dot_N_0 = sum(abs(dot_gamma_sl(s(2:1:-1)))*(stt%rho_mob(s,en)+stt%rho_dip(s,en))) & dot_N_0 = sum(abs_dot_gamma_sl(s(2:1:-1))*(stt%rho_mob(s,en)+stt%rho_dip(s,en)))/(prm%L_tr*3.0_pReal)
/ (prm%L_tr*prm%b_sl(i))
P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau)) P_ncs = 1.0_pReal-exp(-prm%V_cs/(K_B*T)*(tau_r-tau))
dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal) dP_ncs_dtau = prm%V_cs / (K_B * T) * (P_ncs - 1.0_pReal)
V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i) V = PI/4.0_pReal*dst%Lambda_tr(i,en)**2*prm%t_tr(i)
dot_gamma_tr(i) = V*dot_N_0*P_ncs*P dot_gamma_tr(i) = V*dot_N_0*P_ncs*P*prm%gamma_char_tr
if (present(ddot_gamma_dtau_tr)) & if (present(ddot_gamma_dtau_tr)) &
ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau) ddot_gamma_dtau_tr(i) = V*dot_N_0*(P*dP_ncs_dtau + P_ncs*dP_dtau)*prm%gamma_char_tr
else else
dot_gamma_tr(i) = 0.0_pReal dot_gamma_tr(i) = 0.0_pReal
if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal if (present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr(i) = 0.0_pReal