Merge remote-tracking branch 'origin/development' into internal-restructure

This commit is contained in:
Sharan Roongta 2021-02-19 15:19:42 +01:00
commit db25bc947d
11 changed files with 350 additions and 325 deletions

View File

@ -1 +1 @@
v3.0.0-alpha2-421-ge96352b0e
v3.0.0-alpha2-469-gcf3084a5c

View File

@ -1,3 +1,5 @@
import inspect
import numpy as np
from . import Rotation
@ -7,7 +9,7 @@ from . import tensor
_parameter_doc = \
"""lattice : str
Either a crystal family out of [triclinic, monoclinic, orthorhombic, tetragonal, hexagonal, cubic]
or a Bravais lattice out of [aP, mP, mS, oP, oS, oI, oF, tP, tI, hP, cP, cI, cF].
or a Bravais lattice out of [aP, mP, mS, oP, oS, oI, oF, tP, tI, hP, cP, cI, cF].
When specifying a Bravais lattice, additional lattice parameters might be required:
a : float, optional
Length of lattice parameter "a".
@ -107,8 +109,7 @@ class Orientation(Rotation):
lattice = None,
a = None,b = None,c = None,
alpha = None,beta = None,gamma = None,
degrees = False,
**kwargs):
degrees = False):
"""
Initialize orientation object.
@ -263,70 +264,112 @@ class Orientation(Rotation):
raise TypeError('Use "O@b", i.e. matmul, to apply Orientation "O" to object "b"')
@staticmethod
def _split_kwargs(kwargs,target):
"""
Separate keyword arguments in 'kwargs' targeted at 'target' from general keyword arguments of Orientation objects.
Parameters
----------
kwargs : dictionary
Contains all **kwargs.
target: method
Function to scan for kwarg signature.
Returns
-------
rot_kwargs: dictionary
Valid keyword arguments of 'target' function of Rotation class.
ori_kwargs: dictionary
Valid keyword arguments of Orientation object.
"""
kws = ()
for t in (target,Orientation.__init__):
kws += ({key: kwargs[key] for key in set(inspect.signature(t).parameters) & set(kwargs)},)
invalid_keys = set(kwargs)-(set(kws[0])|set(kws[1]))
if invalid_keys:
raise TypeError(f"{inspect.stack()[1][3]}() got an unexpected keyword argument '{invalid_keys.pop()}'")
return kws
@classmethod
@util.extended_docstring(Rotation.from_random,_parameter_doc)
def from_random(cls,**kwargs):
return cls(rotation=Rotation.from_random(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_random)
return cls(rotation=Rotation.from_random(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_quaternion,_parameter_doc)
def from_quaternion(cls,**kwargs):
return cls(rotation=Rotation.from_quaternion(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_quaternion)
return cls(rotation=Rotation.from_quaternion(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_Euler_angles,_parameter_doc)
def from_Euler_angles(cls,**kwargs):
return cls(rotation=Rotation.from_Euler_angles(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Euler_angles)
return cls(rotation=Rotation.from_Euler_angles(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_axis_angle,_parameter_doc)
def from_axis_angle(cls,**kwargs):
return cls(rotation=Rotation.from_axis_angle(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_axis_angle)
return cls(rotation=Rotation.from_axis_angle(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_basis,_parameter_doc)
def from_basis(cls,**kwargs):
return cls(rotation=Rotation.from_basis(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_basis)
return cls(rotation=Rotation.from_basis(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_matrix,_parameter_doc)
def from_matrix(cls,**kwargs):
return cls(rotation=Rotation.from_matrix(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_matrix)
return cls(rotation=Rotation.from_matrix(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_Rodrigues_vector,_parameter_doc)
def from_Rodrigues_vector(cls,**kwargs):
return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_Rodrigues_vector)
return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_homochoric,_parameter_doc)
def from_homochoric(cls,**kwargs):
return cls(rotation=Rotation.from_homochoric(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_homochoric)
return cls(rotation=Rotation.from_homochoric(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_cubochoric,_parameter_doc)
def from_cubochoric(cls,**kwargs):
return cls(rotation=Rotation.from_cubochoric(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_cubochoric)
return cls(rotation=Rotation.from_cubochoric(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_spherical_component,_parameter_doc)
def from_spherical_component(cls,**kwargs):
return cls(rotation=Rotation.from_spherical_component(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_spherical_component)
return cls(rotation=Rotation.from_spherical_component(**kwargs_rot),**kwargs_ori)
@classmethod
@util.extended_docstring(Rotation.from_fiber_component,_parameter_doc)
def from_fiber_component(cls,**kwargs):
return cls(rotation=Rotation.from_fiber_component(**kwargs),**kwargs)
kwargs_rot,kwargs_ori = Orientation._split_kwargs(kwargs,Rotation.from_fiber_component)
return cls(rotation=Rotation.from_fiber_component(**kwargs_rot),**kwargs_ori)
@classmethod
@ -481,26 +524,26 @@ class Orientation(Rotation):
if self.family is None:
raise ValueError('Missing crystal symmetry')
rho_abs = np.abs(self.as_Rodrigues_vector(compact=True))
rho_abs = np.abs(self.as_Rodrigues_vector(compact=True))*(1.-1.e-9)
with np.errstate(invalid='ignore'):
# using '*'/prod for 'and'
if self.family == 'cubic':
return (np.prod(np.sqrt(2)-1. >= rho_abs,axis=-1) *
(1. >= np.sum(rho_abs,axis=-1))).astype(np.bool)
(1. >= np.sum(rho_abs,axis=-1))).astype(bool)
elif self.family == 'hexagonal':
return (np.prod(1. >= rho_abs,axis=-1) *
(2. >= np.sqrt(3)*rho_abs[...,0] + rho_abs[...,1]) *
(2. >= np.sqrt(3)*rho_abs[...,1] + rho_abs[...,0]) *
(2. >= np.sqrt(3) + rho_abs[...,2])).astype(np.bool)
(2. >= np.sqrt(3) + rho_abs[...,2])).astype(bool)
elif self.family == 'tetragonal':
return (np.prod(1. >= rho_abs[...,:2],axis=-1) *
(np.sqrt(2) >= rho_abs[...,0] + rho_abs[...,1]) *
(np.sqrt(2) >= rho_abs[...,2] + 1.)).astype(np.bool)
(np.sqrt(2) >= rho_abs[...,2] + 1.)).astype(bool)
elif self.family == 'orthorhombic':
return (np.prod(1. >= rho_abs,axis=-1)).astype(np.bool)
return (np.prod(1. >= rho_abs,axis=-1)).astype(bool)
elif self.family == 'monoclinic':
return (1. >= rho_abs[...,1]).astype(np.bool)
return (1. >= rho_abs[...,1]).astype(bool)
else:
return np.all(np.isfinite(rho_abs),axis=-1)
@ -524,28 +567,28 @@ class Orientation(Rotation):
if self.family is None:
raise ValueError('Missing crystal symmetry')
rho = self.as_Rodrigues_vector(compact=True)
rho = self.as_Rodrigues_vector(compact=True)*(1.0-1.0e-9)
with np.errstate(invalid='ignore'):
if self.family == 'cubic':
return ((rho[...,0] >= rho[...,1]) &
(rho[...,1] >= rho[...,2]) &
(rho[...,2] >= 0)).astype(np.bool)
(rho[...,2] >= 0)).astype(bool)
elif self.family == 'hexagonal':
return ((rho[...,0] >= rho[...,1]*np.sqrt(3)) &
(rho[...,1] >= 0) &
(rho[...,2] >= 0)).astype(np.bool)
(rho[...,2] >= 0)).astype(bool)
elif self.family == 'tetragonal':
return ((rho[...,0] >= rho[...,1]) &
(rho[...,1] >= 0) &
(rho[...,2] >= 0)).astype(np.bool)
(rho[...,2] >= 0)).astype(bool)
elif self.family == 'orthorhombic':
return ((rho[...,0] >= 0) &
(rho[...,1] >= 0) &
(rho[...,2] >= 0)).astype(np.bool)
(rho[...,2] >= 0)).astype(bool)
elif self.family == 'monoclinic':
return ((rho[...,1] >= 0) &
(rho[...,2] >= 0)).astype(np.bool)
(rho[...,2] >= 0)).astype(bool)
else:
return np.ones_like(rho[...,0],dtype=bool)

View File

@ -502,8 +502,8 @@ class Rotation:
Returns
-------
c : numpy.ndarray of shape (...,3)
Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3).
x : numpy.ndarray of shape (...,3)
Cubochoric vector: (x_1, x_2, x_3), max(x_i) < 1/2*π^(2/3).
"""
return Rotation._qu2cu(self.quaternion)
@ -514,8 +514,7 @@ class Rotation:
@staticmethod
def from_quaternion(q,
accept_homomorph = False,
P = -1,
**kwargs):
P = -1):
"""
Initialize from quaternion.
@ -550,8 +549,7 @@ class Rotation:
@staticmethod
def from_Euler_angles(phi,
degrees = False,
**kwargs):
degrees = False):
"""
Initialize from Bunge-Euler angles.
@ -578,8 +576,7 @@ class Rotation:
def from_axis_angle(axis_angle,
degrees = False,
normalize = False,
P = -1,
**kwargs):
P = -1):
"""
Initialize from Axis angle pair.
@ -616,8 +613,7 @@ class Rotation:
@staticmethod
def from_basis(basis,
orthonormal = True,
reciprocal = False,
**kwargs):
reciprocal = False):
"""
Initialize from lattice basis vectors.
@ -651,7 +647,7 @@ class Rotation:
return Rotation(Rotation._om2qu(om))
@staticmethod
def from_matrix(R,**kwargs):
def from_matrix(R):
"""
Initialize from rotation matrix.
@ -695,8 +691,7 @@ class Rotation:
@staticmethod
def from_Rodrigues_vector(rho,
normalize = False,
P = -1,
**kwargs):
P = -1):
"""
Initialize from Rodrigues-Frank vector (angle separated from axis).
@ -727,8 +722,7 @@ class Rotation:
@staticmethod
def from_homochoric(h,
P = -1,
**kwargs):
P = -1):
"""
Initialize from homochoric vector.
@ -754,21 +748,20 @@ class Rotation:
return Rotation(Rotation._ho2qu(ho))
@staticmethod
def from_cubochoric(c,
P = -1,
**kwargs):
def from_cubochoric(x,
P = -1):
"""
Initialize from cubochoric vector.
Parameters
----------
c : numpy.ndarray of shape (...,3)
Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3).
x : numpy.ndarray of shape (...,3)
Cubochoric vector: (x_1, x_2, x_3), max(x_i) < 1/2*π^(2/3).
P : int {-1,1}, optional
Convention used. Defaults to -1.
"""
cu = np.array(c,dtype=float)
cu = np.array(x,dtype=float)
if cu.shape[:-2:-1] != (3,):
raise ValueError('Invalid shape.')
if abs(P) != 1:
@ -784,8 +777,7 @@ class Rotation:
@staticmethod
def from_random(shape = None,
rng_seed = None,
**kwargs):
rng_seed = None):
"""
Draw random rotation.
@ -878,8 +870,7 @@ class Rotation:
sigma,
N = 500,
degrees = True,
rng_seed = None,
**kwargs):
rng_seed = None):
"""
Calculate set of rotations with Gaussian distribution around center.
@ -915,8 +906,7 @@ class Rotation:
sigma = 0.0,
N = 500,
degrees = True,
rng_seed = None,
**kwargs):
rng_seed = None):
"""
Calculate set of rotations with Gaussian distribution around direction.

View File

@ -7,6 +7,7 @@ from damask import Orientation
from damask import Table
from damask import lattice
from damask import util
from damask import grid_filters
@pytest.fixture
@ -118,7 +119,7 @@ class TestOrientation:
== np.eye(3))
def test_from_cubochoric(self):
assert np.all(Orientation.from_cubochoric(c=np.zeros(3),lattice='triclinic').as_matrix()
assert np.all(Orientation.from_cubochoric(x=np.zeros(3),lattice='triclinic').as_matrix()
== np.eye(3))
def test_from_spherical_component(self):
@ -141,7 +142,7 @@ class TestOrientation:
dict(lattice='hP',a=1.0 ),
dict(lattice='cI',a=1.0, ),
])
def test_from_direction(self,kwargs):
def test_from_directions(self,kwargs):
for a,b in np.random.random((10,2,3)):
c = np.cross(b,a)
if np.all(np.isclose(c,0)): continue
@ -151,6 +152,21 @@ class TestOrientation:
assert np.isclose(np.dot(x/np.linalg.norm(x),np.array([1,0,0])),1) \
and np.isclose(np.dot(z/np.linalg.norm(z),np.array([0,0,1])),1)
@pytest.mark.parametrize('function',[Orientation.from_random,
Orientation.from_quaternion,
Orientation.from_Euler_angles,
Orientation.from_axis_angle,
Orientation.from_basis,
Orientation.from_matrix,
Orientation.from_Rodrigues_vector,
Orientation.from_homochoric,
Orientation.from_cubochoric,
Orientation.from_spherical_component,
Orientation.from_fiber_component,
Orientation.from_directions])
def test_invalid_from(self,function):
with pytest.raises(TypeError):
function(c=.1,degrees=True,invalid=66)
def test_negative_angle(self):
with pytest.raises(ValueError):
@ -221,6 +237,16 @@ class TestOrientation:
for r, theO in zip(o.reduced.flatten(),o.flatten()):
assert r == theO.reduced
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
def test_reduced_corner_cases(self,lattice):
# test whether there is always a sym-eq rotation that falls into the FZ
N = np.random.randint(10,40)
size = np.ones(3)*np.pi**(2./3.)
grid = grid_filters.coordinates0_node([N+1,N+1,N+1],size,-size*.5)
evenly_distributed = Orientation.from_cubochoric(x=grid[:-2,:-2,:-2],lattice=lattice)
assert evenly_distributed.shape == evenly_distributed.reduced.shape
@pytest.mark.parametrize('lattice',Orientation.crystal_families)
@pytest.mark.parametrize('shape',[(1),(2,3),(4,3,2)])
@pytest.mark.parametrize('vector',np.array([[1,0,0],[1,2,3],[-1,1,-1]]))

View File

@ -78,7 +78,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
p, i, &
ph, i, &
Nconstituents, &
sizeState, sizeDotState, &
startIndex, endIndex
@ -114,15 +114,13 @@ module function plastic_dislotungsten_init() result(myPlasticity)
allocate(dependentState(phases%length))
do p = 1, phases%length
if(.not. myPlasticity(p)) cycle
phase => phases%get(p)
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
i = p
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i), &
dst => dependentState(i))
pl => mech%get('plasticity')
#if defined (__GFORTRAN__)
@ -132,7 +130,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
#endif
! This data is read in already in lattice
prm%mu = lattice_mu(p)
prm%mu = lattice_mu(ph)
!--------------------------------------------------------------------------------------------------
! slip related parameters
@ -222,41 +220,41 @@ module function plastic_dislotungsten_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob = spread(rho_mob_0,2,Nconstituents)
dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
dot%rho_mob => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip => plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip = spread(rho_dip_0,2,Nconstituents)
dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
dot%rho_dip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl => plasticState(p)%state(startIndex:endIndex,:)
dot%gamma_sl => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = 1.0e-2_pReal
stt%gamma_sl => plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents), source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate

View File

@ -126,7 +126,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
p, i, &
ph, i, &
Nconstituents, &
sizeState, sizeDotState, &
startIndex, endIndex
@ -167,15 +167,13 @@ module function plastic_dislotwin_init() result(myPlasticity)
allocate(dependentState(phases%length))
do p = 1, phases%length
if(.not. myPlasticity(p)) cycle
phase => phases%get(p)
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), dst => dependentState(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
i = p
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i), &
dst => dependentState(i))
pl => mech%get('plasticity')
#if defined (__GFORTRAN__)
@ -185,9 +183,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
#endif
! This data is read in already in lattice
prm%mu = lattice_mu(p)
prm%nu = lattice_nu(p)
prm%C66 = lattice_C66(1:6,1:6,p)
prm%mu = lattice_mu(ph)
prm%nu = lattice_nu(ph)
prm%C66 = lattice_C66(1:6,1:6,ph)
!--------------------------------------------------------------------------------------------------
! slip related parameters
@ -204,8 +202,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
prm%n0_sl = lattice_slip_normal(N_sl,phase%get_asString('lattice'),&
phase%get_asFloat('c/a',defaultVal=0.0_pReal))
prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == lattice_FCC_ID) &
.and. (N_sl(1) == 12)
prm%fccTwinTransNucleation = lattice_structure(ph) == lattice_FCC_ID .and. (N_sl(1) == 12)
if(prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR
rho_mob_0 = pl%get_asFloats('rho_mob_0', requiredSize=size(N_sl))
@ -234,7 +231,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex)
! details: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981
prm%omega = pl%get_asFloat('omega', defaultVal = 1000.0_pReal) &
* merge(12.0_pReal,8.0_pReal,any(lattice_structure(p) == [lattice_FCC_ID,lattice_HEX_ID]))
* merge(12.0_pReal,8.0_pReal,any(lattice_structure(ph) == [lattice_FCC_ID,lattice_HEX_ID]))
! expand: family => system
rho_mob_0 = math_expand(rho_mob_0, N_sl)
@ -342,7 +339,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
pl%get_asFloat('a_cI', defaultVal=0.0_pReal), &
pl%get_asFloat('a_cF', defaultVal=0.0_pReal))
if (lattice_structure(p) /= lattice_FCC_ID) then
if (lattice_structure(ph) /= lattice_FCC_ID) then
prm%dot_N_0_tr = pl%get_asFloats('dot_N_0_tr')
prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,N_tr)
endif
@ -357,7 +354,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
if ( prm%i_tr < 0.0_pReal) extmsg = trim(extmsg)//' i_tr'
if (any(prm%t_tr < 0.0_pReal)) extmsg = trim(extmsg)//' t_tr'
if (any(prm%s < 0.0_pReal)) extmsg = trim(extmsg)//' p_tr'
if (lattice_structure(p) /= lattice_FCC_ID) then
if (lattice_structure(ph) /= lattice_FCC_ID) then
if (any(prm%dot_N_0_tr < 0.0_pReal)) extmsg = trim(extmsg)//' dot_N_0_tr'
endif
else transActive
@ -408,53 +405,53 @@ module function plastic_dislotwin_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol
startIndex = 1
endIndex = prm%sum_N_sl
stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob=>plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_mob= spread(rho_mob_0,2,Nconstituents)
dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
dot%rho_mob=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip=>plasticState(ph)%state(startIndex:endIndex,:)
stt%rho_dip= spread(rho_dip_0,2,Nconstituents)
dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
dot%rho_dip=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:)
dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = 1.0e-2_pReal
stt%gamma_sl=>plasticState(ph)%state(startIndex:endIndex,:)
dot%gamma_sl=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = 1.0e-2_pReal
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%f_tw=>plasticState(p)%state(startIndex:endIndex,:)
dot%f_tw=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin'
stt%f_tw=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tw=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_twin',defaultVal=1.0e-7_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_twin'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tr
stt%f_tr=>plasticState(p)%state(startIndex:endIndex,:)
dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans'
stt%f_tr=>plasticState(ph)%state(startIndex:endIndex,:)
dot%f_tr=>plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal)
if (any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans'
allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
@ -469,7 +466,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
allocate(dst%tau_r_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
allocate(dst%V_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate

View File

@ -22,8 +22,6 @@ submodule(phase:plastic) isotropic
c_4, &
c_3, &
c_2
integer :: &
of_debug = 0
logical :: &
dilatation
character(len=pStringLen), allocatable, dimension(:) :: &
@ -53,8 +51,7 @@ module function plastic_isotropic_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
p, &
i, &
ph, &
Nconstituents, &
sizeState, sizeDotState
real(pReal) :: &
@ -82,16 +79,14 @@ module function plastic_isotropic_init() result(myPlasticity)
allocate(state(phases%length))
allocate(dotState(phases%length))
do p = 1, phases%length
if(.not. myPlasticity(p)) cycle
phase => phases%get(p)
mech => phase%get('mechanics')
i = p
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i))
pl => mech%get('plasticity')
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
pl => mech%get('plasticity')
#if defined (__GFORTRAN__)
prm%output = output_asStrings(pl)
@ -99,11 +94,6 @@ module function plastic_isotropic_init() result(myPlasticity)
prm%output = pl%get_asStrings('output',defaultVal=emptyStringArray)
#endif
#ifdef DEBUG
if (p==material_phaseAt(debugConstitutive%grain,debugConstitutive%element)) &
prm%of_debug = material_phasememberAt(debugConstitutive%grain,debugConstitutive%ip,debugConstitutive%element)
#endif
xi_0 = pl%get_asFloat('xi_0')
prm%xi_inf = pl%get_asFloat('xi_inf')
prm%dot_gamma_0 = pl%get_asFloat('dot_gamma_0')
@ -129,28 +119,28 @@ module function plastic_isotropic_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
stt%xi => plasticState(p)%state (1,:)
stt%xi => plasticState(ph)%state (1,:)
stt%xi = xi_0
dot%xi => plasticState(p)%dotState(1,:)
plasticState(p)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if (plasticState(p)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
dot%xi => plasticState(ph)%dotState(1,:)
plasticState(ph)%atol(1) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if (plasticState(ph)%atol(1) < 0.0_pReal) extmsg = trim(extmsg)//' atol_xi'
stt%gamma => plasticState(p)%state (2,:)
dot%gamma => plasticState(p)%dotState(2,:)
plasticState(p)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if (plasticState(p)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
stt%gamma => plasticState(ph)%state (2,:)
dot%gamma => plasticState(ph)%dotState(2,:)
plasticState(ph)%atol(2) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if (plasticState(ph)%atol(2) < 0.0_pReal) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(2:2,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(2:2,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate
@ -240,13 +230,11 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me)
tr=math_trace33(math_spherical33(Mi))
if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero
if (prm%dilatation .and. abs(tr) > 0.0_pReal) then ! no stress or J2 plasticity --> Li and its derivative are zero
Li = math_I3 &
* prm%dot_gamma_0/prm%M * (3.0_pReal*prm%M*stt%xi(me))**(-prm%n) &
* tr * abs(tr)**(prm%n-1.0_pReal)
forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
else
Li = 0.0_pReal
dLi_dMi = 0.0_pReal

View File

@ -25,8 +25,7 @@ submodule(phase:plastic) kinehardening
nonSchmid_pos, &
nonSchmid_neg
integer :: &
sum_N_sl, & !< total number of active slip system
of_debug = 0
sum_N_sl
logical :: &
nonSchmidActive = .false.
character(len=pStringLen), allocatable, dimension(:) :: &
@ -62,7 +61,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
p, i, o, &
ph, o, &
Nconstituents, &
sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex
@ -93,15 +92,13 @@ module function plastic_kinehardening_init() result(myPlasticity)
allocate(deltaState(phases%length))
do p = 1, phases%length
if(.not. myPlasticity(p)) cycle
phase => phases%get(p)
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), dlt => deltaState(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
i = p
associate(prm => param(i), &
dot => dotState(i), &
dlt => deltaState(i), &
stt => state(i))
pl => mech%get('plasticity')
#if defined (__GFORTRAN__)
@ -110,12 +107,6 @@ module function plastic_kinehardening_init() result(myPlasticity)
prm%output = pl%get_asStrings('output',defaultVal=emptyStringArray)
#endif
#ifdef DEBUG
if (p==material_phaseAt(debugConstitutive%grain,debugConstitutive%element)) then
prm%of_debug = material_phasememberAt(debugConstitutive%grain,debugConstitutive%ip,debugConstitutive%element)
endif
#endif
!--------------------------------------------------------------------------------------------------
! slip related parameters
N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
@ -174,55 +165,55 @@ module function plastic_kinehardening_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeState = sizeDotState + sizeDeltaState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss => plasticState(ph)%state (startIndex:endIndex,:)
stt%crss = spread(xi_0, 2, Nconstituents)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%crss_back => plasticState(p)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
stt%crss_back => plasticState(ph)%state (startIndex:endIndex,:)
dot%crss_back => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%accshear => plasticState(p)%state (startIndex:endIndex,:)
dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
stt%accshear => plasticState(ph)%state (startIndex:endIndex,:)
dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
o = plasticState(p)%offsetDeltaState
o = plasticState(ph)%offsetDeltaState
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%sense => plasticState(p)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
stt%sense => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%sense => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
stt%chi0 => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%chi0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:)
stt%gamma0 => plasticState(ph)%state (startIndex :endIndex ,:)
dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate

View File

@ -16,8 +16,7 @@ module function plastic_none_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
p, &
Nconstituents
ph
class(tNode), pointer :: &
phases
@ -29,10 +28,9 @@ module function plastic_none_init() result(myPlasticity)
print'(a,i0)', ' # phases: ',count(myPlasticity); flush(IO_STDOUT)
phases => config_material%get('phase')
do p = 1, phases%length
if(.not. myPlasticity(p)) cycle
Nconstituents = count(material_phaseAt2 == p)
call phase_allocateState(plasticState(p),Nconstituents,0,0,0)
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
call phase_allocateState(plasticState(ph),count(material_phaseAt2 == ph),0,0,0)
enddo
end function plastic_none_init

View File

@ -176,7 +176,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
Ninstances, &
p, i, &
ph, &
Nconstituents, &
sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
s1, s2, &
@ -220,21 +220,17 @@ module function plastic_nonlocal_init() result(myPlasticity)
allocate(deltaState(phases%length))
allocate(microstructure(phases%length))
do p = 1, phases%length
if(.not. myPlasticity(p)) cycle
phase => phases%get(p)
mech => phase%get('mechanics')
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
i = p
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i), &
st0 => state0(i), &
del => deltaState(i), &
dst => microstructure(i))
associate(prm => param(ph), dot => dotState(ph), stt => state(ph), &
st0 => state0(ph), del => deltaState(ph), dst => microstructure(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
pl => mech%get('plasticity')
phase_localPlasticity(p) = .not. pl%contains('nonlocal')
phase_localPlasticity(ph) = .not. pl%contains('nonlocal')
#if defined (__GFORTRAN__)
prm%output = output_asStrings(pl)
@ -245,8 +241,8 @@ module function plastic_nonlocal_init() result(myPlasticity)
prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0e4_pReal)
! This data is read in already in lattice
prm%mu = lattice_mu(p)
prm%nu = lattice_nu(p)
prm%mu = lattice_mu(ph)
prm%nu = lattice_nu(ph)
ini%N_sl = pl%get_asInts('N_sl',defaultVal=emptyIntArray)
prm%sum_N_sl = sum(abs(ini%N_sl))
@ -402,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -416,101 +412,101 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
allocate(geom(p)%V_0(Nconstituents))
call storeGeometry(p)
allocate(geom(ph)%V_0(Nconstituents))
call storeGeometry(ph)
plasticState(p)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) &
plasticState(ph)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(ph)%nonlocal .and. .not. allocated(IPneighborhood)) &
call IO_error(212,ext_msg='IPneighborhood does not exist')
plasticState(p)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
plasticState(ph)%offsetDeltaState = 0 ! ToDo: state structure does not follow convention
st0%rho => plasticState(p)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho => plasticState(p)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho => plasticState(p)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho => plasticState(p)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
plasticState(p)%atol(1:10*prm%sum_N_sl) = prm%atol_rho
st0%rho => plasticState(ph)%state0 (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho => plasticState(ph)%state (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho => plasticState(ph)%dotState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho => plasticState(ph)%deltaState (0*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
plasticState(ph)%atol(1:10*prm%sum_N_sl) = prm%atol_rho
stt%rhoSgl => plasticState(p)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSgl => plasticState(p)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSgl => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoSgl => plasticState(ph)%state (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSgl => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSgl => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoSglMobile => plasticState(p)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rhoSglMobile => plasticState(p)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rhoSglMobile => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rhoSglMobile => plasticState(ph)%state (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rhoSglMobile => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rhoSglMobile => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_pos => plasticState(p)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_pos => plasticState(p)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_pos => plasticState(p)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_pos => plasticState(ph)%state (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_pos => plasticState(ph)%dotState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_pos => plasticState(ph)%deltaState (0*prm%sum_N_sl+1: 1*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_neg => plasticState(p)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_neg => plasticState(p)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_neg => plasticState(p)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
stt%rho_sgl_mob_edg_neg => plasticState(ph)%state (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
dot%rho_sgl_mob_edg_neg => plasticState(ph)%dotState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
del%rho_sgl_mob_edg_neg => plasticState(ph)%deltaState (1*prm%sum_N_sl+1: 2*prm%sum_N_sl,:)
stt%rho_sgl_mob_scr_pos => plasticState(p)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_pos => plasticState(p)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_pos => plasticState(p)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
stt%rho_sgl_mob_scr_pos => plasticState(ph)%state (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_pos => plasticState(ph)%dotState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_pos => plasticState(ph)%deltaState (2*prm%sum_N_sl+1: 3*prm%sum_N_sl,:)
stt%rho_sgl_mob_scr_neg => plasticState(p)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_neg => plasticState(p)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_neg => plasticState(p)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rho_sgl_mob_scr_neg => plasticState(ph)%state (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
dot%rho_sgl_mob_scr_neg => plasticState(ph)%dotState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
del%rho_sgl_mob_scr_neg => plasticState(ph)%deltaState (3*prm%sum_N_sl+1: 4*prm%sum_N_sl,:)
stt%rhoSglImmobile => plasticState(p)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSglImmobile => plasticState(p)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSglImmobile => plasticState(p)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoSglImmobile => plasticState(ph)%state (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rhoSglImmobile => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rhoSglImmobile => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_pos => plasticState(p)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_pos => plasticState(p)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_pos => plasticState(p)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_pos => plasticState(ph)%state (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_pos => plasticState(ph)%dotState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_pos => plasticState(ph)%deltaState (4*prm%sum_N_sl+1: 5*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_neg => plasticState(p)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_neg => plasticState(p)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_neg => plasticState(p)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
stt%rho_sgl_imm_edg_neg => plasticState(ph)%state (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
dot%rho_sgl_imm_edg_neg => plasticState(ph)%dotState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
del%rho_sgl_imm_edg_neg => plasticState(ph)%deltaState (5*prm%sum_N_sl+1: 6*prm%sum_N_sl,:)
stt%rho_sgl_imm_scr_pos => plasticState(p)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_pos => plasticState(p)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_pos => plasticState(p)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
stt%rho_sgl_imm_scr_pos => plasticState(ph)%state (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_pos => plasticState(ph)%dotState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_pos => plasticState(ph)%deltaState (6*prm%sum_N_sl+1: 7*prm%sum_N_sl,:)
stt%rho_sgl_imm_scr_neg => plasticState(p)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_neg => plasticState(p)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_neg => plasticState(p)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rho_sgl_imm_scr_neg => plasticState(ph)%state (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
dot%rho_sgl_imm_scr_neg => plasticState(ph)%dotState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
del%rho_sgl_imm_scr_neg => plasticState(ph)%deltaState (7*prm%sum_N_sl+1: 8*prm%sum_N_sl,:)
stt%rhoDip => plasticState(p)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rhoDip => plasticState(p)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rhoDip => plasticState(p)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rhoDip => plasticState(ph)%state (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rhoDip => plasticState(ph)%dotState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rhoDip => plasticState(ph)%deltaState (8*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho_dip_edg => plasticState(p)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
dot%rho_dip_edg => plasticState(p)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
del%rho_dip_edg => plasticState(p)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
stt%rho_dip_edg => plasticState(ph)%state (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
dot%rho_dip_edg => plasticState(ph)%dotState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
del%rho_dip_edg => plasticState(ph)%deltaState (8*prm%sum_N_sl+1: 9*prm%sum_N_sl,:)
stt%rho_dip_scr => plasticState(p)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho_dip_scr => plasticState(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho_dip_scr => plasticState(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%rho_dip_scr => plasticState(ph)%state (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
dot%rho_dip_scr => plasticState(ph)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho_dip_scr => plasticState(ph)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%gamma => plasticState(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
dot%gamma => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
del%gamma => plasticState(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
if(any(plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
stt%gamma => plasticState(ph)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
dot%gamma => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
del%gamma => plasticState(ph)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
if(any(plasticState(ph)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
plasticState(ph)%slipRate => plasticState(ph)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
stt%rho_forest => plasticState(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents)
stt%v => plasticState(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_pos => plasticState(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_neg => plasticState(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_pos => plasticState(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_neg => plasticState(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
stt%rho_forest => plasticState(ph)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents)
stt%v => plasticState(ph)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_pos => plasticState(ph)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_neg => plasticState(ph)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_pos => plasticState(ph)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_neg => plasticState(ph)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
end associate
if (Nconstituents > 0) call stateInit(ini,p,Nconstituents)
plasticState(p)%state0 = plasticState(p)%state
if (Nconstituents > 0) call stateInit(ini,ph,Nconstituents)
plasticState(ph)%state0 = plasticState(ph)%state
!--------------------------------------------------------------------------------------------------
! exit if any parameter is out of range
@ -526,35 +522,34 @@ module function plastic_nonlocal_init() result(myPlasticity)
allocate(iV(maxval(param%sum_N_sl),4,phases%length), source=0)
allocate(iD(maxval(param%sum_N_sl),2,phases%length), source=0)
do p = 1, phases%length
phase => phases%get(p)
do ph = 1, phases%length
if(.not. myPlasticity(p)) cycle
i = p
if(.not. myPlasticity(ph)) cycle
Nconstituents = count(material_phaseAt2 == p)
phase => phases%get(ph)
Nconstituents = count(material_phaseAt2 == ph)
l = 0
do t = 1,4
do s = 1,param(i)%sum_N_sl
do s = 1,param(ph)%sum_N_sl
l = l + 1
iRhoU(s,t,i) = l
iRhoU(s,t,ph) = l
enddo
enddo
l = l + (4+2+1+1)*param(i)%sum_N_sl ! immobile(4), dipole(2), shear, forest
l = l + (4+2+1+1)*param(ph)%sum_N_sl ! immobile(4), dipole(2), shear, forest
do t = 1,4
do s = 1,param(i)%sum_N_sl
do s = 1,param(ph)%sum_N_sl
l = l + 1
iV(s,t,i) = l
iV(s,t,ph) = l
enddo
enddo
do t = 1,2
do s = 1,param(i)%sum_N_sl
do s = 1,param(ph)%sum_N_sl
l = l + 1
iD(s,t,i) = l
iD(s,t,ph) = l
enddo
enddo
if (iD(param(i)%sum_N_sl,2,i) /= plasticState(p)%sizeState) &
call IO_error(0, ext_msg = 'state indices not properly set (nonlocal)')
if (iD(param(ph)%sum_N_sl,2,ph) /= plasticState(ph)%sizeState) &
error stop 'state indices not properly set (nonlocal)'
enddo
end function plastic_nonlocal_init

View File

@ -70,7 +70,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity
integer :: &
p, i, &
ph, i, &
Nconstituents, &
sizeState, sizeDotState, &
startIndex, endIndex
@ -101,14 +101,13 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
allocate(state(phases%length))
allocate(dotState(phases%length))
do p = 1, phases%length
if(.not. myPlasticity(p)) cycle
phase => phases%get(p)
do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle
associate(prm => param(ph), dot => dotState(ph), stt => state(ph))
phase => phases%get(ph)
mech => phase%get('mechanics')
i = p
associate(prm => param(i), &
dot => dotState(i), &
stt => state(i))
pl => mech%get('plasticity')
!--------------------------------------------------------------------------------------------------
@ -135,7 +134,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
xi_0_sl = pl%get_asFloats('xi_0_sl', requiredSize=size(N_sl))
prm%xi_inf_sl = pl%get_asFloats('xi_inf_sl', requiredSize=size(N_sl))
prm%h_int = pl%get_asFloats('h_int', requiredSize=size(N_sl), &
defaultVal=[(0.0_pReal,i=1,size(N_sl))])
defaultVal=[(0.0_pReal,i=1,size(N_sl))])
prm%dot_gamma_0_sl = pl%get_asFloat('dot_gamma_0_sl')
prm%n_sl = pl%get_asFloat('n_sl')
@ -224,49 +223,49 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!--------------------------------------------------------------------------------------------------
! allocate state arrays
Nconstituents = count(material_phaseAt2 == p)
Nconstituents = count(material_phaseAt2 == ph)
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState
call phase_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
call phase_allocateState(plasticState(ph),Nconstituents,sizeState,sizeDotState,0)
!--------------------------------------------------------------------------------------------------
! state aliases and initialization
startIndex = 1
endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:)
stt%xi_slip = spread(xi_0_sl, 2, Nconstituents)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:)
stt%xi_twin = spread(xi_0_tw, 2, Nconstituents)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
stt%gamma_slip => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(ph)%slipRate => plasticState(ph)%dotState(startIndex:endIndex,:)
startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw
stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:)
dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally
plasticState(ph)%state0 = plasticState(ph)%state ! ToDo: this could be done centrally
end associate