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

Merge branch 'development' into misc-improvements

This commit is contained in:
Martin Diehl 2020-04-15 21:03:51 +02:00
parent 69857176b2 875b3b6321
commit e3958263e3
8 changed files with 776 additions and 408 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
PRIVATE

@ -1 +1 @@
Subproject commit aa0b9fe992ce0bfc172989ab006893ddb8939c1e Subproject commit 232a094c715bcbbd1c6652c4dc4a4a50d402b82f

2
VERSION
View File

@ -1 +1 @@
v2.0.3-2255-gc74ffae8 v2.0.3-2303-g2a6132b7

34
python/damask/_Lambert.py
View File

@ -51,19 +51,18 @@ def cube_to_ball(cube):
https://doi.org/10.1088/0965-0393/22/7/075013 https://doi.org/10.1088/0965-0393/22/7/075013
""" """
if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5: cube_ = np.clip(cube,None,np.pi**(2./3.) * 0.5) if np.isclose(np.abs(np.max(cube)),np.pi**(2./3.) * 0.5,atol=1e-6) else cube
raise ValueError
# transform to the sphere grid via the curved square, and intercept the zero point # transform to the sphere grid via the curved square, and intercept the zero point
if np.allclose(cube,0.0,rtol=0.0,atol=1.0e-300): if np.allclose(cube_,0.0,rtol=0.0,atol=1.0e-16):
ball = np.zeros(3) ball = np.zeros(3)
else: else:
# get pyramide and scale by grid parameter ratio # get pyramide and scale by grid parameter ratio
p = _get_order(cube) p = _get_order(cube_)
XYZ = cube[p] * sc XYZ = cube_[p[0]] * sc
# intercept all the points along the z-axis # intercept all the points along the z-axis
if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-300): if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-16):
ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]]) ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]])
else: else:
order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1] order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1]
@ -82,7 +81,7 @@ def cube_to_ball(cube):
ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ]) ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ])
# reverse the coordinates back to the regular order according to the original pyramid number # reverse the coordinates back to the regular order according to the original pyramid number
ball = ball[p] ball = ball[p[1]]
return ball return ball
@ -102,15 +101,14 @@ def ball_to_cube(ball):
https://doi.org/10.1088/0965-0393/22/7/075013 https://doi.org/10.1088/0965-0393/22/7/075013
""" """
rs = np.linalg.norm(ball) ball_ = ball/np.linalg.norm(ball)*R1 if np.isclose(np.linalg.norm(ball),R1,atol=1e-6) else ball
if rs > R1: rs = np.linalg.norm(ball_)
raise ValueError
if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300): if np.allclose(ball_,0.0,rtol=0.0,atol=1.0e-16):
cube = np.zeros(3) cube = np.zeros(3)
else: else:
p = _get_order(ball) p = _get_order(ball_)
xyz3 = ball[p] xyz3 = ball_[p[0]]
# inverse M_3 # inverse M_3
xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) ) xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) )
@ -118,7 +116,7 @@ def ball_to_cube(ball):
# inverse M_2 # inverse M_2
qxy = np.sum(xyz2**2) qxy = np.sum(xyz2**2)
if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-300): if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-16):
Tinv = np.zeros(2) Tinv = np.zeros(2)
else: else:
q2 = qxy + np.max(np.abs(xyz2))**2 q2 = qxy + np.max(np.abs(xyz2))**2
@ -132,7 +130,7 @@ def ball_to_cube(ball):
# inverse M_1 # inverse M_1
cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc
# reverse the coordinates back to the regular order according to the original pyramid number # reverse the coordinates back to the regular order according to the original pyramid number
cube = cube[p] cube = cube[p[1]]
return cube return cube
@ -157,10 +155,10 @@ def _get_order(xyz):
""" """
if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \ if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \
(abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]): (abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]):
return [0,1,2] return [[0,1,2],[0,1,2]]
elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \ elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \
(abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]): (abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]):
return [1,2,0] return [[1,2,0],[2,0,1]]
elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \ elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \
(abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]): (abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]):
return [2,0,1] return [[2,0,1],[1,2,0]]

530
python/damask/_rotation.py
View File

@ -2,7 +2,7 @@ import numpy as np
from ._Lambert import ball_to_cube, cube_to_ball from ._Lambert import ball_to_cube, cube_to_ball
P = -1 _P = -1
def iszero(a): def iszero(a):
return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0)
@ -89,7 +89,7 @@ class Rotation:
other_q = other.quaternion[0] other_q = other.quaternion[0]
other_p = other.quaternion[1:] other_p = other.quaternion[1:]
R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p), R = self.__class__(np.append(self_q*other_q - np.dot(self_p,other_p),
self_q*other_p + other_q*self_p + P * np.cross(self_p,other_p))) self_q*other_p + other_q*self_p + _P * np.cross(self_p,other_p)))
return R.standardize() return R.standardize()
elif isinstance(other, (tuple,np.ndarray)): elif isinstance(other, (tuple,np.ndarray)):
if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid if isinstance(other,tuple) or other.shape == (3,): # rotate a single (3)-vector or meshgrid
@ -97,7 +97,7 @@ class Rotation:
B = 2.0 * ( self.quaternion[1]*other[0] B = 2.0 * ( self.quaternion[1]*other[0]
+ self.quaternion[2]*other[1] + self.quaternion[2]*other[1]
+ self.quaternion[3]*other[2]) + self.quaternion[3]*other[2])
C = 2.0 * P*self.quaternion[0] C = 2.0 * _P*self.quaternion[0]
return np.array([ return np.array([
A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]), A*other[0] + B*self.quaternion[1] + C*(self.quaternion[2]*other[2] - self.quaternion[3]*other[1]),
@ -206,7 +206,7 @@ class Rotation:
""" """
ax = Rotation.qu2ax(self.quaternion) ax = Rotation.qu2ax(self.quaternion)
if degrees: ax[3] = np.degrees(ax[3]) if degrees: ax[3] = np.degrees(ax[3])
return (ax[:3],np.degrees(ax[3])) if pair else ax return (ax[:3],ax[3]) if pair else ax
def asMatrix(self): def asMatrix(self):
"""Rotation matrix.""" """Rotation matrix."""
@ -262,9 +262,9 @@ class Rotation:
if acceptHomomorph: if acceptHomomorph:
qu *= -1. qu *= -1.
else: else:
raise ValueError('Quaternion has negative first component.\n{}'.format(qu[0])) raise ValueError('Quaternion has negative first component: {}.'.format(qu[0]))
if not np.isclose(np.linalg.norm(qu), 1.0): if not np.isclose(np.linalg.norm(qu), 1.0):
raise ValueError('Quaternion is not of unit length.\n{} {} {} {}'.format(*qu)) raise ValueError('Quaternion is not of unit length: {} {} {} {}.'.format(*qu))
return Rotation(qu) return Rotation(qu)
@ -276,7 +276,7 @@ class Rotation:
else np.array(eulers,dtype=float) else np.array(eulers,dtype=float)
eu = np.radians(eu) if degrees else eu eu = np.radians(eu) if degrees else eu
if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi: if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi:
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eu)) raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]: {} {} {}.'.format(*eu))
return Rotation(Rotation.eu2qu(eu)) return Rotation(Rotation.eu2qu(eu))
@ -292,9 +292,9 @@ class Rotation:
if degrees: ax[ 3] = np.radians(ax[3]) if degrees: ax[ 3] = np.radians(ax[3])
if normalise: ax[0:3] /= np.linalg.norm(ax[0:3]) if normalise: ax[0:3] /= np.linalg.norm(ax[0:3])
if ax[3] < 0.0 or ax[3] > np.pi: if ax[3] < 0.0 or ax[3] > np.pi:
raise ValueError('Axis angle rotation angle outside of [0..π].\n{}'.format(ax[3])) raise ValueError('Axis angle rotation angle outside of [0..π]: {}.'.format(ax[3]))
if not np.isclose(np.linalg.norm(ax[0:3]), 1.0): if not np.isclose(np.linalg.norm(ax[0:3]), 1.0):
raise ValueError('Axis angle rotation axis is not of unit length.\n{} {} {}'.format(*ax[0:3])) raise ValueError('Axis angle rotation axis is not of unit length: {} {} {}.'.format(*ax[0:3]))
return Rotation(Rotation.ax2qu(ax)) return Rotation(Rotation.ax2qu(ax))
@ -312,11 +312,11 @@ class Rotation:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition (U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.dot(U,Vh) om = np.dot(U,Vh)
if not np.isclose(np.linalg.det(om),1.0): if not np.isclose(np.linalg.det(om),1.0):
raise ValueError('matrix is not a proper rotation.\n{}'.format(om)) raise ValueError('matrix is not a proper rotation: {}.'.format(om))
if not np.isclose(np.dot(om[0],om[1]), 0.0) \ if not np.isclose(np.dot(om[0],om[1]), 0.0) \
or not np.isclose(np.dot(om[1],om[2]), 0.0) \ or not np.isclose(np.dot(om[1],om[2]), 0.0) \
or not np.isclose(np.dot(om[2],om[0]), 0.0): or not np.isclose(np.dot(om[2],om[0]), 0.0):
raise ValueError('matrix is not orthogonal.\n{}'.format(om)) raise ValueError('matrix is not orthogonal: {}.'.format(om))
return Rotation(Rotation.om2qu(om)) return Rotation(Rotation.om2qu(om))
@ -336,9 +336,9 @@ class Rotation:
if P > 0: ro[0:3] *= -1 # convert from P=1 to P=-1 if P > 0: ro[0:3] *= -1 # convert from P=1 to P=-1
if normalise: ro[0:3] /= np.linalg.norm(ro[0:3]) if normalise: ro[0:3] /= np.linalg.norm(ro[0:3])
if not np.isclose(np.linalg.norm(ro[0:3]), 1.0): if not np.isclose(np.linalg.norm(ro[0:3]), 1.0):
raise ValueError('Rodrigues rotation axis is not of unit length.\n{} {} {}'.format(*ro[0:3])) raise ValueError('Rodrigues rotation axis is not of unit length: {} {} {}.'.format(*ro[0:3]))
if ro[3] < 0.0: if ro[3] < 0.0:
raise ValueError('Rodrigues rotation angle not positive.\n{}'.format(ro[3])) raise ValueError('Rodrigues rotation angle not positive: {}.'.format(ro[3]))
return Rotation(Rotation.ro2qu(ro)) return Rotation(Rotation.ro2qu(ro))
@ -350,6 +350,9 @@ class Rotation:
else np.array(homochoric,dtype=float) else np.array(homochoric,dtype=float)
if P > 0: ho *= -1 # convert from P=1 to P=-1 if P > 0: ho *= -1 # convert from P=1 to P=-1
if np.linalg.norm(ho) > (3.*np.pi/4.)**(1./3.)+1e-9:
raise ValueError('Coordinate outside of the sphere: {} {} {}.'.format(ho))
return Rotation(Rotation.ho2qu(ho)) return Rotation(Rotation.ho2qu(ho))
@staticmethod @staticmethod
@ -358,6 +361,10 @@ class Rotation:
cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \ cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \
else np.array(cubochoric,dtype=float) else np.array(cubochoric,dtype=float)
if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9:
raise ValueError('Coordinate outside of the cube: {} {} {}.'.format(*cu))
ho = Rotation.cu2ho(cu) ho = Rotation.cu2ho(cu)
if P > 0: ho *= -1 # convert from P=1 to P=-1 if P > 0: ho *= -1 # convert from P=1 to P=-1
@ -383,7 +390,7 @@ class Rotation:
""" """
if not all(isinstance(item, Rotation) for item in rotations): if not all(isinstance(item, Rotation) for item in rotations):
raise TypeError("Only instances of Rotation can be averaged.") raise TypeError('Only instances of Rotation can be averaged.')
N = len(rotations) N = len(rotations)
if not weights: if not weights:
@ -441,34 +448,69 @@ class Rotation:
#---------- Quaternion ---------- #---------- Quaternion ----------
@staticmethod @staticmethod
def qu2om(qu): def qu2om(qu):
"""Quaternion to rotation matrix.""" if len(qu.shape) == 1:
qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) """Quaternion to rotation matrix."""
om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2)
om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2)
om[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) om[0,1] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3])
om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) om[1,0] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3])
om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) om[1,2] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1])
om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) om[2,1] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1])
om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) om[2,0] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2])
om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) om[0,2] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2])
return om if P > 0.0 else om.T else:
qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2)
om = np.block([qq + 2.0*qu[...,1:2]**2,
2.0*(qu[...,2:3]*qu[...,1:2]+qu[...,0:1]*qu[...,3:4]),
2.0*(qu[...,3:4]*qu[...,1:2]-qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,1:2]*qu[...,2:3]-qu[...,0:1]*qu[...,3:4]),
qq + 2.0*qu[...,2:3]**2,
2.0*(qu[...,3:4]*qu[...,2:3]+qu[...,0:1]*qu[...,1:2]),
2.0*(qu[...,1:2]*qu[...,3:4]+qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,2:3]*qu[...,3:4]-qu[...,0:1]*qu[...,1:2]),
qq + 2.0*qu[...,3:4]**2,
]).reshape(qu.shape[:-1]+(3,3))
return om if _P < 0.0 else np.swapaxes(om,(-1,-2))
@staticmethod @staticmethod
def qu2eu(qu): def qu2eu(qu):
"""Quaternion to Bunge-Euler angles.""" """Quaternion to Bunge-Euler angles."""
q03 = qu[0]**2+qu[3]**2 if len(qu.shape) == 1:
q12 = qu[1]**2+qu[2]**2 q03 = qu[0]**2+qu[3]**2
chi = np.sqrt(q03*q12) q12 = qu[1]**2+qu[2]**2
chi = np.sqrt(q03*q12)
if iszero(chi): if np.abs(q12) < 1.e-8:
eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) if iszero(q12) else \ eu = np.array([np.arctan2(-_P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0])
np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) elif np.abs(q03) < 1.e-8:
eu = np.array([np.arctan2( 2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0])
else:
eu = np.array([np.arctan2((-_P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]-qu[2]*qu[3])*chi ),
np.arctan2( 2.0*chi, q03-q12 ),
np.arctan2(( _P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-_P*qu[0]*qu[1]+qu[2]*qu[3])*chi )])
else: else:
eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), q02 = qu[...,0:1]*qu[...,2:3]
np.arctan2( 2.0*chi, q03-q12 ), q13 = qu[...,1:2]*qu[...,3:4]
np.arctan2(( P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]+qu[2]*qu[3])*chi )]) q01 = qu[...,0:1]*qu[...,1:2]
q23 = qu[...,2:3]*qu[...,3:4]
q03_s = qu[...,0:1]**2+qu[...,3:4]**2
q12_s = qu[...,1:2]**2+qu[...,2:3]**2
chi = np.sqrt(q03_s*q12_s)
# reduce Euler angles to definition range, i.e a lower limit of 0.0 eu = np.where(np.abs(q12_s) < 1.0e-8,
np.block([np.arctan2(-_P*2.0*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2),
np.zeros(qu.shape[:-1]+(2,))]),
np.where(np.abs(q03_s) < 1.0e-8,
np.block([np.arctan2( 2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2),
np.broadcast_to(np.pi,qu.shape[:-1]+(1,)),
np.zeros(qu.shape[:-1]+(1,))]),
np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi),
np.arctan2( 2.0*chi, q03_s-q12_s ),
np.arctan2(( _P*q02+q13)*chi, (-_P*q01+q23)*chi)])
)
)
# reduce Euler angles to definition range
eu[np.abs(eu)<1.e-6] = 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu return eu
@ -479,38 +521,65 @@ class Rotation:
Modified version of the original formulation, should be numerically more stable Modified version of the original formulation, should be numerically more stable
""" """
if iszero(qu[1]**2+qu[2]**2+qu[3]**2): # set axis to [001] if the angle is 0/360 if len(qu.shape) == 1:
ax = [ 0.0, 0.0, 1.0, 0.0 ] if np.abs(np.sum(qu[1:4]**2)) < 1.e-6: # set axis to [001] if the angle is 0/360
elif not iszero(qu[0]): ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) elif qu[0] > 1.e-6:
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2)
ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ] omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))
ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ])
else:
ax = ax = np.array([ qu[1], qu[2], qu[3], np.pi])
else: else:
ax = [ qu[1], qu[2], qu[3], np.pi] with np.errstate(invalid='ignore',divide='ignore'):
return np.array(ax) s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2)
omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0))
ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-6,qu.shape),
np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]),
np.block([qu[...,1:4]*s,omega]))
ax[np.sum(np.abs(qu[...,1:4])**2,axis=-1) < 1.0e-6,] = [0.0, 0.0, 1.0, 0.0]
return ax
@staticmethod @staticmethod
def qu2ro(qu): def qu2ro(qu):
"""Quaternion to Rodrigues-Frank vector.""" """Quaternion to Rodrigues-Frank vector."""
if iszero(qu[0]): if len(qu.shape) == 1:
ro = [qu[1], qu[2], qu[3], np.inf] if iszero(qu[0]):
ro = np.array([qu[1], qu[2], qu[3], np.inf])
else:
s = np.linalg.norm(qu[1:4])
ro = np.array([0.0,0.0,_P,0.0] if iszero(s) else \
[ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))])
else: else:
s = np.linalg.norm([qu[1],qu[2],qu[3]]) with np.errstate(invalid='ignore',divide='ignore'):
ro = [0.0,0.0,P,0.0] if iszero(s) else \ s = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True)
[ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))] ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape),
return np.array(ro) np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu.shape[:-1]+(1,))]),
np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s,
np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0)))
])
)
ro[np.abs(s).squeeze(-1) < 1.0e-12] = [0.0,0.0,_P,0.0]
return ro
@staticmethod @staticmethod
def qu2ho(qu): def qu2ho(qu):
"""Quaternion to homochoric vector.""" """Quaternion to homochoric vector."""
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) if len(qu.shape) == 1:
omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0))
if iszero(omega): if np.abs(omega) < 1.0e-12:
ho = np.array([ 0.0, 0.0, 0.0 ]) ho = np.zeros(3)
else:
ho = np.array([qu[1], qu[2], qu[3]])
f = 0.75 * ( omega - np.sin(omega) )
ho = ho/np.linalg.norm(ho) * f**(1./3.)
else: else:
ho = np.array([qu[1], qu[2], qu[3]]) with np.errstate(invalid='ignore'):
f = 0.75 * ( omega - np.sin(omega) ) omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0))
ho = ho/np.linalg.norm(ho) * f**(1./3.) ho = np.where(np.abs(omega) < 1.0e-12,
np.zeros(3),
qu[...,1:4]/np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True) \
* (0.75*(omega - np.sin(omega)))**(1./3.))
return ho return ho
@staticmethod @staticmethod
@ -532,37 +601,71 @@ class Rotation:
@staticmethod @staticmethod
def om2eu(om): def om2eu(om):
"""Rotation matrix to Bunge-Euler angles.""" """Rotation matrix to Bunge-Euler angles."""
if abs(om[2,2]) < 1.0: if len(om.shape) == 2:
zeta = 1.0/np.sqrt(1.0-om[2,2]**2) if not np.isclose(np.abs(om[2,2]),1.0,1.e-4):
eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
np.arccos(om[2,2]), eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) np.arccos(om[2,2]),
np.arctan2(om[0,2]*zeta, om[1,2]*zeta)])
else:
eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation
else: else:
eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation with np.errstate(invalid='ignore',divide='ignore'):
zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2)
# reduce Euler angles to definition range, i.e a lower limit of 0.0 eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-4),
np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]),
np.pi*0.5*(1-om[...,2,2:3]),
np.zeros(om.shape[:-2]+(1,)),
]),
np.block([np.arctan2(om[...,2,0:1]*zeta,-om[...,2,1:2]*zeta),
np.arccos(om[...,2,2:3]),
np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta)
])
)
eu[np.abs(eu)<1.e-6] = 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu return eu
@staticmethod @staticmethod
def om2ax(om): def om2ax(om):
"""Rotation matrix to axis angle pair.""" """Rotation matrix to axis angle pair."""
ax=np.empty(4) if len(om.shape) == 2:
ax=np.empty(4)
# first get the rotation angle # first get the rotation angle
t = 0.5*(om.trace() -1.0) t = 0.5*(om.trace() -1.0)
ax[3] = np.arccos(np.clip(t,-1.0,1.0)) ax[3] = np.arccos(np.clip(t,-1.0,1.0))
if np.abs(ax[3])<1.e-6:
if iszero(ax[3]): ax = np.array([ 0.0, 0.0, 1.0, 0.0])
ax = [ 0.0, 0.0, 1.0, 0.0] else:
w,vr = np.linalg.eig(om)
# next, find the eigenvalue (1,0j)
i = np.where(np.isclose(w,1.0+0.0j))[0][0]
ax[0:3] = np.real(vr[0:3,i])
diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]])
diagDelta[np.abs(diagDelta)<1.e-6] = 1.0
ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta))
else: else:
diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2],
om[...,2,0:1]-om[...,0,2:3],
om[...,0,1:2]-om[...,1,0:1]
])
diag_delta[np.abs(diag_delta)<1.e-6] = 1.0
t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,))
w,vr = np.linalg.eig(om) w,vr = np.linalg.eig(om)
# next, find the eigenvalue (1,0j) # mask duplicated real eigenvalues
i = np.where(np.isclose(w,1.0+0.0j))[0][0] w[np.isclose(w[...,0],1.0+0.0j),1:] = 0.
ax[0:3] = np.real(vr[0:3,i]) w[np.isclose(w[...,1],1.0+0.0j),2:] = 0.
diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) vr = np.swapaxes(vr,-1,-2)
ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) ax = np.where(np.abs(diag_delta)<0,
return np.array(ax) np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)),
np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \
*np.sign(diag_delta))
ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))])
ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0]
return ax
@staticmethod @staticmethod
def om2ro(om): def om2ro(om):
@ -584,57 +687,103 @@ class Rotation:
@staticmethod @staticmethod
def eu2qu(eu): def eu2qu(eu):
"""Bunge-Euler angles to quaternion.""" """Bunge-Euler angles to quaternion."""
ee = 0.5*eu if len(eu.shape) == 1:
cPhi = np.cos(ee[1]) ee = 0.5*eu
sPhi = np.sin(ee[1]) cPhi = np.cos(ee[1])
qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), sPhi = np.sin(ee[1])
-P*sPhi*np.cos(ee[0]-ee[2]), qu = np.array([ cPhi*np.cos(ee[0]+ee[2]),
-P*sPhi*np.sin(ee[0]-ee[2]), -_P*sPhi*np.cos(ee[0]-ee[2]),
-P*cPhi*np.sin(ee[0]+ee[2]) ]) -_P*sPhi*np.sin(ee[0]-ee[2]),
if qu[0] < 0.0: qu*=-1 -_P*cPhi*np.sin(ee[0]+ee[2]) ])
if qu[0] < 0.0: qu*=-1
else:
ee = 0.5*eu
cPhi = np.cos(ee[...,1:2])
sPhi = np.sin(ee[...,1:2])
qu = np.block([ cPhi*np.cos(ee[...,0:1]+ee[...,2:3]),
-_P*sPhi*np.cos(ee[...,0:1]-ee[...,2:3]),
-_P*sPhi*np.sin(ee[...,0:1]-ee[...,2:3]),
-_P*cPhi*np.sin(ee[...,0:1]+ee[...,2:3])])
qu[qu[...,0]<0.0]*=-1
return qu return qu
@staticmethod @staticmethod
def eu2om(eu): def eu2om(eu):
"""Bunge-Euler angles to rotation matrix.""" """Bunge-Euler angles to rotation matrix."""
c = np.cos(eu) if len(eu.shape) == 1:
s = np.sin(eu) c = np.cos(eu)
s = np.sin(eu)
om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]], om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]],
[-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]], [-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]],
[+s[0]*s[1], -c[0]*s[1], +c[1] ]]) [+s[0]*s[1], -c[0]*s[1], +c[1] ]])
else:
om[np.where(iszero(om))] = 0.0 c = np.cos(eu)
s = np.sin(eu)
om = np.block([+c[...,0:1]*c[...,2:3]-s[...,0:1]*s[...,2:3]*c[...,1:2],
+s[...,0:1]*c[...,2:3]+c[...,0:1]*s[...,2:3]*c[...,1:2],
+s[...,2:3]*s[...,1:2],
-c[...,0:1]*s[...,2:3]-s[...,0:1]*c[...,2:3]*c[...,1:2],
-s[...,0:1]*s[...,2:3]+c[...,0:1]*c[...,2:3]*c[...,1:2],
+c[...,2:3]*s[...,1:2],
+s[...,0:1]*s[...,1:2],
-c[...,0:1]*s[...,1:2],
+c[...,1:2]
]).reshape(eu.shape[:-1]+(3,3))
om[np.abs(om)<1.e-12] = 0.0
return om return om
@staticmethod @staticmethod
def eu2ax(eu): def eu2ax(eu):
"""Bunge-Euler angles to axis angle pair.""" """Bunge-Euler angles to axis angle pair."""
t = np.tan(eu[1]*0.5) if len(eu.shape) == 1:
sigma = 0.5*(eu[0]+eu[2]) t = np.tan(eu[1]*0.5)
delta = 0.5*(eu[0]-eu[2]) sigma = 0.5*(eu[0]+eu[2])
tau = np.linalg.norm([t,np.sin(sigma)]) delta = 0.5*(eu[0]-eu[2])
alpha = np.pi if iszero(np.cos(sigma)) else \ tau = np.linalg.norm([t,np.sin(sigma)])
2.0*np.arctan(tau/np.cos(sigma)) alpha = np.pi if iszero(np.cos(sigma)) else \
2.0*np.arctan(tau/np.cos(sigma))
if iszero(alpha): if np.abs(alpha)<1.e-6:
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
else:
ax = -_P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front
ax = np.append(ax,alpha)
if alpha < 0.0: ax *= -1.0 # ensure alpha is positive
else: else:
ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis angle pair so a minus sign in front t = np.tan(eu[...,1:2]*0.5)
ax = np.append(ax,alpha) sigma = 0.5*(eu[...,0:1]+eu[...,2:3])
if alpha < 0.0: ax *= -1.0 # ensure alpha is positive delta = 0.5*(eu[...,0:1]-eu[...,2:3])
tau = np.linalg.norm(np.block([t,np.sin(sigma)]),axis=-1,keepdims=True)
alpha = np.where(np.abs(np.cos(sigma))<1.e-12,np.pi,2.0*np.arctan(tau/np.cos(sigma)))
with np.errstate(invalid='ignore',divide='ignore'):
ax = np.where(np.broadcast_to(np.abs(alpha)<1.0e-12,eu.shape[:-1]+(4,)),
[0.0,0.0,1.0,0.0],
np.block([-_P/tau*t*np.cos(delta),
-_P/tau*t*np.sin(delta),
-_P/tau* np.sin(sigma),
alpha
]))
ax[(alpha<0.0).squeeze()] *=-1
return ax return ax
@staticmethod @staticmethod
def eu2ro(eu): def eu2ro(eu):
"""Bunge-Euler angles to Rodrigues-Frank vector.""" """Bunge-Euler angles to Rodrigues-Frank vector."""
ro = Rotation.eu2ax(eu) # convert to axis angle pair representation if len(eu.shape) == 1:
if ro[3] >= np.pi: # Differs from original implementation. check convention 5 ro = Rotation.eu2ax(eu) # convert to axis angle pair representation
ro[3] = np.inf if ro[3] >= np.pi: # Differs from original implementation. check convention 5
elif iszero(ro[3]): ro[3] = np.inf
ro = np.array([ 0.0, 0.0, P, 0.0 ]) elif iszero(ro[3]):
ro = np.array([ 0.0, 0.0, _P, 0.0 ])
else:
ro[3] = np.tan(ro[3]*0.5)
else: else:
ro[3] = np.tan(ro[3]*0.5) ax = Rotation.eu2ax(eu)
ro = np.block([ax[...,:3],np.tan(ax[...,3:4]*.5)])
ro[ax[...,3]>=np.pi,3] = np.inf
ro[np.abs(ax[...,3])<1.e-16] = [ 0.0, 0.0, _P, 0.0 ]
return ro return ro
@staticmethod @staticmethod
@ -652,27 +801,47 @@ class Rotation:
@staticmethod @staticmethod
def ax2qu(ax): def ax2qu(ax):
"""Axis angle pair to quaternion.""" """Axis angle pair to quaternion."""
if iszero(ax[3]): if len(ax.shape) == 1:
qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) if np.abs(ax[3])<1.e-6:
qu = np.array([ 1.0, 0.0, 0.0, 0.0 ])
else:
c = np.cos(ax[3]*0.5)
s = np.sin(ax[3]*0.5)
qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ])
return qu
else: else:
c = np.cos(ax[3]*0.5) c = np.cos(ax[...,3:4]*.5)
s = np.sin(ax[3]*0.5) s = np.sin(ax[...,3:4]*.5)
qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s]))
return qu return qu
@staticmethod @staticmethod
def ax2om(ax): def ax2om(ax):
"""Axis angle pair to rotation matrix.""" """Axis angle pair to rotation matrix."""
c = np.cos(ax[3]) if len(ax.shape) == 1:
s = np.sin(ax[3]) c = np.cos(ax[3])
omc = 1.0-c s = np.sin(ax[3])
om=np.diag(ax[0:3]**2*omc + c) omc = 1.0-c
om=np.diag(ax[0:3]**2*omc + c)
for idx in [[0,1,2],[1,2,0],[2,0,1]]: for idx in [[0,1,2],[1,2,0],[2,0,1]]:
q = omc*ax[idx[0]] * ax[idx[1]] q = omc*ax[idx[0]] * ax[idx[1]]
om[idx[0],idx[1]] = q + s*ax[idx[2]] om[idx[0],idx[1]] = q + s*ax[idx[2]]
om[idx[1],idx[0]] = q - s*ax[idx[2]] om[idx[1],idx[0]] = q - s*ax[idx[2]]
return om if P < 0.0 else om.T else:
c = np.cos(ax[...,3:4])
s = np.sin(ax[...,3:4])
omc = 1. -c
om = np.block([c+omc*ax[...,0:1]**2,
omc*ax[...,0:1]*ax[...,1:2] + s*ax[...,2:3],
omc*ax[...,0:1]*ax[...,2:3] - s*ax[...,1:2],
omc*ax[...,0:1]*ax[...,1:2] - s*ax[...,2:3],
c+omc*ax[...,1:2]**2,
omc*ax[...,1:2]*ax[...,2:3] + s*ax[...,0:1],
omc*ax[...,0:1]*ax[...,2:3] + s*ax[...,1:2],
omc*ax[...,1:2]*ax[...,2:3] - s*ax[...,0:1],
c+omc*ax[...,2:3]**2]).reshape(ax.shape[:-1]+(3,3))
return om if _P < 0.0 else np.swapaxes(om,(-1,-2))
@staticmethod @staticmethod
def ax2eu(ax): def ax2eu(ax):
@ -682,21 +851,35 @@ class Rotation:
@staticmethod @staticmethod
def ax2ro(ax): def ax2ro(ax):
"""Axis angle pair to Rodrigues-Frank vector.""" """Axis angle pair to Rodrigues-Frank vector."""
if iszero(ax[3]): if len(ax.shape) == 1:
ro = [ 0.0, 0.0, P, 0.0 ] if np.abs(ax[3])<1.e-6:
ro = [ 0.0, 0.0, _P, 0.0 ]
else:
ro = [ax[0], ax[1], ax[2]]
# 180 degree case
ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \
[np.tan(ax[3]*0.5)]
return np.array(ro)
else: else:
ro = [ax[0], ax[1], ax[2]] ro = np.block([ax[...,:3],
# 180 degree case np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0),
ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ np.inf,
[np.tan(ax[3]*0.5)] np.tan(ax[...,3:4]*0.5))
return np.array(ro) ])
ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0]
return ro
@staticmethod @staticmethod
def ax2ho(ax): def ax2ho(ax):
"""Axis angle pair to homochoric vector.""" """Axis angle pair to homochoric vector."""
f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) if len(ax.shape) == 1:
ho = ax[0:3] * f f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0)
return ho ho = ax[0:3] * f
return ho
else:
f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0)
ho = ax[...,:3] * f
return ho
@staticmethod @staticmethod
def ax2cu(ax): def ax2cu(ax):
@ -723,27 +906,38 @@ class Rotation:
@staticmethod @staticmethod
def ro2ax(ro): def ro2ax(ro):
"""Rodrigues-Frank vector to axis angle pair.""" """Rodrigues-Frank vector to axis angle pair."""
ta = ro[3] if len(ro.shape) == 1:
if np.abs(ro[3]) < 1.e-6:
if iszero(ta): ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
ax = [ 0.0, 0.0, 1.0, 0.0 ] elif not np.isfinite(ro[3]):
elif not np.isfinite(ta): ax = np.array([ ro[0], ro[1], ro[2], np.pi ])
ax = [ ro[0], ro[1], ro[2], np.pi ] else:
angle = 2.0*np.arctan(ro[3])
ta = np.linalg.norm(ro[0:3])
ax = np.array([ ro[0]*ta, ro[1]*ta, ro[2]*ta, angle ])
else: else:
angle = 2.0*np.arctan(ta) with np.errstate(invalid='ignore',divide='ignore'):
ta = 1.0/np.linalg.norm(ro[0:3]) ax = np.where(np.isfinite(ro[...,3:4]),
ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ] np.block([ro[...,0:3]*np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True),2.*np.arctan(ro[...,3:4])]),
return np.array(ax) np.block([ro[...,0:3],np.broadcast_to(np.pi,ro[...,3:4].shape)]))
ax[np.abs(ro[...,3]) < 1.e-6] = np.array([ 0.0, 0.0, 1.0, 0.0 ])
return ax
@staticmethod @staticmethod
def ro2ho(ro): def ro2ho(ro):
"""Rodrigues-Frank vector to homochoric vector.""" """Rodrigues-Frank vector to homochoric vector."""
if iszero(np.sum(ro[0:3]**2.0)): if len(ro.shape) == 1:
ho = [ 0.0, 0.0, 0.0 ] if np.sum(ro[0:3]**2.0) < 1.e-6:
ho = np.zeros(3)
else:
f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi
ho = ro[0:3] * (0.75*f)**(1.0/3.0)
else: else:
f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi)
ho = ro[0:3] * (0.75*f)**(1.0/3.0) ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape),
return np.array(ho) np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0))
return ho
@staticmethod @staticmethod
def ro2cu(ro): def ro2cu(ro):
@ -778,19 +972,31 @@ class Rotation:
+0.0001703481934140054, -0.00012062065004116828, +0.0001703481934140054, -0.00012062065004116828,
+0.000059719705868660826, -0.00001980756723965647, +0.000059719705868660826, -0.00001980756723965647,
+0.000003953714684212874, -0.00000036555001439719544]) +0.000003953714684212874, -0.00000036555001439719544])
# normalize h and store the magnitude if len(ho.shape) == 1:
hmag_squared = np.sum(ho**2.) # normalize h and store the magnitude
if iszero(hmag_squared): hmag_squared = np.sum(ho**2.)
ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) if iszero(hmag_squared):
else: ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
hm = hmag_squared else:
hm = hmag_squared
# convert the magnitude to the rotation angle # convert the magnitude to the rotation angle
s = tfit[0] + tfit[1] * hmag_squared
for i in range(2,16):
hm *= hmag_squared
s += tfit[i] * hm
ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0)))
else:
hmag_squared = np.sum(ho**2.,axis=-1,keepdims=True)
hm = hmag_squared.copy()
s = tfit[0] + tfit[1] * hmag_squared s = tfit[0] + tfit[1] * hmag_squared
for i in range(2,16): for i in range(2,16):
hm *= hmag_squared hm *= hmag_squared
s += tfit[i] * hm s += tfit[i] * hm
ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) with np.errstate(invalid='ignore'):
ax = np.where(np.broadcast_to(np.abs(hmag_squared)<1.e-6,ho.shape[:-1]+(4,)),
[ 0.0, 0.0, 1.0, 0.0 ],
np.block([ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))]))
return ax return ax
@staticmethod @staticmethod
@ -801,7 +1007,10 @@ class Rotation:
@staticmethod @staticmethod
def ho2cu(ho): def ho2cu(ho):
"""Homochoric vector to cubochoric vector.""" """Homochoric vector to cubochoric vector."""
return ball_to_cube(ho) if len(ho.shape) == 1:
return ball_to_cube(ho)
else:
raise NotImplementedError
#---------- Cubochoric ---------- #---------- Cubochoric ----------
@ -833,4 +1042,7 @@ class Rotation:
@staticmethod @staticmethod
def cu2ho(cu): def cu2ho(cu):
"""Cubochoric vector to homochoric vector.""" """Cubochoric vector to homochoric vector."""
return cube_to_ball(cu) if len(cu.shape) == 1:
return cube_to_ball(cu)
else:
raise NotImplementedError

222
python/tests/test_Rotation.py
View File

@ -4,13 +4,83 @@ import pytest
import numpy as np import numpy as np
from damask import Rotation from damask import Rotation
n = 1000 n = 1100
atol=1.e-4
scatter=1.e-2
@pytest.fixture @pytest.fixture
def default(): def default():
"""A set of n random rotations.""" """A set of n random rotations."""
return [Rotation.fromRandom() for r in range(n)] specials = np.array(
[np.array([ 1.0, 0.0, 0.0, 0.0]),
#-----------------------------------------------
np.array([0.0, 1.0, 0.0, 0.0]),
np.array([0.0, 0.0, 1.0, 0.0]),
np.array([0.0, 0.0, 0.0, 1.0]),
np.array([0.0,-1.0, 0.0, 0.0]),
np.array([0.0, 0.0,-1.0, 0.0]),
np.array([0.0, 0.0, 0.0,-1.0]),
#-----------------------------------------------
np.array([1.0, 1.0, 0.0, 0.0])/np.sqrt(2.),
np.array([1.0, 0.0, 1.0, 0.0])/np.sqrt(2.),
np.array([1.0, 0.0, 0.0, 1.0])/np.sqrt(2.),
np.array([0.0, 1.0, 1.0, 0.0])/np.sqrt(2.),
np.array([0.0, 1.0, 0.0, 1.0])/np.sqrt(2.),
np.array([0.0, 0.0, 1.0, 1.0])/np.sqrt(2.),
#-----------------------------------------------
np.array([1.0,-1.0, 0.0, 0.0])/np.sqrt(2.),
np.array([1.0, 0.0,-1.0, 0.0])/np.sqrt(2.),
np.array([1.0, 0.0, 0.0,-1.0])/np.sqrt(2.),
np.array([0.0, 1.0,-1.0, 0.0])/np.sqrt(2.),
np.array([0.0, 1.0, 0.0,-1.0])/np.sqrt(2.),
np.array([0.0, 0.0, 1.0,-1.0])/np.sqrt(2.),
#-----------------------------------------------
np.array([0.0, 1.0,-1.0, 0.0])/np.sqrt(2.),
np.array([0.0, 1.0, 0.0,-1.0])/np.sqrt(2.),
np.array([0.0, 0.0, 1.0,-1.0])/np.sqrt(2.),
#-----------------------------------------------
np.array([0.0,-1.0,-1.0, 0.0])/np.sqrt(2.),
np.array([0.0,-1.0, 0.0,-1.0])/np.sqrt(2.),
np.array([0.0, 0.0,-1.0,-1.0])/np.sqrt(2.),
#-----------------------------------------------
np.array([1.0, 1.0, 1.0, 0.0])/np.sqrt(3.),
np.array([1.0, 1.0, 0.0, 1.0])/np.sqrt(3.),
np.array([1.0, 0.0, 1.0, 1.0])/np.sqrt(3.),
np.array([1.0,-1.0, 1.0, 0.0])/np.sqrt(3.),
np.array([1.0,-1.0, 0.0, 1.0])/np.sqrt(3.),
np.array([1.0, 0.0,-1.0, 1.0])/np.sqrt(3.),
np.array([1.0, 1.0,-1.0, 0.0])/np.sqrt(3.),
np.array([1.0, 1.0, 0.0,-1.0])/np.sqrt(3.),
np.array([1.0, 0.0, 1.0,-1.0])/np.sqrt(3.),
np.array([1.0,-1.0,-1.0, 0.0])/np.sqrt(3.),
np.array([1.0,-1.0, 0.0,-1.0])/np.sqrt(3.),
np.array([1.0, 0.0,-1.0,-1.0])/np.sqrt(3.),
#-----------------------------------------------
np.array([0.0, 1.0, 1.0, 1.0])/np.sqrt(3.),
np.array([0.0, 1.0,-1.0, 1.0])/np.sqrt(3.),
np.array([0.0, 1.0, 1.0,-1.0])/np.sqrt(3.),
np.array([0.0,-1.0, 1.0, 1.0])/np.sqrt(3.),
np.array([0.0,-1.0,-1.0, 1.0])/np.sqrt(3.),
np.array([0.0,-1.0, 1.0,-1.0])/np.sqrt(3.),
np.array([0.0,-1.0,-1.0,-1.0])/np.sqrt(3.),
#-----------------------------------------------
np.array([1.0, 1.0, 1.0, 1.0])/2.,
np.array([1.0,-1.0, 1.0, 1.0])/2.,
np.array([1.0, 1.0,-1.0, 1.0])/2.,
np.array([1.0, 1.0, 1.0,-1.0])/2.,
np.array([1.0,-1.0,-1.0, 1.0])/2.,
np.array([1.0,-1.0, 1.0,-1.0])/2.,
np.array([1.0, 1.0,-1.0,-1.0])/2.,
np.array([1.0,-1.0,-1.0,-1.0])/2.,
])
specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape)
specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1)
specials_scatter[specials_scatter[:,0]<0]*=-1
return [Rotation.fromQuaternion(s) for s in specials] + \
[Rotation.fromQuaternion(s) for s in specials_scatter] + \
[Rotation.fromRandom() for _ in range(n-len(specials)-len(specials_scatter))]
@pytest.fixture @pytest.fixture
def reference_dir(reference_dir_base): def reference_dir(reference_dir_base):
@ -22,35 +92,151 @@ class TestRotation:
def test_Eulers(self,default): def test_Eulers(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asQuaternion(), m = rot.asQuaternion()
Rotation.fromEulers(rot.asEulers()).asQuaternion()) o = Rotation.fromEulers(rot.asEulers()).asQuaternion()
ok = np.allclose(m,o,atol=atol)
if np.isclose(rot.asQuaternion()[0],0.0,atol=atol):
ok = ok or np.allclose(m*-1.,o,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and np.isclose(np.linalg.norm(o),1.0)
def test_AxisAngle(self,default): def test_AxisAngle(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asEulers(), m = rot.asEulers()
Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers()) o = Rotation.fromAxisAngle(rot.asAxisAngle()).asEulers()
u = np.array([np.pi*2,np.pi,np.pi*2])
ok = np.allclose(m,o,atol=atol)
ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol)
if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol):
sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]])
ok = ok or np.isclose(sum_phi[0],sum_phi[1],atol=atol)
print(m,o,rot.asQuaternion())
assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all()
def test_Matrix(self,default): def test_Matrix(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asAxisAngle(), m = rot.asAxisAngle()
Rotation.fromMatrix(rot.asMatrix()).asAxisAngle()) o = Rotation.fromAxisAngle(rot.asAxisAngle()).asAxisAngle()
ok = np.allclose(m,o,atol=atol)
if np.isclose(m[3],np.pi,atol=atol):
ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9
def test_Rodriques(self,default): def test_Rodrigues(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asMatrix(), m = rot.asMatrix()
Rotation.fromRodrigues(rot.asRodrigues()).asMatrix()) o = Rotation.fromRodrigues(rot.asRodrigues()).asMatrix()
ok = np.allclose(m,o,atol=atol)
print(m,o)
assert ok and np.isclose(np.linalg.det(o),1.0)
def test_Homochoric(self,default): def test_Homochoric(self,default):
cutoff = np.tan(np.pi*.5*(1.-1e-4))
for rot in default: for rot in default:
assert np.allclose(rot.asRodrigues(), m = rot.asRodrigues()
Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues(),rtol=1.e-4) o = Rotation.fromHomochoric(rot.asHomochoric()).asRodrigues()
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
ok = ok or np.isclose(m[3],0.0,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0)
def test_Cubochoric(self,default): def test_Cubochoric(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asHomochoric(), m = rot.asHomochoric()
Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric()) o = Rotation.fromCubochoric(rot.asCubochoric()).asHomochoric()
ok = np.allclose(m,o,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and np.linalg.norm(o) < (3.*np.pi/4.)**(1./3.) + 1.e-9
def test_Quaternion(self,default): def test_Quaternion(self,default):
for rot in default: for rot in default:
assert np.allclose(rot.asCubochoric(), m = rot.asCubochoric()
Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric()) o = Rotation.fromQuaternion(rot.asQuaternion()).asCubochoric()
ok = np.allclose(m,o,atol=atol)
print(m,o,rot.asQuaternion())
assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
@pytest.mark.parametrize('conversion',[Rotation.qu2om,
Rotation.qu2eu,
Rotation.qu2ax,
Rotation.qu2ro,
Rotation.qu2ho])
def test_quaternion_vectorization(self,default,conversion):
qu = np.array([rot.asQuaternion() for rot in default])
conversion(qu.reshape(qu.shape[0]//2,-1,4))
co = conversion(qu)
for q,c in zip(qu,co):
print(q,c)
assert np.allclose(conversion(q),c)
@pytest.mark.parametrize('conversion',[Rotation.om2qu,
Rotation.om2eu,
Rotation.om2ax,
Rotation.om2ro,
Rotation.om2ho,
])
def test_matrix_vectorization(self,default,conversion):
om = np.array([rot.asMatrix() for rot in default])
conversion(om.reshape(om.shape[0]//2,-1,3,3))
co = conversion(om)
for o,c in zip(om,co):
print(o,c)
assert np.allclose(conversion(o),c)
@pytest.mark.parametrize('conversion',[Rotation.eu2qu,
Rotation.eu2om,
Rotation.eu2ax,
Rotation.eu2ro,
Rotation.eu2ho,
])
def test_Euler_vectorization(self,default,conversion):
eu = np.array([rot.asEulers() for rot in default])
conversion(eu.reshape(eu.shape[0]//2,-1,3))
co = conversion(eu)
for e,c in zip(eu,co):
print(e,c)
assert np.allclose(conversion(e),c)
@pytest.mark.parametrize('conversion',[Rotation.ax2qu,
Rotation.ax2om,
Rotation.ax2eu,
Rotation.ax2ro,
Rotation.ax2ho,
])
def test_axisAngle_vectorization(self,default,conversion):
ax = np.array([rot.asAxisAngle() for rot in default])
conversion(ax.reshape(ax.shape[0]//2,-1,4))
co = conversion(ax)
for a,c in zip(ax,co):
print(a,c)
assert np.allclose(conversion(a),c)
@pytest.mark.parametrize('conversion',[Rotation.ro2qu,
Rotation.ro2om,
Rotation.ro2eu,
Rotation.ro2ax,
Rotation.ro2ho,
])
def test_Rodrigues_vectorization(self,default,conversion):
ro = np.array([rot.asRodrigues() for rot in default])
conversion(ro.reshape(ro.shape[0]//2,-1,4))
co = conversion(ro)
for r,c in zip(ro,co):
print(r,c)
assert np.allclose(conversion(r),c)
@pytest.mark.parametrize('conversion',[Rotation.ho2qu,
Rotation.ho2om,
Rotation.ho2eu,
Rotation.ho2ax,
Rotation.ho2ro,
])
def test_homochoric_vectorization(self,default,conversion):
ho = np.array([rot.asHomochoric() for rot in default])
conversion(ho.reshape(ho.shape[0]//2,-1,3))
co = conversion(ho)
for h,c in zip(ho,co):
print(h,c)
assert np.allclose(conversion(h),c)

42
src/Lambert.f90
View File

@ -70,13 +70,13 @@ contains
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
pure function Lambert_CubeToBall(cube) result(ball) pure function Lambert_CubeToBall(cube) result(ball)
real(pReal), intent(in), dimension(3) :: cube real(pReal), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ real(pReal), dimension(3) :: ball, LamXYZ, XYZ
real(pReal), dimension(2) :: T real(pReal), dimension(2) :: T
real(pReal) :: c, s, q real(pReal) :: c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal real(pReal), parameter :: eps = 1.0e-8_pReal
integer, dimension(3) :: p integer, dimension(3,2) :: p
integer, dimension(2) :: order integer, dimension(2) :: order
if (maxval(abs(cube)) > AP/2.0+eps) then if (maxval(abs(cube)) > AP/2.0+eps) then
ball = IEEE_value(cube,IEEE_positive_inf) ball = IEEE_value(cube,IEEE_positive_inf)
@ -89,7 +89,7 @@ pure function Lambert_CubeToBall(cube) result(ball)
else center else center
! get pyramide and scale by grid parameter ratio ! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cube) p = GetPyramidOrder(cube)
XYZ = cube(p) * sc XYZ = cube(p(:,1)) * sc
! intercept all the points along the z-axis ! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then special: if (all(dEq0(XYZ(1:2)))) then
@ -112,7 +112,7 @@ pure function Lambert_CubeToBall(cube) result(ball)
endif special endif special
! reverse the coordinates back to order according to the original pyramid number ! reverse the coordinates back to order according to the original pyramid number
ball = LamXYZ(p) ball = LamXYZ(p(:,2))
endif center endif center
@ -126,11 +126,11 @@ end function Lambert_CubeToBall
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
pure function Lambert_BallToCube(xyz) result(cube) pure function Lambert_BallToCube(xyz) result(cube)
real(pReal), intent(in), dimension(3) :: xyz real(pReal), intent(in), dimension(3) :: xyz
real(pReal), dimension(3) :: cube, xyz1, xyz3 real(pReal), dimension(3) :: cube, xyz1, xyz3
real(pReal), dimension(2) :: Tinv, xyz2 real(pReal), dimension(2) :: Tinv, xyz2
real(pReal) :: rs, qxy, q2, sq2, q, tt real(pReal) :: rs, qxy, q2, sq2, q, tt
integer, dimension(3) :: p integer, dimension(3,2) :: p
rs = norm2(xyz) rs = norm2(xyz)
if (rs > R1) then if (rs > R1) then
@ -142,7 +142,7 @@ pure function Lambert_BallToCube(xyz) result(cube)
cube = 0.0_pReal cube = 0.0_pReal
else center else center
p = GetPyramidOrder(xyz) p = GetPyramidOrder(xyz)
xyz3 = xyz(p) xyz3 = xyz(p(:,1))
! inverse M_3 ! inverse M_3
xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) ) xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) )
@ -166,7 +166,7 @@ pure function Lambert_BallToCube(xyz) result(cube)
xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / pref ] /sc
! reverse the coordinates back to order according to the original pyramid number ! reverse the coordinates back to order according to the original pyramid number
cube = xyz1(p) cube = xyz1(p(:,2))
endif center endif center
@ -180,18 +180,18 @@ end function Lambert_BallToCube
!-------------------------------------------------------------------------- !--------------------------------------------------------------------------
pure function GetPyramidOrder(xyz) pure function GetPyramidOrder(xyz)
real(pReal),intent(in),dimension(3) :: xyz real(pReal),intent(in),dimension(3) :: xyz
integer, dimension(3) :: GetPyramidOrder integer, dimension(3,2) :: GetPyramidOrder
if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. & if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. &
((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then
GetPyramidOrder = [1,2,3] GetPyramidOrder = reshape([[1,2,3],[1,2,3]],[3,2])
else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. & else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. &
((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then
GetPyramidOrder = [2,3,1] GetPyramidOrder = reshape([[2,3,1],[3,1,2]],[3,2])
else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. & else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. &
((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then
GetPyramidOrder = [3,1,2] GetPyramidOrder = reshape([[3,1,2],[2,3,1]],[3,2])
else else
GetPyramidOrder = -1 ! should be impossible, but might simplify debugging GetPyramidOrder = -1 ! should be impossible, but might simplify debugging
end if end if

341
src/homogenization.f90
View File

@ -27,33 +27,22 @@ module homogenization
implicit none implicit none
private private
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
logical, public :: & logical, public :: &
terminallyIll = .false. !< at least one material point is terminally ill terminallyIll = .false. !< at least one material point is terminally ill
real(pReal), dimension(:,:,:,:), allocatable, public :: &
materialpoint_F0, & !< def grad of IP at start of FE increment
materialpoint_F, & !< def grad of IP to be reached at end of FE increment
materialpoint_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public :: &
materialpoint_dPdF !< tangent of first P--K stress at IP
real(pReal), dimension(:,:,:,:), allocatable :: & !--------------------------------------------------------------------------------------------------
materialpoint_subF0, & !< def grad of IP at beginning of homogenization increment ! General variables for the homogenization at a material point
materialpoint_subF !< def grad of IP to be reached at end of homog inc real(pReal), dimension(:,:,:,:), allocatable, public :: &
real(pReal), dimension(:,:), allocatable :: & materialpoint_F0, & !< def grad of IP at start of FE increment
materialpoint_subFrac, & materialpoint_F !< def grad of IP to be reached at end of FE increment
materialpoint_subStep, & real(pReal), dimension(:,:,:,:), allocatable, public, protected :: &
materialpoint_subdt materialpoint_P !< first P--K stress of IP
logical, dimension(:,:), allocatable :: & real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: &
materialpoint_requested, & materialpoint_dPdF !< tangent of first P--K stress at IP
materialpoint_converged
logical, dimension(:,:,:), allocatable :: &
materialpoint_doneAndHappy
type :: tNumerics type :: tNumerics
integer :: & integer :: &
nMPstate !< materialpoint state loop limit nMPstate !< materialpoint state loop limit
real(pReal) :: & real(pReal) :: &
subStepMinHomog, & !< minimum (relative) size of sub-step allowed during cutback in homogenization subStepMinHomog, & !< minimum (relative) size of sub-step allowed during cutback in homogenization
subStepSizeHomog, & !< size of first substep when cutback in homogenization subStepSizeHomog, & !< size of first substep when cutback in homogenization
@ -161,15 +150,7 @@ subroutine homogenization_init
allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity
materialpoint_F = materialpoint_F0 ! initialize to identity materialpoint_F = materialpoint_F0 ! initialize to identity
allocate(materialpoint_subF0(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_subF(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_subFrac(discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_subStep(discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_subdt(discretization_nIP,discretization_nElem), source=0.0_pReal)
allocate(materialpoint_requested(discretization_nIP,discretization_nElem), source=.false.)
allocate(materialpoint_converged(discretization_nIP,discretization_nElem), source=.true.)
allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem), source=.true.)
write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6)
@ -203,6 +184,16 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
e, & !< element number e, & !< element number
mySource, & mySource, &
myNgrains myNgrains
real(pReal), dimension(3,3) :: &
subF
real(pReal), dimension(discretization_nIP,discretization_nElem) :: &
subFrac, &
subStep
logical, dimension(discretization_nIP,discretization_nElem) :: &
requested, &
converged
logical, dimension(2,discretization_nIP,discretization_nElem) :: &
doneAndHappy
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then
@ -216,7 +207,7 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! initialize restoration points of ... ! initialize restoration points
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
do i = FEsolving_execIP(1),FEsolving_execIP(2); do i = FEsolving_execIP(1),FEsolving_execIP(2);
@ -238,74 +229,60 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
enddo enddo
subFrac(i,e) = 0.0_pReal
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_F0(1:3,1:3,i,e) converged(i,e) = .false. ! pretend failed step ...
materialpoint_subFrac(i,e) = 0.0_pReal subStep(i,e) = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
materialpoint_subStep(i,e) = 1.0_pReal/num%subStepSizeHomog ! <<added to adopt flexibility in cutback size>> requested(i,e) = .true. ! everybody requires calculation
materialpoint_converged(i,e) = .false. ! pretend failed step of twice the required size
materialpoint_requested(i,e) = .true. ! everybody requires calculation
if (homogState(material_homogenizationAt(e))%sizeState > 0) & if (homogState(material_homogenizationAt(e))%sizeState > 0) &
homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & homogState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal homogenization state homogState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e))
if (thermalState(material_homogenizationAt(e))%sizeState > 0) & if (thermalState(material_homogenizationAt(e))%sizeState > 0) &
thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & thermalState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
thermalState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal thermal state thermalState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e))
if (damageState(material_homogenizationAt(e))%sizeState > 0) & if (damageState(material_homogenizationAt(e))%sizeState > 0) &
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e)) ! ...internal damage state damageState(material_homogenizationAt(e))%State0( :,material_homogenizationMemberAt(i,e))
enddo enddo
enddo enddo
NiterationHomog = 0 NiterationHomog = 0
cutBackLooping: do while (.not. terminallyIll .and. & cutBackLooping: do while (.not. terminallyIll .and. &
any(materialpoint_subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > num%subStepMinHomog)) any(subStep(:,FEsolving_execELem(1):FEsolving_execElem(2)) > num%subStepMinHomog))
!$OMP PARALLEL DO PRIVATE(myNgrains) !$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
converged: if (materialpoint_converged(i,e)) then if (converged(i,e)) then
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 & if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 &
.and. ((e == debug_e .and. i == debug_i) & .and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0)) then .or. .not. iand(debug_level(debug_homogenization),debug_levelSelective) /= 0)) then
write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', & write(6,'(a,1x,f12.8,1x,a,1x,f12.8,1x,a,i8,1x,i2/)') '<< HOMOG >> winding forward from', &
materialpoint_subFrac(i,e), 'to current materialpoint_subFrac', & subFrac(i,e), 'to current subFrac', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),'in materialpoint_stressAndItsTangent at el ip',e,i subFrac(i,e)+subStep(i,e),'in materialpoint_stressAndItsTangent at el ip',e,i
endif endif
#endif #endif
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! calculate new subStep and new subFrac ! calculate new subStep and new subFrac
materialpoint_subFrac(i,e) = materialpoint_subFrac(i,e) + materialpoint_subStep(i,e) subFrac(i,e) = subFrac(i,e) + subStep(i,e)
materialpoint_subStep(i,e) = min(1.0_pReal-materialpoint_subFrac(i,e), & subStep(i,e) = min(1.0_pReal-subFrac(i,e),num%stepIncreaseHomog*subStep(i,e)) ! introduce flexibility for step increase/acceleration
num%stepIncreaseHomog*materialpoint_subStep(i,e)) ! introduce flexibility for step increase/acceleration
steppingNeeded: if (materialpoint_subStep(i,e) > num%subStepMinHomog) then steppingNeeded: if (subStep(i,e) > num%subStepMinHomog) then
! wind forward grain starting point of... ! wind forward grain starting point
crystallite_partionedF0 (1:3,1:3,1:myNgrains,i,e) = & crystallite_partionedF0 (1:3,1:3,1:myNgrains,i,e) = crystallite_partionedF(1:3,1:3,1:myNgrains,i,e)
crystallite_partionedF(1:3,1:3,1:myNgrains,i,e) crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fp (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) = crystallite_Lp (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedFp0 (1:3,1:3,1:myNgrains,i,e) = & crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e) = crystallite_Fi (1:3,1:3,1:myNgrains,i,e)
crystallite_Fp (1:3,1:3,1:myNgrains,i,e) crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e) = crystallite_Li (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedS0 (1:3,1:3,1:myNgrains,i,e) = crystallite_S (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedLp0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Lp (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedFi0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Fi (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedLi0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_Li (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedS0 (1:3,1:3,1:myNgrains,i,e) = &
crystallite_S (1:3,1:3,1:myNgrains,i,e)
do g = 1,myNgrains do g = 1,myNgrains
plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) = & plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) = &
@ -326,15 +303,12 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e)) damageState(material_homogenizationAt(e))%State (:,material_homogenizationMemberAt(i,e))
materialpoint_subF0(1:3,1:3,i,e) = materialpoint_subF(1:3,1:3,i,e)
endif steppingNeeded endif steppingNeeded
else converged else
if ( (myNgrains == 1 .and. materialpoint_subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite if ( (myNgrains == 1 .and. subStep(i,e) <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
num%subStepSizeHomog * materialpoint_subStep(i,e) <= num%subStepMinHomog ) then ! would require too small subStep num%subStepSizeHomog * subStep(i,e) <= num%subStepMinHomog ) then ! would require too small subStep
! cutback makes no sense ! cutback makes no sense
!$OMP FLUSH(terminallyIll)
if (.not. terminallyIll) then ! so first signals terminally ill... if (.not. terminallyIll) then ! so first signals terminally ill...
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill' write(6,*) 'Integration point ', i,' at element ', e, ' terminally ill'
@ -342,32 +316,27 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
endif endif
terminallyIll = .true. ! ...and kills all others terminallyIll = .true. ! ...and kills all others
else ! cutback makes sense else ! cutback makes sense
materialpoint_subStep(i,e) = num%subStepSizeHomog * materialpoint_subStep(i,e) ! crystallite had severe trouble, so do a significant cutback subStep(i,e) = num%subStepSizeHomog * subStep(i,e) ! crystallite had severe trouble, so do a significant cutback
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 & if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0 &
.and. ((e == debug_e .and. i == debug_i) & .and. ((e == debug_e .and. i == debug_i) &
.or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0)) then .or. .not. iand(debug_level(debug_homogenization), debug_levelSelective) /= 0)) then
write(6,'(a,1x,f12.8,a,i8,1x,i2/)') & write(6,'(a,1x,f12.8,a,i8,1x,i2/)') &
'<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep:',& '<< HOMOG >> cutback step in materialpoint_stressAndItsTangent with new subStep:',&
materialpoint_subStep(i,e),' at el ip',e,i subStep(i,e),' at el ip',e,i
endif endif
#endif #endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! restore... ! restore
if (materialpoint_subStep(i,e) < 1.0_pReal) then ! protect against fake cutback from \Delta t = 2 to 1. Maybe that "trick" is not necessary anymore at all? I.e. start with \Delta t = 1 if (subStep(i,e) < 1.0_pReal) then ! protect against fake cutback from \Delta t = 2 to 1. Maybe that "trick" is not necessary anymore at all? I.e. start with \Delta t = 1
crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = & crystallite_Lp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e)
crystallite_partionedLp0(1:3,1:3,1:myNgrains,i,e) crystallite_Li(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e)
crystallite_Li(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedLi0(1:3,1:3,1:myNgrains,i,e)
endif ! maybe protecting everything from overwriting (not only L) makes even more sense endif ! maybe protecting everything from overwriting (not only L) makes even more sense
crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = & crystallite_Fp(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e)
crystallite_partionedFp0(1:3,1:3,1:myNgrains,i,e) crystallite_Fi(1:3,1:3,1:myNgrains,i,e) = crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e)
crystallite_Fi(1:3,1:3,1:myNgrains,i,e) = & crystallite_S (1:3,1:3,1:myNgrains,i,e) = crystallite_partionedS0 (1:3,1:3,1:myNgrains,i,e)
crystallite_partionedFi0(1:3,1:3,1:myNgrains,i,e)
crystallite_S(1:3,1:3,1:myNgrains,i,e) = &
crystallite_partionedS0(1:3,1:3,1:myNgrains,i,e)
do g = 1, myNgrains do g = 1, myNgrains
plasticState (material_phaseAt(g,e))%state( :,material_phasememberAt(g,i,e)) = & plasticState (material_phaseAt(g,e))%state( :,material_phasememberAt(g,i,e)) = &
plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e)) plasticState (material_phaseAt(g,e))%partionedState0(:,material_phasememberAt(g,i,e))
@ -386,15 +355,11 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = & damageState(material_homogenizationAt(e))%State( :,material_homogenizationMemberAt(i,e)) = &
damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e)) damageState(material_homogenizationAt(e))%subState0(:,material_homogenizationMemberAt(i,e))
endif endif
endif converged endif
if (materialpoint_subStep(i,e) > num%subStepMinHomog) then if (subStep(i,e) > num%subStepMinHomog) then
materialpoint_requested(i,e) = .true. requested(i,e) = .true.
materialpoint_subF(1:3,1:3,i,e) = materialpoint_subF0(1:3,1:3,i,e) & doneAndHappy(1:2,i,e) = [.false.,.true.]
+ materialpoint_subStep(i,e) * (materialpoint_F(1:3,1:3,i,e) &
- materialpoint_F0(1:3,1:3,i,e))
materialpoint_subdt(i,e) = materialpoint_subStep(i,e) * dt
materialpoint_doneAndHappy(1:2,i,e) = [.false.,.true.]
endif endif
enddo IpLooping1 enddo IpLooping1
enddo elementLooping1 enddo elementLooping1
@ -403,8 +368,8 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
NiterationMPstate = 0 NiterationMPstate = 0
convergenceLooping: do while (.not. terminallyIll .and. & convergenceLooping: do while (.not. terminallyIll .and. &
any( materialpoint_requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) & any( requested(:,FEsolving_execELem(1):FEsolving_execElem(2)) &
.and. .not. materialpoint_doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) & .and. .not. doneAndHappy(1,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
) .and. & ) .and. &
NiterationMPstate < num%nMPstate) NiterationMPstate < num%nMPstate)
NiterationMPstate = NiterationMPstate + 1 NiterationMPstate = NiterationMPstate + 1
@ -413,14 +378,15 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
! deformation partitioning ! deformation partitioning
! based on materialpoint_subF0,.._subF,crystallite_partionedF0, and homogenization_state, ! based on materialpoint_subF0,.._subF,crystallite_partionedF0, and homogenization_state,
! results in crystallite_partionedF ! results in crystallite_partionedF
!$OMP PARALLEL DO PRIVATE(myNgrains) !$OMP PARALLEL DO PRIVATE(myNgrains,subF)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Ngrains(material_homogenizationAt(e))
IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if ( materialpoint_requested(i,e) .and. & ! process requested but... if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done
.not. materialpoint_doneAndHappy(1,i,e)) then ! ...not yet done material points subF = materialpoint_F0(1:3,1:3,i,e) &
call partitionDeformation(i,e) ! partition deformation onto constituents + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))*(subStep(i,e)+subFrac(i,e))
crystallite_dt(1:myNgrains,i,e) = materialpoint_subdt(i,e) ! propagate materialpoint dt to grains call partitionDeformation(subF,i,e) ! partition deformation onto constituents
crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains
crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents crystallite_requested(1:myNgrains,i,e) = .true. ! request calculation for constituents
else else
crystallite_requested(1:myNgrains,i,e) = .false. ! calculation for constituents not required anymore crystallite_requested(1:myNgrains,i,e) = .false. ! calculation for constituents not required anymore
@ -434,20 +400,21 @@ subroutine materialpoint_stressAndItsTangent(updateJaco,dt)
! based on crystallite_partionedF0,.._partionedF ! based on crystallite_partionedF0,.._partionedF
! incrementing by crystallite_dt ! incrementing by crystallite_dt
materialpoint_converged = crystallite_stress() !ToDo: MD not sure if that is the best logic converged = crystallite_stress() !ToDo: MD not sure if that is the best logic
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state update ! state update
!$OMP PARALLEL DO !$OMP PARALLEL DO PRIVATE(subF)
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if ( materialpoint_requested(i,e) .and. & if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then
.not. materialpoint_doneAndHappy(1,i,e)) then if (.not. converged(i,e)) then
if (.not. materialpoint_converged(i,e)) then doneAndHappy(1:2,i,e) = [.true.,.false.]
materialpoint_doneAndHappy(1:2,i,e) = [.true.,.false.]
else else
materialpoint_doneAndHappy(1:2,i,e) = updateState(i,e) subF = materialpoint_F0(1:3,1:3,i,e) &
materialpoint_converged(i,e) = all(materialpoint_doneAndHappy(1:2,i,e)) ! converged if done and happy + (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))*(subStep(i,e)+subFrac(i,e))
doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e),subF,i,e)
converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy
endif endif
endif endif
enddo IpLooping3 enddo IpLooping3
@ -481,29 +448,31 @@ end subroutine materialpoint_stressAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief partition material point def grad onto constituents !> @brief partition material point def grad onto constituents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine partitionDeformation(ip,el) subroutine partitionDeformation(subF,ip,el)
integer, intent(in) :: & real(pReal), intent(in), dimension(3,3) :: &
ip, & !< integration point subF
el !< element number integer, intent(in) :: &
ip, & !< integration point
el !< element number
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
crystallite_partionedF(1:3,1:3,1,ip,el) = materialpoint_subF(1:3,1:3,ip,el) crystallite_partionedF(1:3,1:3,1,ip,el) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(& call mech_isostrain_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(&
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
materialpoint_subF(1:3,1:3,ip,el)) subF,&
ip, &
case (HOMOGENIZATION_RGC_ID) chosenHomogenization el)
call mech_RGC_partitionDeformation(& end select chosenHomogenization
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
materialpoint_subF(1:3,1:3,ip,el),&
ip, &
el)
end select chosenHomogenization
end subroutine partitionDeformation end subroutine partitionDeformation
@ -512,45 +481,49 @@ end subroutine partitionDeformation
!> @brief update the internal state of the homogenization scheme and tell whether "done" and !> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> "happy" with result !> "happy" with result
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function updateState(ip,el) function updateState(subdt,subF,ip,el)
integer, intent(in) :: & real(pReal), intent(in) :: &
ip, & !< integration point subdt !< current time step
el !< element number real(pReal), intent(in), dimension(3,3) :: &
logical, dimension(2) :: updateState subF
integer, intent(in) :: &
ip, & !< integration point
el !< element number
logical, dimension(2) :: updateState
updateState = .true. updateState = .true.
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
updateState = & updateState = &
updateState .and. & updateState .and. &
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),& crystallite_partionedF0(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el),&
materialpoint_subF(1:3,1:3,ip,el),& subF,&
materialpoint_subdt(ip,el), & subdt, &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
ip, & ip, &
el) el)
end select chosenHomogenization end select chosenHomogenization
chosenThermal: select case (thermal_type(material_homogenizationAt(el))) chosenThermal: select case (thermal_type(material_homogenizationAt(el)))
case (THERMAL_adiabatic_ID) chosenThermal case (THERMAL_adiabatic_ID) chosenThermal
updateState = & updateState = &
updateState .and. & updateState .and. &
thermal_adiabatic_updateState(materialpoint_subdt(ip,el), & thermal_adiabatic_updateState(subdt, &
ip, & ip, &
el) el)
end select chosenThermal end select chosenThermal
chosenDamage: select case (damage_type(material_homogenizationAt(el))) chosenDamage: select case (damage_type(material_homogenizationAt(el)))
case (DAMAGE_local_ID) chosenDamage case (DAMAGE_local_ID) chosenDamage
updateState = & updateState = &
updateState .and. & updateState .and. &
damage_local_updateState(materialpoint_subdt(ip,el), & damage_local_updateState(subdt, &
ip, & ip, &
el) el)
end select chosenDamage end select chosenDamage
end function updateState end function updateState
@ -560,31 +533,31 @@ end function updateState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine averageStressAndItsTangent(ip,el) subroutine averageStressAndItsTangent(ip,el)
integer, intent(in) :: & integer, intent(in) :: &
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el)
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_dPdF(1:3,1:3,1:3,1:3,1,ip,el) materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_dPdF(1:3,1:3,1:3,1:3,1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_averageStressAndItsTangent(& call mech_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_averageStressAndItsTangent(& call mech_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & materialpoint_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_dPdF(1:3,1:3,1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization end select chosenHomogenization
end subroutine averageStressAndItsTangent end subroutine averageStressAndItsTangent

11
src/rotations.f90
View File

@ -432,18 +432,17 @@ pure function qu2eu(qu) result(eu)
real(pReal), intent(in), dimension(4) :: qu real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3) :: eu real(pReal), dimension(3) :: eu
real(pReal) :: q12, q03, chi, chiInv real(pReal) :: q12, q03, chi
q03 = qu(1)**2+qu(4)**2 q03 = qu(1)**2+qu(4)**2
q12 = qu(2)**2+qu(3)**2 q12 = qu(2)**2+qu(3)**2
chi = sqrt(q03*q12) chi = sqrt(q03*q12)
degenerated: if (dEq0(chi)) then degenerated: if (dEq0(q12)) then
eu = merge([atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal], & eu = [atan2(-P*2.0_pReal*qu(1)*qu(4),qu(1)**2-qu(4)**2), 0.0_pReal, 0.0_pReal]
[atan2( 2.0_pReal*qu(2)*qu(3),qu(2)**2-qu(3)**2), PI, 0.0_pReal], & elseif (dEq0(q03)) then
dEq0(q12)) eu = [atan2( 2.0_pReal*qu(2)*qu(3),qu(2)**2-qu(3)**2), PI, 0.0_pReal]
else degenerated else degenerated
chiInv = 1.0_pReal/chi
eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), & eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), &
atan2( 2.0_pReal*chi, q03-q12 ), & atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )] atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]