Merge branch 'development' into 'python-vtk-improvements'

# Conflicts:
#   python/damask/_vtk.py
This commit is contained in:
Philip Eisenlohr 2022-02-23 14:20:52 +00:00
commit 320ab70270
44 changed files with 618 additions and 486 deletions

View File

@ -232,7 +232,9 @@ update_revision:
- git pull
- export VERSION=$(git describe ${CI_COMMIT_SHA})
- echo ${VERSION} > python/damask/VERSION
- git commit python/damask/VERSION -m "[skip ci] updated version information after successful test of $VERSION"
- >
git diff-index --quiet HEAD ||
git commit python/damask/VERSION -m "[skip ci] updated version information after successful test of $VERSION"
- if [ ${CI_COMMIT_SHA} == $(git rev-parse HEAD^) ]; then git push origin HEAD:master HEAD:development; fi
only:
- development

@ -1 +1 @@
Subproject commit e663a548b10ee3f9ced35c5a5105bd6824d3eebf
Subproject commit 4c8116ba3b9e9fbb325a580705028e8310139117

View File

@ -0,0 +1 @@
N_constituents: 8

View File

@ -1,4 +0,0 @@
[Parallel3]
mech isostrain
nconstituents 3
mapping sum # or 'parallel'

View File

@ -0,0 +1 @@
N_constituents: 2

View File

@ -0,0 +1,3 @@
# For single point calculations, requires N_constituents = 1
type: pass
output: ['T']

View File

@ -0,0 +1 @@
N_constituents: 1

View File

@ -1,10 +1,8 @@
8Grains:
N_constituents: 8
mechanical:
type: RGC
D_alpha: [4.0e-06, 4.0e-06, 2.0e-06]
a_g: [0.0, 0.0, 0.0]
c_alpha: 2.0
cluster_size: [2, 2, 2]
output: [M, Delta_V, avg_dot_a, max_dot_a]
xi_alpha: 10.0
# For Relaxed Grain Cluster homogenization, requires N_constituents = 8
type: RGC
D_alpha: [4.0e-06, 4.0e-06, 2.0e-06]
a_g: [0.0, 0.0, 0.0]
c_alpha: 2.0
cluster_size: [2, 2, 2]
output: [M, Delta_V, avg_dot_a, max_dot_a]
xi_alpha: 10.0

View File

@ -1,3 +0,0 @@
Taylor2:
N_constituents: 2
mechanical: {type: isostrain}

View File

@ -0,0 +1,3 @@
# For Taylor homogenization with N_constituents > 1
type: isostrain
output: ['F', 'P']

View File

@ -0,0 +1,3 @@
# For single point calculations, requires N_constituents = 1
type: pass
output: ['F', 'P']

View File

@ -0,0 +1,3 @@
# For homogenization with N_constituents > 1
type: isotemperature
output: ['T']

View File

@ -0,0 +1,3 @@
# For single point calculations, requires N_constituents = 1
type: pass
output: ['T']

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 3.4.1 (RRR=1000, T_min=200K, T_max=900K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
output: [T]
K_11: 2.380e+2
K_11,T: 2.068e-3
K_11,T^2: -7.765e-5

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 2.4.1 (RRR=1000, T_min=200K, T_max=1000K)
- https://www.mit.edu/~6.777/matprops/copper.htm
output: [T]
K_11: 4.039e+2
K_11,T: -8.119e-2
K_11,T^2: 1.454e-5

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 4.4.1 (RRR=300, T_min=200K, T_max=1000K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
output: [T]
K_11: 8.055e+1
K_11,T: -1.051e-1
K_11,T^2: 5.464e-5

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 35R (T_min=150K, T_max=500K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
output: [T]
K_11: 9.132e+1
K_11,T: -1.525e-1
K_11,T^2: 3.053e-4

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 61R (T_min=100K, T_max=400K)
- https://www.engineeringtoolbox.com/specific-heat-metals-d_152.html
output: [T]
K_11: 7.414e+1
K_11,T: -6.465e-2
K_11,T^2: 2.066e-4

View File

@ -5,6 +5,8 @@ references:
fit to Tab. 5.4.1 (RRR=300, T_min=200K, T_max=1000K)
- https://www.mit.edu/~6.777/matprops/tungsten.htm
output: [T]
K_11: 1.758e+2
K_11,T: -1.605e-1
K_11,T^2: 1.160e-4

View File

@ -1 +1 @@
v3.0.0-alpha6
v3.0.0-alpha6-23-gf944975ac

View File

@ -141,7 +141,7 @@ class Colormap(mpl.colors.ListedColormap):
)
if model.lower() not in toMsh:
raise ValueError(f'Invalid color model: {model}.')
raise ValueError(f'invalid color model "{model}"')
low_high = np.vstack((low,high)).astype(float)
out_of_bounds = np.bool_(False)
@ -156,7 +156,7 @@ class Colormap(mpl.colors.ListedColormap):
out_of_bounds = np.any(low_high[:,0]<0)
if out_of_bounds:
raise ValueError(f'{model.upper()} colors {low_high[0]} | {low_high[1]} are out of bounds.')
raise ValueError(f'{model.upper()} colors {low_high[0]} | {low_high[1]} are out of bounds')
low_,high_ = map(toMsh[model.lower()],low_high)
msh = map(functools.partial(Colormap._interpolate_msh,low=low_,high=high_),np.linspace(0,1,N))

View File

@ -29,7 +29,29 @@ lattice_symmetries: Dict[CrystalLattice, CrystalFamily] = {
class Crystal():
"""Crystal lattice."""
"""
Representation of a crystal as (general) crystal family or (more specific) as a scaled Bravais lattice.
Examples
--------
Cubic crystal family:
>>> import damask
>>> cubic = damask.Crystal(family='cubic')
>>> cubic
Crystal family: cubic
Body-centered cubic Bravais lattice with parameters of iron:
>>> import damask
>>> Fe = damask.Crystal(lattice='cI', a=287e-12)
>>> Fe
Crystal family: cubic
Bravais lattice: cI
a=2.87e-10 m, b=2.87e-10 m, c=2.87e-10 m
α=90°, β=90°, γ=90°
"""
def __init__(self, *,
family: CrystalFamily = None,
@ -38,7 +60,7 @@ class Crystal():
alpha: float = None, beta: float = None, gamma: float = None,
degrees: bool = False):
"""
Representation of crystal in terms of crystal family or Bravais lattice.
New representation of a crystal.
Parameters
----------
@ -114,10 +136,10 @@ class Crystal():
"""Give short human-readable summary."""
family = f'Crystal family: {self.family}'
return family if self.lattice is None else \
'\n'.join([family,
f'Bravais lattice: {self.lattice}',
'a={:.5g}m, b={:.5g}m, c={:.5g}m'.format(*self.parameters[:3]),
'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
util.srepr([family,
f'Bravais lattice: {self.lattice}',
'a={:.5g} m, b={:.5g} m, c={:.5g} m'.format(*self.parameters[:3]),
'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
def __eq__(self,
@ -243,6 +265,76 @@ class Crystal():
return _basis.get(self.family, None)
@property
def symmetry_operations(self) -> Rotation:
"""Symmetry operations as Rotations."""
_symmetry_operations: Dict[CrystalFamily, List] = {
'cubic': [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ],
[ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, -0.5 ],
[-0.5, 0.5, -0.5, 0.5 ],
[-0.5, -0.5, 0.5, 0.5 ],
[-0.5, -0.5, 0.5, -0.5 ],
[-0.5, -0.5, -0.5, 0.5 ],
[-0.5, 0.5, -0.5, -0.5 ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ],
[-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ],
],
'hexagonal': [
[ 1.0, 0.0, 0.0, 0.0 ],
[-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ],
[ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[ 0.0, 0.0, 0.0, 1.0 ],
[-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ],
[ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
],
'tetragonal': [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
],
'orthorhombic': [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
],
'monoclinic': [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
],
'triclinic': [
[ 1.0,0.0,0.0,0.0 ],
]}
return Rotation.from_quaternion(_symmetry_operations[self.family],accept_homomorph=True)
@property
def ratio(self):
"""Return axes ratios of own lattice."""
@ -323,7 +415,8 @@ class Crystal():
Parameters
----------
direction|plane : numpy.ndarray, shape (...,3)
Vector along direction or plane normal.
Real space vector along direction or
reciprocal space vector along plane normal.
Returns
-------
@ -344,7 +437,7 @@ class Crystal():
uvw: FloatSequence = None,
hkl: FloatSequence = None) -> np.ndarray:
"""
Calculate crystal frame vector along lattice direction [uvw] or plane normal (hkl).
Calculate crystal frame vector corresponding to lattice direction [uvw] or plane normal (hkl).
Parameters
----------
@ -354,7 +447,24 @@ class Crystal():
Returns
-------
vector : numpy.ndarray, shape (...,3)
Crystal frame vector along [uvw] direction or (hkl) plane normal.
Crystal frame vector in real space along [uvw] direction or
in reciprocal space along (hkl) plane normal.
Examples
--------
Crystal frame vector (real space) of Magnesium corresponding to [1,1,0] direction:
>>> import damask
>>> Mg = damask.Crystal(lattice='hP', a=321e-12, c=521e-12)
>>> Mg.to_frame(uvw=[1, 1, 0])
array([1.60500000e-10, 2.77994155e-10, 0.00000000e+00])
Crystal frame vector (reciprocal space) of Titanium along (1,0,0) plane normal:
>>> import damask
>>> Ti = damask.Crystal(lattice='hP', a=0.295e-9, c=0.469e-9)
>>> Ti.to_frame(hkl=(1, 0, 0))
array([ 3.38983051e+09, 1.95711956e+09, -4.15134508e-07])
"""
if (uvw is not None) ^ (hkl is None):

View File

@ -108,7 +108,7 @@ class Grid:
if len(material.shape) != 3:
raise ValueError(f'invalid material shape {material.shape}')
if material.dtype not in np.sctypes['float'] and material.dtype not in np.sctypes['int']:
raise TypeError(f'invalid material data type {material.dtype}')
raise TypeError(f'invalid material data type "{material.dtype}"')
self._material = np.copy(material)
@ -229,7 +229,7 @@ class Grid:
except ValueError:
header_length,keyword = (-1, 'invalid')
if not keyword.startswith('head') or header_length < 3:
raise TypeError('header length information missing or invalid')
raise TypeError('invalid or missing header length information')
comments = []
content = f.readlines()
@ -259,7 +259,7 @@ class Grid:
i += len(items)
if i != cells.prod():
raise TypeError(f'invalid file: expected {cells.prod()} entries, found {i}')
raise TypeError(f'mismatch between {cells.prod()} expected entries and {i} found')
if not np.any(np.mod(material,1) != 0.0): # no float present
material = material.astype('int') - (1 if material.min() > 0 else 0)
@ -823,7 +823,7 @@ class Grid:
"""
if not set(directions).issubset(valid := ['x', 'y', 'z']):
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
raise ValueError(f'invalid direction "{set(directions).difference(valid)}" specified')
limits: Sequence[Optional[int]] = [None,None] if reflect else [-2,0]
mat = self.material.copy()
@ -859,7 +859,7 @@ class Grid:
"""
if not set(directions).issubset(valid := ['x', 'y', 'z']):
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
raise ValueError(f'invalid direction "{set(directions).difference(valid)}" specified')
mat = np.flip(self.material, [valid.index(d) for d in directions if d in valid])
@ -1196,7 +1196,7 @@ class Grid:
"""
if not set(directions).issubset(valid := ['x', 'y', 'z']):
raise ValueError(f'invalid direction {set(directions).difference(valid)} specified')
raise ValueError(f'invalid direction "{set(directions).difference(valid)}" specified')
o = [[0, self.cells[0]+1, np.prod(self.cells[:2]+1)+self.cells[0]+1, np.prod(self.cells[:2]+1)],
[0, np.prod(self.cells[:2]+1), np.prod(self.cells[:2]+1)+1, 1],

View File

@ -1,6 +1,6 @@
import inspect
import copy
from typing import Union, Callable, List, Dict, Any, Tuple, TypeVar
from typing import Union, Callable, Dict, Any, Tuple, TypeVar
import numpy as np
@ -780,76 +780,6 @@ class Orientation(Rotation,Crystal):
return rgb
@property
def symmetry_operations(self) -> Rotation:
"""Symmetry operations as Rotations."""
_symmetry_operations: Dict[CrystalFamily, List] = {
'cubic': [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ],
[ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, -0.5 ],
[-0.5, 0.5, -0.5, 0.5 ],
[-0.5, -0.5, 0.5, 0.5 ],
[-0.5, -0.5, 0.5, -0.5 ],
[-0.5, -0.5, -0.5, 0.5 ],
[-0.5, 0.5, -0.5, -0.5 ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ],
[-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ],
],
'hexagonal': [
[ 1.0, 0.0, 0.0, 0.0 ],
[-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ],
[ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[ 0.0, 0.0, 0.0, 1.0 ],
[-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ],
[ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
],
'tetragonal': [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
],
'orthorhombic': [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
],
'monoclinic': [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
],
'triclinic': [
[ 1.0,0.0,0.0,0.0 ],
]}
return Rotation.from_quaternion(_symmetry_operations[self.family],accept_homomorph=True)
####################################################################################################
# functions that require lattice, not just family
@ -919,7 +849,7 @@ class Orientation(Rotation,Crystal):
"""
if (N_slip is not None) ^ (N_twin is None):
raise KeyError('Specify either "N_slip" or "N_twin"')
raise KeyError('specify either "N_slip" or "N_twin"')
kinematics,active = (self.kinematics('slip'),N_slip) if N_twin is None else \
(self.kinematics('twin'),N_twin)

View File

@ -32,10 +32,10 @@ def _view_transition(what,datasets,increments,times,phases,homogenizations,field
if (datasets is not None and what is None) or (what is not None and datasets is None):
raise ValueError('"what" and "datasets" need to be used as a pair')
if datasets is not None or what is not None:
warnings.warn('Arguments "what" and "datasets" will be removed in DAMASK v3.0.0-alpha7', DeprecationWarning,2)
warnings.warn('arguments "what" and "datasets" will be removed in DAMASK v3.0.0-alpha7', DeprecationWarning,2)
return what,datasets
if sum(1 for _ in filter(None.__ne__, [increments,times,phases,homogenizations,fields])) > 1:
raise ValueError('Only one out of "increments", "times", "phases", "homogenizations", and "fields" can be used')
raise ValueError('only one out of "increments", "times", "phases", "homogenizations", and "fields" can be used')
else:
if increments is not None: return "increments", increments
if times is not None: return "times", times
@ -115,7 +115,7 @@ class Result:
self.version_minor = f.attrs['DADF5_version_minor']
if self.version_major != 0 or not 12 <= self.version_minor <= 14:
raise TypeError(f'Unsupported DADF5 version {self.version_major}.{self.version_minor}')
raise TypeError(f'unsupported DADF5 version "{self.version_major}.{self.version_minor}"')
if self.version_major == 0 and self.version_minor < 14:
self.export_setup = None
@ -490,7 +490,7 @@ class Result:
"""
if self._protected:
raise PermissionError('Renaming datasets not permitted')
raise PermissionError('rename datasets')
with h5py.File(self.fname,'a') as f:
for inc in self.visible['increments']:
@ -529,7 +529,7 @@ class Result:
"""
if self._protected:
raise PermissionError('Removing datasets not permitted')
raise PermissionError('delete datasets')
with h5py.File(self.fname,'a') as f:
for inc in self.visible['increments']:
@ -639,7 +639,7 @@ class Result:
data = eval(formula)
if not hasattr(data,'shape') or data.shape[0] != kwargs[d]['data'].shape[0]:
raise ValueError("'{}' results in invalid shape".format(kwargs['formula']))
raise ValueError('"{}" results in invalid shape'.format(kwargs['formula']))
return {
'data': data,
@ -939,7 +939,7 @@ class Result:
elif T_sym['meta']['unit'] == 'Pa':
k = 'stress'
if k not in ['stress', 'strain']:
raise ValueError(f'invalid von Mises kind {kind}')
raise ValueError(f'invalid von Mises kind "{kind}"')
return {
'data': (mechanics.equivalent_strain_Mises if k=='strain' else \
@ -1633,7 +1633,7 @@ class Result:
elif mode.lower()=='point':
v = VTK.from_poly_data(self.coordinates0_point)
else:
raise ValueError(f'invalid mode {mode}')
raise ValueError(f'invalid mode "{mode}"')
v.comments = util.execution_stamp('Result','export_VTK')

View File

@ -279,7 +279,7 @@ class Rotation:
p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o)
return self.copy(Rotation(np.block([q,p]))._standardize())
else:
raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
raise TypeError('use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
def __imul__(self: MyType,
other: MyType) -> MyType:
@ -314,7 +314,7 @@ class Rotation:
if isinstance(other,Rotation):
return self*~other
else:
raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
raise TypeError('use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
def __itruediv__(self: MyType,
other: MyType) -> MyType:
@ -365,11 +365,11 @@ class Rotation:
R = self.as_matrix()
return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other)
else:
raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors')
raise ValueError('can only rotate vectors, 2nd order tensors, and 4th order tensors')
elif isinstance(other, Rotation):
raise TypeError('Use "R1*R2", i.e. multiplication, to compose rotations "R1" and "R2"')
raise TypeError('use "R1*R2", i.e. multiplication, to compose rotations "R1" and "R2"')
else:
raise TypeError(f'Cannot rotate {type(other)}')
raise TypeError(f'cannot rotate "{type(other)}"')
apply = __matmul__
@ -731,7 +731,7 @@ class Rotation:
"""
qu = np.array(q,dtype=float)
if qu.shape[:-2:-1] != (4,):
raise ValueError('Invalid shape.')
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
@ -740,9 +740,9 @@ class Rotation:
qu[qu[...,0] < 0.0] *= -1
else:
if np.any(qu[...,0] < 0.0):
raise ValueError('Quaternion with negative first (real) component.')
raise ValueError('quaternion with negative first (real) component')
if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0,rtol=0.0)):
raise ValueError('Quaternion is not of unit length.')
raise ValueError('quaternion is not of unit length')
return Rotation(qu)
@ -767,11 +767,11 @@ class Rotation:
"""
eu = np.array(phi,dtype=float)
if eu.shape[:-2:-1] != (3,):
raise ValueError('Invalid shape.')
raise ValueError('invalid shape')
eu = np.radians(eu) if degrees else eu
if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or np.any(eu[...,1] > np.pi): # ToDo: No separate check for PHI
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].')
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]')
return Rotation(Rotation._eu2qu(eu))
@ -798,7 +798,7 @@ class Rotation:
"""
ax = np.array(axis_angle,dtype=float)
if ax.shape[:-2:-1] != (4,):
raise ValueError('Invalid shape.')
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
@ -806,10 +806,10 @@ class Rotation:
if degrees: ax[..., 3] = np.radians(ax[...,3])
if normalize: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1,keepdims=True)
if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi):
raise ValueError('Axisangle rotation angle outside of [0..π].')
raise ValueError('axisangle rotation angle outside of [0..π]')
if not np.all(np.isclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0)):
print(np.linalg.norm(ax[...,0:3],axis=-1))
raise ValueError('Axisangle rotation axis is not of unit length.')
raise ValueError('axisangle rotation axis is not of unit length')
return Rotation(Rotation._ax2qu(ax))
@ -832,7 +832,7 @@ class Rotation:
"""
om = np.array(basis,dtype=float)
if om.shape[-2:] != (3,3):
raise ValueError('Invalid shape.')
raise ValueError('invalid shape')
if reciprocal:
om = np.linalg.inv(tensor.transpose(om)/np.pi) # transform reciprocal basis set
@ -841,11 +841,11 @@ class Rotation:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.einsum('...ij,...jl',U,Vh)
if not np.all(np.isclose(np.linalg.det(om),1.0)):
raise ValueError('Orientation matrix has determinant ≠ 1.')
raise ValueError('orientation matrix has determinant ≠ 1')
if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \
or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \
or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)):
raise ValueError('Orientation matrix is not orthogonal.')
raise ValueError('orientation matrix is not orthogonal')
return Rotation(Rotation._om2qu(om))
@ -879,7 +879,7 @@ class Rotation:
a_ = np.array(a)
b_ = np.array(b)
if a_.shape[-2:] != (2,3) or b_.shape[-2:] != (2,3) or a_.shape != b_.shape:
raise ValueError('Invalid shape.')
raise ValueError('invalid shape')
am = np.stack([ a_[...,0,:],
a_[...,1,:],
np.cross(a_[...,0,:],a_[...,1,:]) ],axis=-2)
@ -910,16 +910,16 @@ class Rotation:
"""
ro = np.array(rho,dtype=float)
if ro.shape[:-2:-1] != (4,):
raise ValueError('Invalid shape.')
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
ro[...,0:3] *= -P
if normalize: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True)
if np.any(ro[...,3] < 0.0):
raise ValueError('Rodrigues vector rotation angle is negative.')
raise ValueError('Rodrigues vector rotation angle is negative')
if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)):
raise ValueError('Rodrigues vector rotation axis is not of unit length.')
raise ValueError('Rodrigues vector rotation axis is not of unit length')
return Rotation(Rotation._ro2qu(ro))
@ -939,14 +939,14 @@ class Rotation:
"""
ho = np.array(h,dtype=float)
if ho.shape[:-2:-1] != (3,):
raise ValueError('Invalid shape.')
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
ho *= -P
if np.any(np.linalg.norm(ho,axis=-1) >_R1+1e-9):
raise ValueError('Homochoric coordinate outside of the sphere.')
raise ValueError('homochoric coordinate outside of the sphere')
return Rotation(Rotation._ho2qu(ho))
@ -966,12 +966,12 @@ class Rotation:
"""
cu = np.array(x,dtype=float)
if cu.shape[:-2:-1] != (3,):
raise ValueError('Invalid shape.')
raise ValueError('invalid shape')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
if np.abs(np.max(cu)) > np.pi**(2./3.) * 0.5+1e-9:
raise ValueError('Cubochoric coordinate outside of the cube.')
raise ValueError('cubochoric coordinate outside of the cube')
ho = -P * Rotation._cu2ho(cu)

View File

@ -532,7 +532,7 @@ class Table:
"""
if self.shapes != other.shapes or not self.data.columns.equals(other.data.columns):
raise KeyError('Labels or shapes or order do not match')
raise KeyError('mismatch of shapes or labels or their order')
dup = self.copy()
dup.data = dup.data.append(other.data,ignore_index=True)
@ -558,7 +558,7 @@ class Table:
"""
if set(self.shapes) & set(other.shapes) or self.data.shape[0] != other.data.shape[0]:
raise KeyError('Duplicated keys or row count mismatch')
raise KeyError('duplicated keys or row count mismatch')
dup = self.copy()
dup.data = dup.data.join(other.data)

View File

@ -421,8 +421,7 @@ class VTK:
data: np.ndarray,
label: str):
N_data = data.shape[0]
N_data = data.shape[0]
data_ = data.reshape(N_data,-1) \
.astype(np.single if data.dtype in [np.double,np.longdouble] else data.dtype)

View File

@ -300,7 +300,7 @@ def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
origin[_np.where(cells==1)] = 0.0
if cells.prod() != len(coordinates0):
raise ValueError(f'Data count {len(coordinates0)} does not match cells {cells}.')
raise ValueError(f'data count {len(coordinates0)} does not match cells {cells}')
start = origin + delta*.5
end = origin - delta*.5 + size
@ -309,11 +309,11 @@ def cellsSizeOrigin_coordinates0_point(coordinates0: _np.ndarray,
if not (_np.allclose(coords[0],_np.linspace(start[0],end[0],cells[0]),atol=atol) and \
_np.allclose(coords[1],_np.linspace(start[1],end[1],cells[1]),atol=atol) and \
_np.allclose(coords[2],_np.linspace(start[2],end[2],cells[2]),atol=atol)):
raise ValueError('Regular cell spacing violated.')
raise ValueError('non-uniform cell spacing')
if ordered and not _np.allclose(coordinates0.reshape(tuple(cells)+(3,),order='F'),
coordinates0_point(list(cells),size,origin),atol=atol):
raise ValueError('Input data is not ordered (x fast, z slow).')
raise ValueError('input data is not ordered (x fast, z slow)')
return (cells,size,origin)
@ -460,17 +460,17 @@ def cellsSizeOrigin_coordinates0_node(coordinates0: _np.ndarray,
origin = mincorner
if (cells+1).prod() != len(coordinates0):
raise ValueError(f'Data count {len(coordinates0)} does not match cells {cells}.')
raise ValueError(f'data count {len(coordinates0)} does not match cells {cells}')
atol = _np.max(size)*5e-2
if not (_np.allclose(coords[0],_np.linspace(mincorner[0],maxcorner[0],cells[0]+1),atol=atol) and \
_np.allclose(coords[1],_np.linspace(mincorner[1],maxcorner[1],cells[1]+1),atol=atol) and \
_np.allclose(coords[2],_np.linspace(mincorner[2],maxcorner[2],cells[2]+1),atol=atol)):
raise ValueError('Regular cell spacing violated.')
raise ValueError('non-uniform cell spacing')
if ordered and not _np.allclose(coordinates0.reshape(tuple(cells+1)+(3,),order='F'),
coordinates0_node(list(cells),size,origin),atol=atol):
raise ValueError('Input data is not ordered (x fast, z slow).')
raise ValueError('input data is not ordered (x fast, z slow)')
return (cells,size,origin)

View File

@ -258,7 +258,7 @@ def _polar_decomposition(T: _np.ndarray,
Tensor of which the singular values are computed.
requested : sequence of {'R', 'U', 'V'}
Requested outputs: R for the rotation tensor,
V for left stretch tensor and U for right stretch tensor.
V for left stretch tensor, and U for right stretch tensor.
"""
u, _, vh = _np.linalg.svd(T)
@ -273,7 +273,7 @@ def _polar_decomposition(T: _np.ndarray,
output+=[_np.einsum('...ji,...jk',R,T)]
if len(output) == 0:
raise ValueError('output needs to be out of V, R, U')
raise ValueError('output not in {V, R, U}')
return tuple(output)

View File

@ -99,8 +99,8 @@ def from_Poisson_disc(size: _FloatSequence,
s += 1
progress.update(s)
if i == 100:
raise ValueError('Seeding not possible')
if i >= 100:
raise ValueError('seeding not possible')
return coords

View File

@ -285,7 +285,7 @@ def scale_to_coprime(v: FloatSequence) -> np.ndarray:
with np.errstate(invalid='ignore'):
if not np.allclose(np.ma.masked_invalid(v_/m),v_[np.argmax(abs(v_))]/m[np.argmax(abs(v_))]):
raise ValueError(f'Invalid result {m} for input {v_}. Insufficient precision?')
raise ValueError(f'invalid result "{m}" for input "{v_}"')
return m
@ -483,7 +483,7 @@ def shapeshifter(fro: Tuple[int, ...],
assert match
grp = match.groups()
except AssertionError:
raise ValueError(f'Shapes can not be shifted {fro} --> {to}')
raise ValueError(f'shapes cannot be shifted {fro} --> {to}')
fill: Any = ()
for g,d in zip(grp,fro+(None,)):
fill += (1,)*g.count(',')+(d,)
@ -576,7 +576,7 @@ def DREAM3D_base_group(fname: Union[str, Path]) -> str:
base_group = f.visit(lambda path: path.rsplit('/',2)[0] if '_SIMPL_GEOMETRY/SPACING' in path else None)
if base_group is None:
raise ValueError(f'Could not determine base group in file {fname}.')
raise ValueError(f'could not determine base group in file "{fname}"')
return base_group
@ -607,7 +607,7 @@ def DREAM3D_cell_data_group(fname: Union[str, Path]) -> str:
else None)
if cell_data_group is None:
raise ValueError(f'Could not determine cell data group in file {fname}/{base_group}.')
raise ValueError(f'could not determine cell-data group in file "{fname}/{base_group}"')
return cell_data_group
@ -630,7 +630,7 @@ def Bravais_to_Miller(*,
"""
if (uvtw is not None) ^ (hkil is None):
raise KeyError('Specify either "uvtw" or "hkil"')
raise KeyError('specify either "uvtw" or "hkil"')
axis,basis = (np.array(uvtw),np.array([[1,0,-1,0],
[0,1,-1,0],
[0,0, 0,1]])) \
@ -659,7 +659,7 @@ def Miller_to_Bravais(*,
"""
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),np.array([[ 2,-1, 0],
[-1, 2, 0],
[-1,-1, 0],

View File

@ -6,9 +6,10 @@
!--------------------------------------------------------------------------------------------------
module homogenization
use prec
use math
use constants
use IO
use config
use math
use material
use phase
use discretization
@ -32,7 +33,7 @@ module homogenization
HOMOGENIZATION_RGC_ID
end enum
type(tState), allocatable, dimension(:), public :: &
type(tState), allocatable, dimension(:), public :: &
homogState, &
damageState_h
@ -66,15 +67,13 @@ module homogenization
!--------------------------------------------------------------------------------------------------
interface
module subroutine mechanical_init(num_homog)
class(tNode), pointer, intent(in) :: &
num_homog !< pointer to mechanical homogenization numerics data
module subroutine mechanical_init()
end subroutine mechanical_init
module subroutine thermal_init
module subroutine thermal_init()
end subroutine thermal_init
module subroutine damage_init
module subroutine damage_init()
end subroutine damage_init
module subroutine mechanical_partition(subF,ce)
@ -204,15 +203,15 @@ subroutine homogenization_init()
allocate(homogState (size(material_name_homogenization)))
allocate(damageState_h (size(material_name_homogenization)))
call material_parseHomogenization()
call parseHomogenization()
num_homog => config_numerics%get('homogenization',defaultVal=emptyDict)
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
num%nMPstate = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
num%nMPstate = num_homogGeneric%get_asInt('nMPstate',defaultVal=10)
if (num%nMPstate < 1) call IO_error(301,ext_msg='nMPstate')
call mechanical_init(num_homog)
call mechanical_init()
call thermal_init()
call damage_init()
@ -323,13 +322,13 @@ subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolvin
elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip
ho = material_homogenizationID(ce)
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
enddo
call mechanical_homogenize(Delta_t,ce)
enddo IpLooping3
enddo elementLooping3
ho = material_homogenizationID(ce)
do co = 1, homogenization_Nconstituents(ho)
call crystallite_orientations(co,ip,el)
end do
call mechanical_homogenize(Delta_t,ce)
end do IpLooping3
end do elementLooping3
!$OMP END PARALLEL DO
@ -447,7 +446,7 @@ end subroutine homogenization_restartRead
!--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration
!--------------------------------------------------------------------------------------------------
subroutine material_parseHomogenization
subroutine parseHomogenization
class(tNode), pointer :: &
material_homogenization, &
@ -459,8 +458,8 @@ subroutine material_parseHomogenization
material_homogenization => config_material%get('homogenization')
allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID)
allocate(thermal_type(size(material_name_homogenization)),source=THERMAL_isothermal_ID)
allocate(damage_type (size(material_name_homogenization)),source=DAMAGE_none_ID)
do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h)
@ -468,7 +467,7 @@ subroutine material_parseHomogenization
if (homog%contains('thermal')) then
homogThermal => homog%get('thermal')
select case (homogThermal%get_asString('type'))
case('pass')
case('pass','isotemperature')
thermal_type(h) = THERMAL_conduction_ID
case default
call IO_error(500,ext_msg=homogThermal%get_asString('type'))
@ -486,7 +485,7 @@ subroutine material_parseHomogenization
endif
enddo
end subroutine material_parseHomogenization
end subroutine parseHomogenization
end module homogenization

View File

@ -7,15 +7,13 @@ submodule(homogenization) mechanical
interface
module subroutine pass_init
module subroutine pass_init()
end subroutine pass_init
module subroutine isostrain_init
module subroutine isostrain_init()
end subroutine isostrain_init
module subroutine RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data
module subroutine RGC_init()
end subroutine RGC_init
@ -52,6 +50,12 @@ submodule(homogenization) mechanical
end interface
type :: tOutput !< requested output (per phase)
character(len=pStringLen), allocatable, dimension(:) :: &
label
end type tOutput
type(tOutput), allocatable, dimension(:) :: output_mechanical
integer(kind(HOMOGENIZATION_undefined_ID)), dimension(:), allocatable :: &
homogenization_type !< type of each homogenization
@ -60,27 +64,20 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters.
!--------------------------------------------------------------------------------------------------
module subroutine mechanical_init(num_homog)
class(tNode), pointer, intent(in) :: &
num_homog
class(tNode), pointer :: &
num_homogMech
module subroutine mechanical_init()
print'(/,1x,a)', '<<<+- homogenization:mechanical init -+>>>'
call material_parseHomogenization2()
call parseMechanical()
allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity
homogenization_F = homogenization_F0 ! initialize to identity
allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
allocate(homogenization_dPdF(3,3,3,3,discretization_Ncells), source=0.0_pReal)
homogenization_F0 = spread(math_I3,3,discretization_Ncells)
homogenization_F = homogenization_F0
allocate(homogenization_P(3,3,discretization_Ncells),source=0.0_pReal)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call pass_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call RGC_init(num_homogMech)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call pass_init()
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init()
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call RGC_init()
end subroutine mechanical_init
@ -185,8 +182,10 @@ module subroutine mechanical_results(group_base,ho)
character(len=*), intent(in) :: group_base
integer, intent(in) :: ho
integer :: ou
character(len=:), allocatable :: group
group = trim(group_base)//'/mechanical'
call results_closeGroup(results_addGroup(group))
@ -197,12 +196,17 @@ module subroutine mechanical_results(group_base,ho)
end select
!temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'F',&
! 'deformation gradient','1')
!temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchhoff stress','Pa')
do ou = 1, size(output_mechanical(1)%label)
select case (output_mechanical(ho)%label(ou))
case('F')
call results_writeDataset(reshape(homogenization_F,[3,3,discretization_nCells]),group,'F', &
'deformation gradient','1')
case('P')
call results_writeDataset(reshape(homogenization_P,[3,3,discretization_nCells]),group,'P', &
'deformation gradient','1')
end select
end do
end subroutine mechanical_results
@ -210,35 +214,42 @@ end subroutine mechanical_results
!--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration
!--------------------------------------------------------------------------------------------------
subroutine material_parseHomogenization2()
subroutine parseMechanical()
class(tNode), pointer :: &
material_homogenization, &
homog, &
homogMech
mechanical
integer :: ho
integer :: h
material_homogenization => config_material%get('homogenization')
allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(output_mechanical(size(material_name_homogenization)))
do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h)
homogMech => homog%get('mechanical')
select case (homogMech%get_asString('type'))
do ho=1, size(material_name_homogenization)
homog => material_homogenization%get(ho)
mechanical => homog%get('mechanical')
#if defined(__GFORTRAN__)
output_mechanical(ho)%label = output_as1dString(mechanical)
#else
output_mechanical(ho)%label = mechanical%get_as1dString('output',defaultVal=emptyStringArray)
#endif
select case (mechanical%get_asString('type'))
case('pass')
homogenization_type(h) = HOMOGENIZATION_NONE_ID
homogenization_type(ho) = HOMOGENIZATION_NONE_ID
case('isostrain')
homogenization_type(h) = HOMOGENIZATION_ISOSTRAIN_ID
homogenization_type(ho) = HOMOGENIZATION_ISOSTRAIN_ID
case('RGC')
homogenization_type(h) = HOMOGENIZATION_RGC_ID
homogenization_type(ho) = HOMOGENIZATION_RGC_ID
case default
call IO_error(500,ext_msg=homogMech%get_asString('type'))
call IO_error(500,ext_msg=mechanical%get_asString('type'))
end select
end do
end subroutine material_parseHomogenization2
end subroutine parseMechanical
end submodule mechanical

View File

@ -71,10 +71,7 @@ contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all necessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
module subroutine RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data
module subroutine RGC_init()
integer :: &
ho, &
@ -82,6 +79,8 @@ module subroutine RGC_init(num_homogMech)
sizeState, nIntFaceTot
class (tNode), pointer :: &
num_homogenization, &
num_mechanical, &
num_RGC, & ! pointer to RGC numerics data
material_homogenization, &
homog, &
@ -105,7 +104,9 @@ module subroutine RGC_init(num_homogMech)
allocate(state0(material_homogenization%length))
allocate(dependentState(material_homogenization%length))
num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)
num_homogenization => config_numerics%get('homogenization',defaultVal=emptyDict)
num_mechanical => num_homogenization%get('mechanical',defaultVal=emptyDict)
num_RGC => num_mechanical%get('RGC',defaultVal=emptyDict)
num%atol = num_RGC%get_asFloat('atol', defaultVal=1.0e+4_pReal)
num%rtol = num_RGC%get_asFloat('rtol', defaultVal=1.0e-3_pReal)
@ -171,8 +172,8 @@ module subroutine RGC_init(num_homogMech)
allocate(homogState(ho)%state0 (sizeState,Nmembers), source=0.0_pReal)
allocate(homogState(ho)%state (sizeState,Nmembers), source=0.0_pReal)
stt%relaxationVector => homogState(ho)%state(1:nIntFaceTot,:)
st0%relaxationVector => homogState(ho)%state0(1:nIntFaceTot,:)
stt%relaxationVector => homogState(ho)%state(1:nIntFaceTot,:)
st0%relaxationVector => homogState(ho)%state0(1:nIntFaceTot,:)
allocate(dst%volumeDiscrepancy( Nmembers), source=0.0_pReal)
allocate(dst%relaxationRate_avg( Nmembers), source=0.0_pReal)

View File

@ -52,10 +52,11 @@ module subroutine thermal_init()
allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length
allocate(current(ho)%T(count(material_homogenizationID==ho)), source=300.0_pReal)
allocate(current(ho)%T(count(material_homogenizationID==ho)), source=T_ROOM)
allocate(current(ho)%dot_T(count(material_homogenizationID==ho)), source=0.0_pReal)
configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho))
if (configHomogenization%contains('thermal')) then
configHomogenizationThermal => configHomogenization%get('thermal')
#if defined (__GFORTRAN__)
@ -63,13 +64,22 @@ module subroutine thermal_init()
#else
prm%output = configHomogenizationThermal%get_as1dString('output',defaultVal=emptyStringArray)
#endif
select case (configHomogenizationThermal%get_asString('type'))
case ('pass')
call pass_init()
case ('isothermal')
call isotemperature_init()
end select
else
prm%output = emptyStringArray
end if
end associate
end do
call pass_init()
end subroutine thermal_init

View File

@ -1,13 +1,14 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Dummy homogenization scheme for 1 constituent per material point
!> @brief Isotemperature homogenization
!--------------------------------------------------------------------------------------------------
submodule(homogenization:thermal) isotemperature
contains
module subroutine isotemperature_init
module subroutine isotemperature_init()
print'(/,1x,a)', '<<<+- homogenization:thermal:isotemperature init -+>>>'
end subroutine isotemperature_init

View File

@ -7,9 +7,12 @@ submodule(homogenization:thermal) thermal_pass
contains
module subroutine pass_init()
print'(/,1x,a)', '<<<+- homogenization:thermal:pass init -+>>>'
if (homogenization_Nconstituents(1) /= 1) &
call IO_error(211,ext_msg='N_constituents (pass)')
end subroutine pass_init
end submodule thermal_pass

View File

@ -17,27 +17,28 @@ module lattice
private
!--------------------------------------------------------------------------------------------------
! face centered cubic (cF)
integer, dimension(*), parameter :: &
FCC_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for fcc
! cF: face centered cubic (fcc)
integer, dimension(*), parameter :: &
FCC_NTWINSYSTEM = [12] !< # of twin systems per family for fcc
CF_NSLIPSYSTEM = [12, 6] !< # of slip systems per family for cF
integer, dimension(*), parameter :: &
FCC_NTRANSSYSTEM = [12] !< # of transformation systems per family for fcc
CF_NTWINSYSTEM = [12] !< # of twin systems per family for cF
integer, dimension(*), parameter :: &
FCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for fcc
CF_NTRANSSYSTEM = [12] !< # of transformation systems per family for cF
integer, dimension(*), parameter :: &
CF_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for cF
integer, parameter :: &
FCC_NSLIP = sum(FCC_NSLIPSYSTEM), & !< total # of slip systems for fcc
FCC_NTWIN = sum(FCC_NTWINSYSTEM), & !< total # of twin systems for fcc
FCC_NTRANS = sum(FCC_NTRANSSYSTEM), & !< total # of transformation systems for fcc
FCC_NCLEAVAGE = sum(FCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for fcc
CF_NSLIP = sum(CF_NSLIPSYSTEM), & !< total # of slip systems for cF
CF_NTWIN = sum(CF_NTWINSYSTEM), & !< total # of twin systems for cF
CF_NTRANS = sum(CF_NTRANSSYSTEM), & !< total # of transformation systems for cF
CF_NCLEAVAGE = sum(CF_NCLEAVAGESYSTEM) !< total # of cleavage systems for cF
real(pReal), dimension(3+3,FCC_NSLIP), parameter :: &
FCC_SYSTEMSLIP = reshape(real([&
real(pReal), dimension(3+3,CF_NSLIP), parameter :: &
CF_SYSTEMSLIP = reshape(real([&
! <110>{111} systems
0, 1,-1, 1, 1, 1, & ! B2
-1, 0, 1, 1, 1, 1, & ! B4
@ -58,10 +59,10 @@ module lattice
1, 0,-1, 1, 0, 1, &
0, 1, 1, 0, 1,-1, &
0, 1,-1, 0, 1, 1 &
],pReal),shape(FCC_SYSTEMSLIP)) !< fcc slip systems
],pReal),shape(CF_SYSTEMSLIP)) !< cF slip systems
real(pReal), dimension(3+3,FCC_NTWIN), parameter :: &
FCC_SYSTEMTWIN = reshape(real( [&
real(pReal), dimension(3+3,CF_NTWIN), parameter :: &
CF_SYSTEMTWIN = reshape(real( [&
! <112>{111} systems
-2, 1, 1, 1, 1, 1, &
1,-2, 1, 1, 1, 1, &
@ -75,10 +76,10 @@ module lattice
2, 1,-1, -1, 1,-1, &
-1,-2,-1, -1, 1,-1, &
-1, 1, 2, -1, 1,-1 &
],pReal),shape(FCC_SYSTEMTWIN)) !< fcc twin systems
],pReal),shape(CF_SYSTEMTWIN)) !< cF twin systems
integer, dimension(2,FCC_NTWIN), parameter, public :: &
lattice_FCC_TWINNUCLEATIONSLIPPAIR = reshape( [&
integer, dimension(2,CF_NTWIN), parameter, public :: &
lattice_CF_TWINNUCLEATIONSLIPPAIR = reshape( [&
2,3, &
1,3, &
1,2, &
@ -91,34 +92,35 @@ module lattice
11,12, &
10,12, &
10,11 &
],shape(lattice_FCC_TWINNUCLEATIONSLIPPAIR))
],shape(lattice_CF_TWINNUCLEATIONSLIPPAIR))
real(pReal), dimension(3+3,FCC_NCLEAVAGE), parameter :: &
FCC_SYSTEMCLEAVAGE = reshape(real([&
real(pReal), dimension(3+3,CF_NCLEAVAGE), parameter :: &
CF_SYSTEMCLEAVAGE = reshape(real([&
! <001>{001} systems
0, 1, 0, 1, 0, 0, &
0, 0, 1, 0, 1, 0, &
1, 0, 0, 0, 0, 1 &
],pReal),shape(FCC_SYSTEMCLEAVAGE)) !< fcc cleavage systems
],pReal),shape(CF_SYSTEMCLEAVAGE)) !< cF cleavage systems
!--------------------------------------------------------------------------------------------------
! body centered cubic (cI)
integer, dimension(*), parameter :: &
BCC_NSLIPSYSTEM = [12, 12, 24] !< # of slip systems per family for bcc
! cI: body centered cubic (bcc)
integer, dimension(*), parameter :: &
BCC_NTWINSYSTEM = [12] !< # of twin systems per family for bcc
CI_NSLIPSYSTEM = [12, 12, 24] !< # of slip systems per family for cI
integer, dimension(*), parameter :: &
BCC_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for bcc
CI_NTWINSYSTEM = [12] !< # of twin systems per family for cI
integer, dimension(*), parameter :: &
CI_NCLEAVAGESYSTEM = [3] !< # of cleavage systems per family for cI
integer, parameter :: &
BCC_NSLIP = sum(BCC_NSLIPSYSTEM), & !< total # of slip systems for bcc
BCC_NTWIN = sum(BCC_NTWINSYSTEM), & !< total # of twin systems for bcc
BCC_NCLEAVAGE = sum(BCC_NCLEAVAGESYSTEM) !< total # of cleavage systems for bcc
CI_NSLIP = sum(CI_NSLIPSYSTEM), & !< total # of slip systems for cI
CI_NTWIN = sum(CI_NTWINSYSTEM), & !< total # of twin systems for cI
CI_NCLEAVAGE = sum(CI_NCLEAVAGESYSTEM) !< total # of cleavage systems for cI
real(pReal), dimension(3+3,BCC_NSLIP), parameter :: &
BCC_SYSTEMSLIP = reshape(real([&
real(pReal), dimension(3+3,CI_NSLIP), parameter :: &
CI_SYSTEMSLIP = reshape(real([&
! <111>{110} systems
1,-1, 1, 0, 1, 1, & ! D1
-1,-1, 1, 0, 1, 1, & ! C1
@ -170,10 +172,10 @@ module lattice
1, 1, 1, -3, 2, 1, &
1, 1,-1, 3,-2, 1, &
1,-1, 1, 3, 2,-1 &
],pReal),shape(BCC_SYSTEMSLIP)) !< bcc slip systems
],pReal),shape(CI_SYSTEMSLIP)) !< cI slip systems
real(pReal), dimension(3+3,BCC_NTWIN), parameter :: &
BCC_SYSTEMTWIN = reshape(real([&
real(pReal), dimension(3+3,CI_NTWIN), parameter :: &
CI_SYSTEMTWIN = reshape(real([&
! <111>{112} systems
-1, 1, 1, 2, 1, 1, &
1, 1, 1, -2, 1, 1, &
@ -187,30 +189,31 @@ module lattice
1,-1, 1, -1, 1, 2, &
-1, 1, 1, 1,-1, 2, &
1, 1, 1, 1, 1,-2 &
],pReal),shape(BCC_SYSTEMTWIN)) !< bcc twin systems
],pReal),shape(CI_SYSTEMTWIN)) !< cI twin systems
real(pReal), dimension(3+3,BCC_NCLEAVAGE), parameter :: &
BCC_SYSTEMCLEAVAGE = reshape(real([&
real(pReal), dimension(3+3,CI_NCLEAVAGE), parameter :: &
CI_SYSTEMCLEAVAGE = reshape(real([&
! <001>{001} systems
0, 1, 0, 1, 0, 0, &
0, 0, 1, 0, 1, 0, &
1, 0, 0, 0, 0, 1 &
],pReal),shape(BCC_SYSTEMCLEAVAGE)) !< bcc cleavage systems
],pReal),shape(CI_SYSTEMCLEAVAGE)) !< cI cleavage systems
!--------------------------------------------------------------------------------------------------
! hexagonal (hP)
integer, dimension(*), parameter :: &
HEX_NSLIPSYSTEM = [3, 3, 6, 12, 6] !< # of slip systems per family for hex
! hP: hexagonal [close packed] (hex, hcp)
integer, dimension(*), parameter :: &
HEX_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hex
HP_NSLIPSYSTEM = [3, 3, 6, 12, 6] !< # of slip systems per family for hP
integer, dimension(*), parameter :: &
HP_NTWINSYSTEM = [6, 6, 6, 6] !< # of slip systems per family for hP
integer, parameter :: &
HEX_NSLIP = sum(HEX_NSLIPSYSTEM), & !< total # of slip systems for hex
HEX_NTWIN = sum(HEX_NTWINSYSTEM) !< total # of twin systems for hex
HP_NSLIP = sum(HP_NSLIPSYSTEM), & !< total # of slip systems for hP
HP_NTWIN = sum(HP_NTWINSYSTEM) !< total # of twin systems for hP
real(pReal), dimension(4+4,HEX_NSLIP), parameter :: &
HEX_SYSTEMSLIP = reshape(real([&
real(pReal), dimension(4+4,HP_NSLIP), parameter :: &
HP_SYSTEMSLIP = reshape(real([&
! <-1-1.0>{00.1}/basal systems (independent of c/a-ratio)
2, -1, -1, 0, 0, 0, 0, 1, &
-1, 2, -1, 0, 0, 0, 0, 1, &
@ -246,10 +249,10 @@ module lattice
1, 1, -2, 3, -1, -1, 2, 2, &
-1, 2, -1, 3, 1, -2, 1, 2, &
-2, 1, 1, 3, 2, -1, -1, 2 &
],pReal),shape(HEX_SYSTEMSLIP)) !< hex slip systems, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
],pReal),shape(HP_SYSTEMSLIP)) !< hP slip systems, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
real(pReal), dimension(4+4,HEX_NTWIN), parameter :: &
HEX_SYSTEMTWIN = reshape(real([&
real(pReal), dimension(4+4,HP_NTWIN), parameter :: &
HP_SYSTEMTWIN = reshape(real([&
! <-10.1>{10.2} systems, shear = (3-(c/a)^2)/(sqrt(3) c/a)
! tension in Co, Mg, Zr, Ti, and Be; compression in Cd and Zn
-1, 0, 1, 1, 1, 0, -1, 2, & !
@ -282,18 +285,19 @@ module lattice
-1, -1, 2, -3, -1, -1, 2, 2, &
1, -2, 1, -3, 1, -2, 1, 2, &
2, -1, -1, -3, 2, -1, -1, 2 &
],pReal),shape(HEX_SYSTEMTWIN)) !< hex twin systems, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
],pReal),shape(HP_SYSTEMTWIN)) !< hP twin systems, sorted by P. Eisenlohr CCW around <c> starting next to a_1 axis
!--------------------------------------------------------------------------------------------------
! body centered tetragonal (tI)
! tI: body centered tetragonal (bct)
integer, dimension(*), parameter :: &
BCT_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for bct
TI_NSLIPSYSTEM = [2, 2, 2, 4, 2, 4, 2, 2, 4, 8, 4, 8, 8 ] !< # of slip systems per family for tI
integer, parameter :: &
BCT_NSLIP = sum(BCT_NSLIPSYSTEM) !< total # of slip systems for bct
TI_NSLIP = sum(TI_NSLIPSYSTEM) !< total # of slip systems for tI
real(pReal), dimension(3+3,BCT_NSLIP), parameter :: &
BCT_SYSTEMSLIP = reshape(real([&
real(pReal), dimension(3+3,TI_NSLIP), parameter :: &
TI_SYSTEMSLIP = reshape(real([&
! {100)<001] systems
0, 0, 1, 1, 0, 0, &
0, 0, 1, 0, 1, 0, &
@ -359,7 +363,7 @@ module lattice
1,-1, 1, -2,-1, 1, &
-1, 1, 1, -1,-2, 1, &
1, 1, 1, 1,-2, 1 &
],pReal),shape(BCT_SYSTEMSLIP)) !< bct slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x)
],pReal),shape(TI_SYSTEMSLIP)) !< tI slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x)
interface lattice_forestProjection_edge
@ -428,8 +432,8 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
f, & !< index of my family
s !< index of my system in current family
integer, dimension(HEX_NTWIN), parameter :: &
HEX_SHEARTWIN = reshape( [&
integer, dimension(HP_NTWIN), parameter :: &
HP_SHEARTWIN = reshape( [&
1, & ! <-10.1>{10.2}
1, &
1, &
@ -454,7 +458,7 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
4, &
4, &
4 &
],[HEX_NTWIN]) ! indicator to formulas below
],[HP_NTWIN]) !< indicator to formulas below
a = 0
myFamilies: do f = 1,size(Ntwin,1)
@ -466,8 +470,8 @@ function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(character
case('hP')
if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) &
call IO_error(131,ext_msg='lattice_characteristicShear_Twin')
p = sum(HEX_NTWINSYSTEM(1:f-1))+s
select case(HEX_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29
p = sum(HP_NTWINSYSTEM(1:f-1))+s
select case(HP_SHEARTWIN(p)) ! from Christian & Mahajan 1995 p.29
case (1) ! <-10.1>{10.2}
characteristicShear(a) = (3.0_pReal-cOverA**2)/sqrt(3.0_pReal)/CoverA
case (2) ! <11.6>{-1-1.1}
@ -504,13 +508,13 @@ function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
select case(lattice)
case('cF')
coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,&
coordinateSystem = buildCoordinateSystem(Ntwin,CF_NSLIPSYSTEM,CF_SYSTEMTWIN,&
lattice,0.0_pReal)
case('cI')
coordinateSystem = buildCoordinateSystem(Ntwin,BCC_NSLIPSYSTEM,BCC_SYSTEMTWIN,&
coordinateSystem = buildCoordinateSystem(Ntwin,CI_NSLIPSYSTEM,CI_SYSTEMTWIN,&
lattice,0.0_pReal)
case('hP')
coordinateSystem = buildCoordinateSystem(Ntwin,HEX_NSLIPSYSTEM,HEX_SYSTEMTWIN,&
coordinateSystem = buildCoordinateSystem(Ntwin,HP_NSLIPSYSTEM,HP_SYSTEMTWIN,&
lattice,cOverA)
case default
call IO_error(137,ext_msg='lattice_C66_twin: '//trim(lattice))
@ -528,12 +532,12 @@ end function lattice_C66_twin
!> @brief Rotated elasticity matrices for transformation in 6x6-matrix notation
!--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
cOverA_trans,a_fcc,a_bcc)
cOverA_trans,a_cF,a_cI)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol)
real(pReal), dimension(6,6), intent(in) :: C_parent66
real(pReal), optional, intent(in) :: cOverA_trans, a_fcc, a_bcc
real(pReal), optional, intent(in) :: cOverA_trans, a_cF, a_cI
real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans
real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66
@ -562,8 +566,8 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
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 = lattice_symmetrize_C66(C_target_unrotated66,'hP')
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) &
elseif (lattice_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then
if (a_cI <= 0.0_pReal .or. a_cF <= 0.0_pReal) &
call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target))
C_target_unrotated66 = C_parent66
else
@ -575,7 +579,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
call IO_error(135,el=i,ext_msg='matrix diagonal "el"ement in transformation')
enddo
call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_fcc,a_bcc)
call buildTransformationSystem(Q,S,Ntrans,cOverA_trans,a_cF,a_cI)
do i = 1,sum(Ntrans)
call R%fromMatrix(Q(1:3,1:3,i))
@ -586,7 +590,7 @@ function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
!--------------------------------------------------------------------------------------------------
!> @brief Non-schmid projections for bcc with up to 6 coefficients
!> @brief Non-schmid projections for cI with up to 6 coefficients
! https://doi.org/10.1016/j.actamat.2012.03.053, eq. (17)
! https://doi.org/10.1016/j.actamat.2008.07.037, table 1
!--------------------------------------------------------------------------------------------------
@ -605,7 +609,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
if (abs(sense) /= 1) error stop 'Sense in lattice_nonSchmidMatrix'
coordinateSystem = buildCoordinateSystem(Nslip,BCC_NSLIPSYSTEM,BCC_SYSTEMSLIP,'cI',0.0_pReal)
coordinateSystem = buildCoordinateSystem(Nslip,CI_NSLIPSYSTEM,CI_SYSTEMSLIP,'cI',0.0_pReal)
coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip))*real(sense,pReal) ! convert unidirectional coordinate system
nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'cI',0.0_pReal) ! Schmid contribution
@ -636,8 +640,8 @@ end function lattice_nonSchmidMatrix
!--------------------------------------------------------------------------------------------------
!> @brief Slip-slip interaction matrix
!> @details only active slip systems are considered
!> @details https://doi.org/10.1016/j.actamat.2016.12.040 (fcc: Tab S4-1, bcc: Tab S5-1)
!> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hex: Tab 3b)
!> @details https://doi.org/10.1016/j.actamat.2016.12.040 (cF: Tab S4-1, cI: Tab S5-1)
!> @details https://doi.org/10.1016/j.ijplas.2014.06.010 (hP: Tab 3b)
!--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
@ -650,8 +654,8 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
integer, dimension(:,:), allocatable :: interactionTypes
integer, dimension(FCC_NSLIP,FCC_NSLIP), parameter :: &
FCC_INTERACTIONSLIPSLIP = reshape( [&
integer, dimension(CF_NSLIP,CF_NSLIP), parameter :: &
CF_INTERACTIONSLIPSLIP = reshape( [&
1, 2, 2, 4, 7, 5, 3, 5, 5, 4, 6, 7, 10,11,10,11,12,13, & ! -----> acting (forest)
2, 1, 2, 7, 4, 5, 6, 4, 7, 5, 3, 5, 10,11,12,13,10,11, & ! |
2, 2, 1, 5, 5, 3, 6, 7, 4, 7, 6, 4, 12,13,10,11,10,11, & ! |
@ -671,7 +675,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
11,13,11,10,12,10,10,12,10,11,13,11, 9, 9, 8, 1, 9, 9, &
12,10,10,13,11,11,12,10,10,13,11,11, 9, 9, 9, 9, 1, 8, &
13,11,11,12,10,10,13,11,11,12,10,10, 9, 9, 9, 9, 8, 1 &
],shape(FCC_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for fcc / Madec 2017 (https://doi.org/10.1016/j.actamat.2016.12.040)
],shape(CF_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for cF / Madec 2017 (https://doi.org/10.1016/j.actamat.2016.12.040)
!< 1: self interaction --> alpha 0
!< 2: coplanar interaction --> alpha copla
!< 3: collinear interaction --> alpha coli
@ -686,8 +690,8 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
!<12: crossing btw one {110} and one {111} plane
!<13: collinear btw one {110} and one {111} plane
integer, dimension(BCC_NSLIP,BCC_NSLIP), parameter :: &
BCC_INTERACTIONSLIPSLIP = reshape( [&
integer, dimension(CI_NSLIP,CI_NSLIP), parameter :: &
CI_INTERACTIONSLIPSLIP = reshape( [&
1, 3, 6, 6, 7, 5, 4, 2, 4, 2, 7, 5, 18,18,11, 8, 9,13,17,14,13, 9,17,14, 28,25,28,28,25,28,28,28,28,25,28,28,25,28,28,28,28,28,28,25,28,28,28,25, &! -----> acting (forest)
3, 1, 6, 6, 4, 2, 7, 5, 7, 5, 4, 2, 18,18, 8,11,13, 9,14,17, 9,13,14,17, 25,28,28,28,28,25,28,28,25,28,28,28,28,25,28,28,28,28,25,28,28,28,25,28, &! |
6, 6, 1, 3, 5, 7, 2, 4, 5, 7, 2, 4, 11, 8,18,18,17,14, 9,13,17,14,13, 9, 28,28,28,25,28,28,25,28,28,28,28,25,28,28,25,28,28,25,28,28,28,25,28,28, &! |
@ -738,7 +742,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
28,28,25,28,28,28,25,28,28,28,25,28, 28,26,28,28,28,28,26,28,28,28,28,26, 28,28,28,27,28,28,27,28,28,28,28,27,28,28,27,28,28,27,28,28,28, 1,28,28, &
28,25,28,28,28,25,28,28,28,28,28,25, 28,28,26,28,28,26,28,28,26,28,28,28, 27,28,28,28,28,27,28,28,27,28,28,28,28,27,28,28,28,28,27,28,28,28, 1,28, &
25,28,28,28,28,28,28,25,28,25,28,28, 28,28,28,26,26,28,28,28,28,26,28,28, 28,27,28,28,27,28,28,28,28,27,28,28,27,28,28,28,28,28,28,27,28,28,28, 1 &
],shape(BCC_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for bcc / Madec 2017 (https://doi.org/10.1016/j.actamat.2016.12.040)
],shape(CI_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for cI / Madec 2017 (https://doi.org/10.1016/j.actamat.2016.12.040)
!< 1: self interaction --> alpha 0
!< 2: collinear interaction --> alpha 1
!< 3: coplanar interaction --> alpha 2
@ -751,8 +755,8 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
!< 27: {123}-{123} collinear
!< 28: other interaction
integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: &
HEX_INTERACTIONSLIPSLIP = reshape( [&
integer, dimension(HP_NSLIP,HP_NSLIP), parameter :: &
HP_INTERACTIONSLIPSLIP = reshape( [&
! basal prism 1. pyr<a> 1. pyr<c+a> 2. pyr<c+a>
1, 2, 2, 3, 4, 4, 9,10, 9, 9,10, 9, 20,21,22,22,21,20,20,21,22,22,21,20, 47,47,48,47,47,48, & ! -----> acting (forest)
2, 1, 2, 4, 3, 4, 10, 9, 9,10, 9, 9, 22,22,21,20,20,21,22,22,21,20,20,21, 47,48,47,47,48,47, & ! | basal
@ -788,7 +792,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
68,68,69, 66,66,67, 64,64,65,64,64,65, 60,63,63,60,62,62,60,61,61,60,62,62, 59,58,58,57,58,58, & ! 2. pyr<c+a>
68,69,68, 66,67,66, 65,64,64,65,64,64, 62,62,60,63,63,60,62,62,60,61,61,60, 58,59,58,58,57,58, &
69,68,68, 67,66,66, 64,65,64,64,65,64, 61,60,62,62,60,63,63,60,62,62,60,61, 58,58,59,58,58,57 &
],shape(HEX_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hex (onion peel naming scheme)
],shape(HP_INTERACTIONSLIPSLIP)) !< Slip-slip interaction types for hP (onion peel naming scheme)
!< 10.1016/j.ijplas.2014.06.010 table 3
!< 10.1080/14786435.2012.699689 table 2 and 3
!< index & label & description
@ -862,8 +866,8 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
!< 68 & 12 & 2. order pyramidal <c+a>/basal non-collinear
!< 69 & 11 & 2. order pyramidal <c+a>/basal semi-collinear
integer, dimension(BCT_NSLIP,BCT_NSLIP), parameter :: &
BCT_INTERACTIONSLIPSLIP = reshape( [&
integer, dimension(TI_NSLIP,TI_NSLIP), parameter :: &
TI_INTERACTIONSLIPSLIP = reshape( [&
1, 2, 3, 3, 7, 7, 13, 13, 13, 13, 21, 21, 31, 31, 31, 31, 43, 43, 57, 57, 73, 73, 73, 73, 91, 91, 91, 91, 91, 91, 91, 91, 111, 111, 111, 111, 133,133,133,133,133,133,133,133, 157,157,157,157,157,157,157,157, & ! -----> acting
2, 1, 3, 3, 7, 7, 13, 13, 13, 13, 21, 21, 31, 31, 31, 31, 43, 43, 57, 57, 73, 73, 73, 73, 91, 91, 91, 91, 91, 91, 91, 91, 111, 111, 111, 111, 133,133,133,133,133,133,133,133, 157,157,157,157,157,157,157,157, & ! |
! |
@ -928,22 +932,22 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(
182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 169,170,170,170,170,169,170,170, &
182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 169,170,170,170,170,170,169,170, &
182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 169,170,170,170,170,170,170,169 &
],shape(BCT_INTERACTIONSLIPSLIP))
],shape(TI_INTERACTIONSLIPSLIP))
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONSLIPSLIP
NslipMax = FCC_NSLIPSYSTEM
interactionTypes = CF_INTERACTIONSLIPSLIP
NslipMax = CF_NSLIPSYSTEM
case('cI')
interactionTypes = BCC_INTERACTIONSLIPSLIP
NslipMax = BCC_NSLIPSYSTEM
interactionTypes = CI_INTERACTIONSLIPSLIP
NslipMax = CI_NSLIPSYSTEM
case('hP')
interactionTypes = HEX_INTERACTIONSLIPSLIP
NslipMax = HEX_NSLIPSYSTEM
interactionTypes = HP_INTERACTIONSLIPSLIP
NslipMax = HP_NSLIPSYSTEM
case('tI')
interactionTypes = BCT_INTERACTIONSLIPSLIP
NslipMax = BCT_NSLIPSYSTEM
interactionTypes = TI_INTERACTIONSLIPSLIP
NslipMax = TI_NSLIPSYSTEM
case default
call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(lattice))
end select
@ -967,8 +971,8 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
integer, dimension(:), allocatable :: NtwinMax
integer, dimension(:,:), allocatable :: interactionTypes
integer, dimension(FCC_NTWIN,FCC_NTWIN), parameter :: &
FCC_INTERACTIONTWINTWIN = reshape( [&
integer, dimension(CF_NTWIN,CF_NTWIN), parameter :: &
CF_INTERACTIONTWINTWIN = reshape( [&
1,1,1,2,2,2,2,2,2,2,2,2, & ! -----> acting
1,1,1,2,2,2,2,2,2,2,2,2, & ! |
1,1,1,2,2,2,2,2,2,2,2,2, & ! |
@ -981,10 +985,10 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1 &
],shape(FCC_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for fcc
],shape(CF_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for cF
integer, dimension(BCC_NTWIN,BCC_NTWIN), parameter :: &
BCC_INTERACTIONTWINTWIN = reshape( [&
integer, dimension(CI_NTWIN,CI_NTWIN), parameter :: &
CI_INTERACTIONTWINTWIN = reshape( [&
1,3,3,3,3,3,3,2,3,3,2,3, & ! -----> acting
3,1,3,3,3,3,2,3,3,3,3,2, & ! |
3,3,1,3,3,2,3,3,2,3,3,3, & ! |
@ -997,12 +1001,12 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
3,3,3,2,2,3,3,3,3,1,3,3, &
2,3,3,3,3,3,3,2,3,3,1,3, &
3,2,3,3,3,3,2,3,3,3,3,1 &
],shape(BCC_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for bcc
],shape(CI_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for cI
!< 1: self interaction
!< 2: collinear interaction
!< 3: other interaction
integer, dimension(HEX_NTWIN,HEX_NTWIN), parameter :: &
HEX_INTERACTIONTWINTWIN = reshape( [&
integer, dimension(HP_NTWIN,HP_NTWIN), parameter :: &
HP_INTERACTIONTWINTWIN = reshape( [&
! <-10.1>{10.2} <11.6>{-1-1.1} <10.-2>{10.1} <11.-3>{11.2}
1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! -----> acting
2, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! |
@ -1031,18 +1035,18 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17, & ! <11.-3>{11.2}
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,16,17, &
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 &
],shape(HEX_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hex
],shape(HP_INTERACTIONTWINTWIN)) !< Twin-twin interaction types for hP
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONTWINTWIN
NtwinMax = FCC_NTWINSYSTEM
interactionTypes = CF_INTERACTIONTWINTWIN
NtwinMax = CF_NTWINSYSTEM
case('cI')
interactionTypes = BCC_INTERACTIONTWINTWIN
NtwinMax = BCC_NTWINSYSTEM
interactionTypes = CI_INTERACTIONTWINTWIN
NtwinMax = CI_NTWINSYSTEM
case('hP')
interactionTypes = HEX_INTERACTIONTWINTWIN
NtwinMax = HEX_NTWINSYSTEM
interactionTypes = HP_INTERACTIONTWINTWIN
NtwinMax = HP_NTWINSYSTEM
case default
call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(lattice))
end select
@ -1066,8 +1070,8 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu
integer, dimension(:), allocatable :: NtransMax
integer, dimension(:,:), allocatable :: interactionTypes
integer, dimension(FCC_NTRANS,FCC_NTRANS), parameter :: &
FCC_INTERACTIONTRANSTRANS = reshape( [&
integer, dimension(CF_NTRANS,CF_NTRANS), parameter :: &
CF_INTERACTIONTRANSTRANS = reshape( [&
1,1,1,2,2,2,2,2,2,2,2,2, & ! -----> acting
1,1,1,2,2,2,2,2,2,2,2,2, & ! |
1,1,1,2,2,2,2,2,2,2,2,2, & ! |
@ -1080,11 +1084,11 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) resu
2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1 &
],shape(FCC_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for fcc
],shape(CF_INTERACTIONTRANSTRANS)) !< Trans-trans interaction types for cF
if (lattice == 'cF') then
interactionTypes = FCC_INTERACTIONTRANSTRANS
NtransMax = FCC_NTRANSSYSTEM
interactionTypes = CF_INTERACTIONTRANSTRANS
NtransMax = CF_NTRANSSYSTEM
else
call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(lattice))
end if
@ -1110,8 +1114,8 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
NtwinMax
integer, dimension(:,:), allocatable :: interactionTypes
integer, dimension(FCC_NTWIN,FCC_NSLIP), parameter :: &
FCC_INTERACTIONSLIPTWIN = reshape( [&
integer, dimension(CF_NTWIN,CF_NSLIP), parameter :: &
CF_INTERACTIONSLIPTWIN = reshape( [&
1,1,1,3,3,3,2,2,2,3,3,3, & ! -----> twin (acting)
1,1,1,3,3,3,3,3,3,2,2,2, & ! |
1,1,1,2,2,2,3,3,3,3,3,3, & ! |
@ -1131,12 +1135,12 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4 &
],shape(FCC_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for fcc
],shape(CF_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for cF
!< 1: coplanar interaction
!< 2: screw trace between slip system and twin habit plane (easy cross slip)
!< 3: other interaction
integer, dimension(BCC_NTWIN,BCC_NSLIP), parameter :: &
BCC_INTERACTIONSLIPTWIN = reshape( [&
integer, dimension(CI_NTWIN,CI_NSLIP), parameter :: &
CI_INTERACTIONSLIPTWIN = reshape( [&
3,3,3,2,2,3,3,3,3,2,3,3, & ! -----> twin (acting)
3,3,2,3,3,2,3,3,2,3,3,3, & ! |
3,2,3,3,3,3,2,3,3,3,3,2, & ! |
@ -1187,14 +1191,14 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4 &
],shape(BCC_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for bcc
],shape(CI_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for cI
!< 1: coplanar interaction
!< 2: screw trace between slip system and twin habit plane (easy cross slip)
!< 3: other interaction
!< 4: other interaction with slip family {123}
integer, dimension(HEX_NTWIN,HEX_NSLIP), parameter :: &
HEX_INTERACTIONSLIPTWIN = reshape( [&
integer, dimension(HP_NTWIN,HP_NSLIP), parameter :: &
HP_INTERACTIONSLIPTWIN = reshape( [&
! <-10.1>{10.2} <11.6>{-1-1.1} <10.-2>{10.1} <11.-3>{11.2}
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! ----> twin (acting)
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! | basal
@ -1230,21 +1234,21 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) r
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, & ! 2. pyr<c+a>
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20, &
17,17,17,17,17,17, 18,18,18,18,18,18, 19,19,19,19,19,19, 20,20,20,20,20,20 &
],shape(HEX_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hex
],shape(HP_INTERACTIONSLIPTWIN)) !< Slip-twin interaction types for hP
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONSLIPTWIN
NslipMax = FCC_NSLIPSYSTEM
NtwinMax = FCC_NTWINSYSTEM
interactionTypes = CF_INTERACTIONSLIPTWIN
NslipMax = CF_NSLIPSYSTEM
NtwinMax = CF_NTWINSYSTEM
case('cI')
interactionTypes = BCC_INTERACTIONSLIPTWIN
NslipMax = BCC_NSLIPSYSTEM
NtwinMax = BCC_NTWINSYSTEM
interactionTypes = CI_INTERACTIONSLIPTWIN
NslipMax = CI_NSLIPSYSTEM
NtwinMax = CI_NTWINSYSTEM
case('hP')
interactionTypes = HEX_INTERACTIONSLIPTWIN
NslipMax = HEX_NSLIPSYSTEM
NtwinMax = HEX_NTWINSYSTEM
interactionTypes = HP_INTERACTIONSLIPTWIN
NslipMax = HP_NSLIPSYSTEM
NtwinMax = HP_NTWINSYSTEM
case default
call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(lattice))
end select
@ -1270,8 +1274,8 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice)
NtransMax
integer, dimension(:,:), allocatable :: interactionTypes
integer, dimension(FCC_NTRANS,FCC_NSLIP), parameter :: &
FCC_INTERACTIONSLIPTRANS = reshape( [&
integer, dimension(CF_NTRANS,CF_NSLIP), parameter :: &
CF_INTERACTIONSLIPTRANS = reshape( [&
1,1,1,3,3,3,2,2,2,3,3,3, & ! -----> trans (acting)
1,1,1,3,3,3,3,3,3,2,2,2, & ! |
1,1,1,2,2,2,3,3,3,3,3,3, & ! |
@ -1291,13 +1295,13 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice)
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4 &
],shape(FCC_INTERACTIONSLIPTRANS)) !< Slip-trans interaction types for fcc
],shape(CF_INTERACTIONSLIPTRANS)) !< Slip-trans interaction types for cF
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONSLIPTRANS
NslipMax = FCC_NSLIPSYSTEM
NtransMax = FCC_NTRANSSYSTEM
interactionTypes = CF_INTERACTIONSLIPTRANS
NslipMax = CF_NSLIPSYSTEM
NtransMax = CF_NTRANSSYSTEM
case default
call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(lattice))
end select
@ -1323,14 +1327,14 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r
NslipMax
integer, dimension(:,:), allocatable :: interactionTypes
integer, dimension(FCC_NSLIP,FCC_NTWIN), parameter :: &
FCC_INTERACTIONTWINSLIP = 1 !< Twin-slip interaction types for fcc
integer, dimension(CF_NSLIP,CF_NTWIN), parameter :: &
CF_INTERACTIONTWINSLIP = 1 !< Twin-slip interaction types for cF
integer, dimension(BCC_NSLIP,BCC_NTWIN), parameter :: &
BCC_INTERACTIONTWINSLIP = 1 !< Twin-slip interaction types for bcc
integer, dimension(CI_NSLIP,CI_NTWIN), parameter :: &
CI_INTERACTIONTWINSLIP = 1 !< Twin-slip interaction types for cI
integer, dimension(HEX_NSLIP,HEX_NTWIN), parameter :: &
HEX_INTERACTIONTWINSLIP = reshape( [&
integer, dimension(HP_NSLIP,HP_NTWIN), parameter :: &
HP_INTERACTIONTWINSLIP = reshape( [&
! basal prism 1. pyr<a> 1. pyr<c+a> 2. pyr<c+a>
1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! ----> slip (acting)
1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9, 13,13,13,13,13,13,13,13,13,13,13,13, 17,17,17,17,17,17, & ! |
@ -1359,21 +1363,21 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) r
4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20, & ! <11.-3>{11.2}
4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20, &
4, 4, 4, 8, 8, 8, 12,12,12,12,12,12, 16,16,16,16,16,16,16,16,16,16,16,16, 20,20,20,20,20,20 &
],shape(HEX_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hex
],shape(HP_INTERACTIONTWINSLIP)) !< Twin-slip interaction types for hP
select case(lattice)
case('cF')
interactionTypes = FCC_INTERACTIONTWINSLIP
NtwinMax = FCC_NTWINSYSTEM
NslipMax = FCC_NSLIPSYSTEM
interactionTypes = CF_INTERACTIONTWINSLIP
NtwinMax = CF_NTWINSYSTEM
NslipMax = CF_NSLIPSYSTEM
case('cI')
interactionTypes = BCC_INTERACTIONTWINSLIP
NtwinMax = BCC_NTWINSYSTEM
NslipMax = BCC_NSLIPSYSTEM
interactionTypes = CI_INTERACTIONTWINSLIP
NtwinMax = CI_NTWINSYSTEM
NslipMax = CI_NSLIPSYSTEM
case('hP')
interactionTypes = HEX_INTERACTIONTWINSLIP
NtwinMax = HEX_NTWINSYSTEM
NslipMax = HEX_NSLIPSYSTEM
interactionTypes = HP_INTERACTIONTWINSLIP
NtwinMax = HP_NTWINSYSTEM
NslipMax = HP_NSLIPSYSTEM
case default
call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(lattice))
end select
@ -1401,17 +1405,17 @@ function lattice_SchmidMatrix_slip(Nslip,lattice,cOverA) result(SchmidMatrix)
select case(lattice)
case('cF')
NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP
NslipMax = CF_NSLIPSYSTEM
slipSystems = CF_SYSTEMSLIP
case('cI')
NslipMax = BCC_NSLIPSYSTEM
slipSystems = BCC_SYSTEMSLIP
NslipMax = CI_NSLIPSYSTEM
slipSystems = CI_SYSTEMSLIP
case('hP')
NslipMax = HEX_NSLIPSYSTEM
slipSystems = HEX_SYSTEMSLIP
NslipMax = HP_NSLIPSYSTEM
slipSystems = HP_SYSTEMSLIP
case('tI')
NslipMax = BCT_NSLIPSYSTEM
slipSystems = BCT_SYSTEMSLIP
NslipMax = TI_NSLIPSYSTEM
slipSystems = TI_SYSTEMSLIP
case default
allocate(NslipMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(lattice))
@ -1451,14 +1455,14 @@ function lattice_SchmidMatrix_twin(Ntwin,lattice,cOverA) result(SchmidMatrix)
select case(lattice)
case('cF')
NtwinMax = FCC_NTWINSYSTEM
twinSystems = FCC_SYSTEMTWIN
NtwinMax = CF_NTWINSYSTEM
twinSystems = CF_SYSTEMTWIN
case('cI')
NtwinMax = BCC_NTWINSYSTEM
twinSystems = BCC_SYSTEMTWIN
NtwinMax = CI_NTWINSYSTEM
twinSystems = CI_SYSTEMTWIN
case('hP')
NtwinMax = HEX_NTWINSYSTEM
twinSystems = HEX_SYSTEMTWIN
NtwinMax = HP_NTWINSYSTEM
twinSystems = HP_SYSTEMTWIN
case default
allocate(NtwinMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(lattice))
@ -1484,11 +1488,11 @@ end function lattice_SchmidMatrix_twin
!> @brief Schmid matrix for transformation
!> details only active twin systems are considered
!--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_fcc,a_bcc) result(SchmidMatrix)
function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_cF,a_cI) result(SchmidMatrix)
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=2), intent(in) :: lattice_target !< Bravais lattice (Pearson symbol)
real(pReal), optional, intent(in) :: cOverA, a_bcc, a_fcc
real(pReal), optional, intent(in) :: cOverA, a_cI, a_cF
real(pReal), dimension(3,3,sum(Ntrans)) :: SchmidMatrix
real(pReal), dimension(3,3,sum(Ntrans)) :: devNull
@ -1498,10 +1502,10 @@ function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_fcc,a_bcc) re
if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) &
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA=cOverA)
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) &
else if (lattice_target == 'cI' .and. present(a_cF) .and. present(a_cI)) then
if (a_cI <= 0.0_pReal .or. a_cF <= 0.0_pReal) &
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,a_cF=a_cF,a_cI=a_cI)
else
call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
end if
@ -1527,11 +1531,11 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,lattice,cOverA) result(SchmidMa
select case(lattice)
case('cF')
NcleavageMax = FCC_NCLEAVAGESYSTEM
cleavageSystems = FCC_SYSTEMCLEAVAGE
NcleavageMax = CF_NCLEAVAGESYSTEM
cleavageSystems = CF_SYSTEMCLEAVAGE
case('cI')
NcleavageMax = BCC_NCLEAVAGESYSTEM
cleavageSystems = BCC_SYSTEMCLEAVAGE
NcleavageMax = CI_NCLEAVAGESYSTEM
cleavageSystems = CI_SYSTEMCLEAVAGE
case default
allocate(NcleavageMax(0))
call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(lattice))
@ -1623,17 +1627,17 @@ function lattice_labels_slip(Nslip,lattice) result(labels)
select case(lattice)
case('cF')
NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP
NslipMax = CF_NSLIPSYSTEM
slipSystems = CF_SYSTEMSLIP
case('cI')
NslipMax = BCC_NSLIPSYSTEM
slipSystems = BCC_SYSTEMSLIP
NslipMax = CI_NSLIPSYSTEM
slipSystems = CI_SYSTEMSLIP
case('hP')
NslipMax = HEX_NSLIPSYSTEM
slipSystems = HEX_SYSTEMSLIP
NslipMax = HP_NSLIPSYSTEM
slipSystems = HP_SYSTEMSLIP
case('tI')
NslipMax = BCT_NSLIPSYSTEM
slipSystems = BCT_SYSTEMSLIP
NslipMax = TI_NSLIPSYSTEM
slipSystems = TI_SYSTEMSLIP
case default
call IO_error(137,ext_msg='lattice_labels_slip: '//trim(lattice))
end select
@ -1737,14 +1741,14 @@ function lattice_labels_twin(Ntwin,lattice) result(labels)
select case(lattice)
case('cF')
NtwinMax = FCC_NTWINSYSTEM
twinSystems = FCC_SYSTEMTWIN
NtwinMax = CF_NTWINSYSTEM
twinSystems = CF_SYSTEMTWIN
case('cI')
NtwinMax = BCC_NTWINSYSTEM
twinSystems = BCC_SYSTEMTWIN
NtwinMax = CI_NTWINSYSTEM
twinSystems = CI_SYSTEMTWIN
case('hP')
NtwinMax = HEX_NTWINSYSTEM
twinSystems = HEX_SYSTEMTWIN
NtwinMax = HP_NTWINSYSTEM
twinSystems = HP_SYSTEMTWIN
case default
call IO_error(137,ext_msg='lattice_labels_twin: '//trim(lattice))
end select
@ -1823,17 +1827,17 @@ function coordinateSystem_slip(Nslip,lattice,cOverA) result(coordinateSystem)
select case(lattice)
case('cF')
NslipMax = FCC_NSLIPSYSTEM
slipSystems = FCC_SYSTEMSLIP
NslipMax = CF_NSLIPSYSTEM
slipSystems = CF_SYSTEMSLIP
case('cI')
NslipMax = BCC_NSLIPSYSTEM
slipSystems = BCC_SYSTEMSLIP
NslipMax = CI_NSLIPSYSTEM
slipSystems = CI_SYSTEMSLIP
case('hP')
NslipMax = HEX_NSLIPSYSTEM
slipSystems = HEX_SYSTEMSLIP
NslipMax = HP_NSLIPSYSTEM
slipSystems = HP_SYSTEMSLIP
case('tI')
NslipMax = BCT_NSLIPSYSTEM
slipSystems = BCT_SYSTEMSLIP
NslipMax = TI_NSLIPSYSTEM
slipSystems = TI_SYSTEMSLIP
case default
allocate(NslipMax(0))
call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(lattice))
@ -1962,24 +1966,24 @@ end function buildCoordinateSystem
!--------------------------------------------------------------------------------------------------
!> @brief Helper function to define transformation systems
! Needed to calculate Schmid matrix and rotated stiffness matrices.
! @details: set c/a = 0.0 for fcc -> bcc transformation
! set a_Xcc = 0.0 for fcc -> hex transformation
! @details: use c/a for cF -> cI transformation
! use a_cX for cF -> hP transformation
!--------------------------------------------------------------------------------------------------
subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_cF,a_cI)
integer, dimension(:), intent(in) :: &
integer, dimension(:), intent(in) :: &
Ntrans
real(pReal), dimension(3,3,sum(Ntrans)), intent(out) :: &
Q, & !< Total rotation: Q = R*B
S !< Eigendeformation tensor for phase transformation
real(pReal), optional, intent(in) :: &
cOverA, & !< c/a for target hex lattice
a_bcc, & !< lattice parameter a for bcc target lattice
a_fcc !< lattice parameter a for fcc parent lattice
real(pReal), optional, intent(in) :: &
cOverA, & !< c/a for target hP lattice
a_cF, & !< lattice parameter a for cF target lattice
a_cI !< lattice parameter a for cI parent lattice
type(tRotation) :: &
R, & !< Pitsch rotation
B !< Rotation of fcc to Bain coordinate system
B !< Rotation of cF to Bain coordinate system
real(pReal), dimension(3,3) :: &
U, & !< Bain deformation
ss, sd
@ -1987,8 +1991,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
x, y, z
integer :: &
i
real(pReal), dimension(3+3,FCC_NTRANS), parameter :: &
FCCTOHEX_SYSTEMTRANS = reshape(real( [&
real(pReal), dimension(3+3,CF_NTRANS), parameter :: &
CFTOHP_SYSTEMTRANS = reshape(real( [&
-2, 1, 1, 1, 1, 1, &
1,-2, 1, 1, 1, 1, &
1, 1,-2, 1, 1, 1, &
@ -2001,10 +2005,10 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
2, 1,-1, -1, 1,-1, &
-1,-2,-1, -1, 1,-1, &
-1, 1, 2, -1, 1,-1 &
],pReal),shape(FCCTOHEX_SYSTEMTRANS))
],pReal),shape(CFTOHP_SYSTEMTRANS))
real(pReal), dimension(4,fcc_Ntrans), parameter :: &
FCCTOBCC_SYSTEMTRANS = reshape([&
real(pReal), dimension(4,cF_Ntrans), parameter :: &
CFTOCI_SYSTEMTRANS = reshape([&
0.0, 1.0, 0.0, 10.26, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3)
0.0,-1.0, 0.0, 10.26, &
0.0, 0.0, 1.0, 10.26, &
@ -2017,10 +2021,10 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
-1.0, 0.0, 0.0, 10.26, &
0.0, 1.0, 0.0, 10.26, &
0.0,-1.0, 0.0, 10.26 &
],shape(FCCTOBCC_SYSTEMTRANS))
],shape(CFTOCI_SYSTEMTRANS))
integer, dimension(9,fcc_Ntrans), parameter :: &
FCCTOBCC_BAINVARIANT = reshape( [&
integer, dimension(9,cF_Ntrans), parameter :: &
CFTOCI_BAINVARIANT = reshape( [&
1, 0, 0, 0, 1, 0, 0, 0, 1, & ! Pitsch OR (Ma & Hartmaier 2014, Table 3)
1, 0, 0, 0, 1, 0, 0, 0, 1, &
1, 0, 0, 0, 1, 0, 0, 0, 1, &
@ -2033,11 +2037,11 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
0, 0, 1, 1, 0, 0, 0, 1, 0, &
0, 0, 1, 1, 0, 0, 0, 1, 0, &
0, 0, 1, 1, 0, 0, 0, 1, 0 &
],shape(FCCTOBCC_BAINVARIANT))
],shape(CFTOCI_BAINVARIANT))
real(pReal), dimension(4,fcc_Ntrans), parameter :: &
FCCTOBCC_BAINROT = reshape([&
1.0, 0.0, 0.0, 45.0, & ! Rotate fcc austensite to bain variant
real(pReal), dimension(4,cF_Ntrans), parameter :: &
CFTOCI_BAINROT = reshape([&
1.0, 0.0, 0.0, 45.0, & ! Rotate cF austensite to bain variant
1.0, 0.0, 0.0, 45.0, &
1.0, 0.0, 0.0, 45.0, &
1.0, 0.0, 0.0, 45.0, &
@ -2049,19 +2053,17 @@ 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, &
0.0, 0.0, 1.0, 45.0 &
],shape(FCCTOBCC_BAINROT))
],shape(CFTOCI_BAINROT))
if (present(a_bcc) .and. present(a_fcc)) then
if (present(a_cI) .and. present(a_cF)) then
do i = 1,sum(Ntrans)
call R%fromAxisAngle(FCCTOBCC_SYSTEMTRANS(:,i),degrees=.true.,P=1)
call B%fromAxisAngle(FCCTOBCC_BAINROT(:,i), degrees=.true.,P=1)
x = real(FCCTOBCC_BAINVARIANT(1:3,i),pReal)
y = real(FCCTOBCC_BAINVARIANT(4:6,i),pReal)
z = real(FCCTOBCC_BAINVARIANT(7:9,i),pReal)
call R%fromAxisAngle(CFTOCI_SYSTEMTRANS(:,i),degrees=.true.,P=1)
call B%fromAxisAngle(CFTOCI_BAINROT(:,i), degrees=.true.,P=1)
x = real(CFTOCI_BAINVARIANT(1:3,i),pReal)
y = real(CFTOCI_BAINVARIANT(4:6,i),pReal)
z = real(CFTOCI_BAINVARIANT(7:9,i),pReal)
U = (a_bcc/a_fcc)*math_outer(x,x) &
+ (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) &
+ (a_bcc/a_fcc)*math_outer(z,z) * sqrt(2.0_pReal)
U = (a_cI/a_cF) * (math_outer(x,x) + (math_outer(y,y)+math_outer(z,z)) * sqrt(2.0_pReal))
Q(1:3,1:3,i) = matmul(R%asMatrix(),B%asMatrix())
S(1:3,1:3,i) = matmul(R%asMatrix(),U) - MATH_I3
enddo
@ -2072,8 +2074,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
sd(3,3) = cOverA/sqrt(8.0_pReal/3.0_pReal)
do i = 1,sum(Ntrans)
x = FCCTOHEX_SYSTEMTRANS(1:3,i)/norm2(FCCTOHEX_SYSTEMTRANS(1:3,i))
z = FCCTOHEX_SYSTEMTRANS(4:6,i)/norm2(FCCTOHEX_SYSTEMTRANS(4:6,i))
x = CFTOHP_SYSTEMTRANS(1:3,i)/norm2(CFTOHP_SYSTEMTRANS(1:3,i))
z = CFTOHP_SYSTEMTRANS(4:6,i)/norm2(CFTOHP_SYSTEMTRANS(4:6,i))
y = -math_cross(x,z)
Q(1:3,1,i) = x
Q(1:3,2,i) = y

View File

@ -91,6 +91,11 @@ module phase
integer, intent(in) :: ph
end subroutine damage_results
module subroutine thermal_results(group,ph)
character(len=*), intent(in) :: group
integer, intent(in) :: ph
end subroutine thermal_results
module subroutine mechanical_forward()
end subroutine mechanical_forward
@ -487,6 +492,7 @@ subroutine phase_results()
call mechanical_results(group,ph)
call damage_results(group,ph)
call thermal_results(group,ph)
end do

View File

@ -336,7 +336,7 @@ end subroutine damage_results
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!> @brief Constitutive equation for calculating the rate of change of microstructure.
!--------------------------------------------------------------------------------------------------
function phase_damage_collectDotState(ph,en) result(broken)

View File

@ -180,11 +180,12 @@ submodule(phase) mechanical
end function elastic_nu
end interface
type :: tOutput !< new requested output (per phase)
type :: tOutput !< requested output (per phase)
character(len=pStringLen), allocatable, dimension(:) :: &
label
end type tOutput
type(tOutput), allocatable, dimension(:) :: output_constituent
type(tOutput), allocatable, dimension(:) :: output_mechanical
procedure(integrateStateFPI), pointer :: integrateState
@ -216,7 +217,7 @@ module subroutine mechanical_init(phases)
print'(/,1x,a)', '<<<+- phase:mechanical init -+>>>'
!-------------------------------------------------------------------------------------------------
allocate(output_constituent(phases%length))
allocate(output_mechanical(phases%length))
allocate(phase_mechanical_Fe(phases%length))
allocate(phase_mechanical_Fi(phases%length))
@ -251,12 +252,12 @@ module subroutine mechanical_init(phases)
allocate(phase_mechanical_P(ph)%data(3,3,Nmembers),source=0.0_pReal)
allocate(phase_mechanical_S0(ph)%data(3,3,Nmembers),source=0.0_pReal)
phase => phases%get(ph)
mech => phase%get('mechanical')
phase => phases%get(ph)
mech => phase%get('mechanical')
#if defined(__GFORTRAN__)
output_constituent(ph)%label = output_as1dString(mech)
output_mechanical(ph)%label = output_as1dString(mech)
#else
output_constituent(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
output_mechanical(ph)%label = mech%get_as1dString('output',defaultVal=emptyStringArray)
#endif
enddo
@ -330,7 +331,7 @@ module subroutine mechanical_results(group,ph)
integer, intent(in) :: ph
call crystallite_results(group,ph)
call results(group,ph)
select case(phase_plasticity(ph))
@ -879,9 +880,9 @@ end function integrateStateRK
!--------------------------------------------------------------------------------------------------
!> @brief writes crystallite results to HDF5 output file
!> @brief Write mechanical results to HDF5 output file.
!--------------------------------------------------------------------------------------------------
subroutine crystallite_results(group,ph)
subroutine results(group,ph)
character(len=*), intent(in) :: group
integer, intent(in) :: ph
@ -891,9 +892,9 @@ subroutine crystallite_results(group,ph)
call results_closeGroup(results_addGroup(group//'/mechanical'))
do ou = 1, size(output_constituent(ph)%label)
do ou = 1, size(output_mechanical(ph)%label)
select case (output_constituent(ph)%label(ou))
select case (output_mechanical(ph)%label(ou))
case('F')
call results_writeDataset(phase_mechanical_F(ph)%data,group//'/mechanical/','F',&
'deformation gradient','1')
@ -919,13 +920,13 @@ subroutine crystallite_results(group,ph)
call results_writeDataset(phase_mechanical_S(ph)%data,group//'/mechanical/','S', &
'second Piola-Kirchhoff stress','Pa')
case('O')
call results_writeDataset(to_quaternion(phase_O(ph)%data),group//'/mechanical',output_constituent(ph)%label(ou),&
call results_writeDataset(to_quaternion(phase_O(ph)%data),group//'/mechanical',output_mechanical(ph)%label(ou),&
'crystal orientation as quaternion','q_0 (q_1 q_2 q_3)')
call results_addAttribute('lattice',phase_lattice(ph),group//'/mechanical/'//output_constituent(ph)%label(ou))
call results_addAttribute('lattice',phase_lattice(ph),group//'/mechanical/'//output_mechanical(ph)%label(ou))
if (any(phase_lattice(ph) == ['hP', 'tI'])) &
call results_addAttribute('c/a',phase_cOverA(ph),group//'/mechanical/'//output_constituent(ph)%label(ou))
call results_addAttribute('c/a',phase_cOverA(ph),group//'/mechanical/'//output_mechanical(ph)%label(ou))
end select
enddo
end do
contains
@ -947,7 +948,7 @@ subroutine crystallite_results(group,ph)
end function to_quaternion
end subroutine crystallite_results
end subroutine results
!--------------------------------------------------------------------------------------------------
@ -1335,5 +1336,4 @@ module subroutine phase_set_F(F,co,ce)
end subroutine phase_set_F
end submodule mechanical

View File

@ -200,7 +200,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%n0_sl = lattice_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph))
prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. (N_sl(1) == 12)
if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR
if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_CF_TWINNUCLEATIONSLIPPAIR
rho_mob_0 = pl%get_as1dFloat('rho_mob_0', requiredSize=size(N_sl))
rho_dip_0 = pl%get_as1dFloat('rho_dip_0', requiredSize=size(N_sl))

View File

@ -6,6 +6,7 @@ submodule(phase) thermal
type :: tThermalParameters
real(pReal) :: C_p = 0.0_pReal !< heat capacity
real(pReal), dimension(3,3) :: K = 0.0_pReal !< thermal conductivity
character(len=pStringLen), allocatable, dimension(:) :: output
end type tThermalParameters
integer, dimension(:), allocatable :: &
@ -108,6 +109,11 @@ module subroutine thermal_init(phases)
if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%K(3,3) = thermal%get_asFloat('K_33')
param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph))
#if defined(__GFORTRAN__)
param(ph)%output = output_as1dString(thermal)
#else
param(ph)%output = thermal%get_as1dString('output',defaultVal=emptyStringArray)
#endif
sources => thermal%get('source',defaultVal=emptyList)
thermal_Nsources(ph) = sources%length
else
@ -381,4 +387,35 @@ function thermal_active(source_label,src_length) result(active_source)
end function thermal_active
!----------------------------------------------------------------------------------------------
!< @brief writes damage sources results to HDF5 output file
!----------------------------------------------------------------------------------------------
module subroutine thermal_results(group,ph)
character(len=*), intent(in) :: group
integer, intent(in) :: ph
integer :: ou
if (allocated(param(ph)%output)) then
call results_closeGroup(results_addGroup(group//'thermal'))
else
return
endif
do ou = 1, size(param(ph)%output)
select case(trim(param(ph)%output(ou)))
case ('T')
call results_writeDataset(current(ph)%T,group//'thermal','T', 'temperature','T')
end select
end do
end subroutine thermal_results
end submodule thermal