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 56-parallel-hdf5

This commit is contained in:
Martin Diehl 2018-12-12 07:46:02 +01:00
parent fd2d4d856b 0ff9e9fa06
commit 3bab08fdeb
13 changed files with 202 additions and 254 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
.gitlab-ci.yml
View File

@ -158,12 +158,12 @@ Post_AverageDown:
- master - master
- release - release
Post_General: #Post_General:
stage: postprocessing # stage: postprocessing
script: PostProcessing/test.py # script: PostProcessing/test.py
except: # except:
- master # - master
- release # - release
Post_GeometryReconstruction: Post_GeometryReconstruction:
stage: postprocessing stage: postprocessing
@ -373,12 +373,12 @@ Phenopowerlaw_singleSlip:
- master - master
- release - release
TextureComponents: #TextureComponents:
stage: spectral # stage: spectral
script: TextureComponents/test.py # script: TextureComponents/test.py
except: # except:
- master # - master
- release # - release
################################################################################################### ###################################################################################################

3
DAMASK_prerequisites.sh
View File

@ -119,6 +119,9 @@ for executable in mpirun mpiexec; do
getDetails $executable '--version' getDetails $executable '--version'
done done
firstLevel "CMake"
getDetails cmake --version
firstLevel "Abaqus" firstLevel "Abaqus"
cd installation/mods_Abaqus # to have the right environment file cd installation/mods_Abaqus # to have the right environment file
for executable in abaqus abq2017 abq2018; do for executable in abaqus abq2017 abq2018; do

2
VERSION
View File

@ -1 +1 @@
v2.0.2-1107-g17716b4f v2.0.2-1122-g2349341e

360
lib/damask/orientation.py
View File

@ -33,7 +33,7 @@ class Quaternion:
All methods and naming conventions based on Rowenhorst_etal2015 All methods and naming conventions based on Rowenhorst_etal2015
Convention 1: coordinate frames are right-handed Convention 1: coordinate frames are right-handed
Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation
when viewing from the end point of the rotation axis unit vector towards the origin when viewing from the end point of the rotation axis towards the origin
Convention 3: rotations will be interpreted in the passive sense Convention 3: rotations will be interpreted in the passive sense
Convention 4: Euler angle triplets are implemented using the Bunge convention, Convention 4: Euler angle triplets are implemented using the Bunge convention,
with the angular ranges as [0, 2π],[0, π],[0, 2π] with the angular ranges as [0, 2π],[0, π],[0, 2π]
@ -48,206 +48,138 @@ class Quaternion:
""" """
def __init__(self, def __init__(self,
quatArray = [1.0,0.0,0.0,0.0]): quat = None,
"""Initializes to identity if not given""" q = 1.0,
self.w, \ p = np.zeros(3,dtype=float)):
self.x, \ """Initializes to identity unless specified"""
self.y, \ self.q = quat[0] if quat is not None else q
self.z = quatArray self.p = np.array(quat[1:4]) if quat is not None else p
self.homomorph() self.homomorph()
def __iter__(self): def __iter__(self):
"""Components""" """Components"""
return iter([self.w,self.x,self.y,self.z]) return iter(self.asList())
def __copy__(self): def __copy__(self):
"""Create copy""" """Copy"""
Q = Quaternion([self.w,self.x,self.y,self.z]) return self.__class__(q=self.q,p=self.p.copy())
return Q
copy = __copy__ copy = __copy__
def __repr__(self): def __repr__(self):
"""Readbable string""" """Readable string"""
return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \ return 'Quaternion(real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p)
(self.w, self.x, self.y, self.z)
def __pow__(self, exponent): def __pow__(self, exponent):
"""Power""" """Power"""
omega = math.acos(self.w) omega = math.acos(self.q)
vRescale = math.sin(exponent*omega)/math.sin(omega) return self.__class__(q= math.cos(exponent*omega),
Q = Quaternion() p=self.p * math.sin(exponent*omega)/math.sin(omega))
Q.w = math.cos(exponent*omega)
Q.x = self.x * vRescale
Q.y = self.y * vRescale
Q.z = self.z * vRescale
return Q
def __ipow__(self, exponent): def __ipow__(self, exponent):
"""In-place power""" """In-place power"""
omega = math.acos(self.w) omega = math.acos(self.q)
vRescale = math.sin(exponent*omega)/math.sin(omega) self.q = math.cos(exponent*omega)
self.w = np.cos(exponent*omega) self.p *= math.sin(exponent*omega)/math.sin(omega)
self.x *= vRescale
self.y *= vRescale
self.z *= vRescale
return self return self
def __mul__(self, other): def __mul__(self, other):
"""Multiplication""" """Multiplication"""
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
try: # quaternion try: # quaternion
Aw = self.w return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
Ax = self.x p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
Ay = self.y
Az = self.z
Bw = other.w
Bx = other.x
By = other.y
Bz = other.z
Q = Quaternion()
Q.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw
Q.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx
Q.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By
Q.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz
return Q
except: pass except: pass
try: # vector (perform active rotation, i.e. q*v*q.conjugated) try: # vector (perform passive rotation)
w = self.w ( x, y, z) = self.p
x = self.x (Vx,Vy,Vz) = other[0:3]
y = self.y A = self.q*self.q - np.dot(self.p,self.p)
z = self.z B = 2.0 * (x*Vx + y*Vy + z*Vz)
Vx = other[0] C = 2.0 * P*self.q
Vy = other[1]
Vz = other[2]
return np.array([\ return np.array([
w * w * Vx + 2 * y * w * Vz - 2 * z * w * Vy + \ A*Vx + B*x + C*(y*Vz - z*Vy),
x * x * Vx + 2 * y * x * Vy + 2 * z * x * Vz - \ A*Vy + B*y + C*(z*Vx - x*Vz),
z * z * Vx - y * y * Vx, A*Vz + B*z + C*(x*Vy - y*Vx),
2 * x * y * Vx + y * y * Vy + 2 * z * y * Vz + \ ])
2 * w * z * Vx - z * z * Vy + w * w * Vy - \
2 * x * w * Vz - x * x * Vy,
2 * x * z * Vx + 2 * y * z * Vy + \
z * z * Vz - 2 * w * y * Vx - y * y * Vz + \
2 * w * x * Vy - x * x * Vz + w * w * Vz ])
except: pass except: pass
try: # scalar try: # scalar
Q = self.copy() return self.__class__(q=self.q*other,
Q.w *= other p=self.p*other)
Q.x *= other
Q.y *= other
Q.z *= other
return Q
except: except:
return self.copy() return self.copy()
def __imul__(self, other): def __imul__(self, other):
"""In-place multiplication""" """In-place multiplication"""
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
try: # Quaternion try: # Quaternion
Aw = self.w self.q = self.q*other.q - np.dot(self.p,other.p)
Ax = self.x self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
Ay = self.y
Az = self.z
Bw = other.w
Bx = other.x
By = other.y
Bz = other.z
self.w = - Ax * Bx - Ay * By - Az * Bz + Aw * Bw
self.x = + Ax * Bw + Ay * Bz - Az * By + Aw * Bx
self.y = - Ax * Bz + Ay * Bw + Az * Bx + Aw * By
self.z = + Ax * By - Ay * Bx + Az * Bw + Aw * Bz
except: pass except: pass
return self return self
def __div__(self, other): def __div__(self, other):
"""Division""" """Division"""
if isinstance(other, (int,float)): if isinstance(other, (int,float)):
w = self.w / other return self.__class__(q=self.q / other,
x = self.x / other p=self.p / other)
y = self.y / other
z = self.z / other
return self.__class__([w,x,y,z])
else: else:
return NotImplemented return NotImplemented
def __idiv__(self, other): def __idiv__(self, other):
"""In-place division""" """In-place division"""
if isinstance(other, (int,float)): if isinstance(other, (int,float)):
self.w /= other self.q /= other
self.x /= other self.p /= other
self.y /= other
self.z /= other
return self return self
def __add__(self, other): def __add__(self, other):
"""Addition""" """Addition"""
if isinstance(other, Quaternion): if isinstance(other, Quaternion):
w = self.w + other.w return self.__class__(q=self.q + other.q,
x = self.x + other.x p=self.p + other.p)
y = self.y + other.y
z = self.z + other.z
return self.__class__([w,x,y,z])
else: else:
return NotImplemented return NotImplemented
def __iadd__(self, other): def __iadd__(self, other):
"""In-place addition""" """In-place addition"""
if isinstance(other, Quaternion): if isinstance(other, Quaternion):
self.w += other.w self.q += other.q
self.x += other.x self.p += other.p
self.y += other.y
self.z += other.z
return self return self
def __sub__(self, other): def __sub__(self, other):
"""Subtraction""" """Subtraction"""
if isinstance(other, Quaternion): if isinstance(other, Quaternion):
Q = self.copy() return self.__class__(q=self.q - other.q,
Q.w -= other.w p=self.p - other.p)
Q.x -= other.x
Q.y -= other.y
Q.z -= other.z
return Q
else: else:
return self.copy() return NotImplemented
def __isub__(self, other): def __isub__(self, other):
"""In-place subtraction""" """In-place subtraction"""
if isinstance(other, Quaternion): if isinstance(other, Quaternion):
self.w -= other.w self.q -= other.q
self.x -= other.x self.p -= other.p
self.y -= other.y
self.z -= other.z
return self return self
def __neg__(self): def __neg__(self):
"""Additive inverse""" """Additive inverse"""
self.w = -self.w self.q = -self.q
self.x = -self.x self.p = -self.p
self.y = -self.y
self.z = -self.z
return self return self
def __abs__(self): def __abs__(self):
"""Norm""" """Norm"""
return math.sqrt(self.w ** 2 + \ return math.sqrt(self.q ** 2 + np.dot(self.p,self.p))
self.x ** 2 + \
self.y ** 2 + \
self.z ** 2)
magnitude = __abs__ magnitude = __abs__
def __eq__(self,other): def __eq__(self,other):
"""Equal at e-8 precision""" """Equal at e-8 precision"""
return (abs(self.w-other.w) < 1e-8 and \ return (self-other).magnitude() < 1e-8 or (-self-other).magnitude() < 1e-8
abs(self.x-other.x) < 1e-8 and \
abs(self.y-other.y) < 1e-8 and \
abs(self.z-other.z) < 1e-8) \
or \
(abs(-self.w-other.w) < 1e-8 and \
abs(-self.x-other.x) < 1e-8 and \
abs(-self.y-other.y) < 1e-8 and \
abs(-self.z-other.z) < 1e-8)
def __ne__(self,other): def __ne__(self,other):
"""Not equal at e-8 precision""" """Not equal at e-8 precision"""
@ -255,31 +187,26 @@ class Quaternion:
def __cmp__(self,other): def __cmp__(self,other):
"""Linear ordering""" """Linear ordering"""
return (self.Rodrigues()>other.Rodrigues()) - (self.Rodrigues()<other.Rodrigues()) return (1 if np.linalg.norm(self.asRodrigues()) > np.linalg.norm(other.asRodrigues()) else 0) \
- (1 if np.linalg.norm(self.asRodrigues()) < np.linalg.norm(other.asRodrigues()) else 0)
def magnitude_squared(self): def magnitude_squared(self):
return self.w ** 2 + \ return self.q ** 2 + np.dot(self.p,self.p)
self.x ** 2 + \
self.y ** 2 + \
self.z ** 2
def identity(self): def identity(self):
self.w = 1. self.q = 1.
self.x = 0. self.p = np.zeros(3,dtype=float)
self.y = 0.
self.z = 0.
return self return self
def normalize(self): def normalize(self):
d = self.magnitude() d = self.magnitude()
if d > 0.0: if d > 0.0:
self /= d self.q /= d
self.p /= d
return self return self
def conjugate(self): def conjugate(self):
self.x = -self.x self.p = -self.p
self.y = -self.y
self.z = -self.z
return self return self
def inverse(self): def inverse(self):
@ -290,11 +217,9 @@ class Quaternion:
return self return self
def homomorph(self): def homomorph(self):
if self.w < 0.0: if self.q < 0.0:
self.w = -self.w self.q = -self.q
self.x = -self.x self.p = -self.p
self.y = -self.y
self.z = -self.z
return self return self
def normalized(self): def normalized(self):
@ -310,27 +235,35 @@ class Quaternion:
return self.copy().homomorph() return self.copy().homomorph()
def asList(self): def asList(self):
return [i for i in self] return [self.q]+list(self.p)
def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.) def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.)
return np.outer([i for i in self],[i for i in self]) return np.outer(self.asList(),self.asList())
def asMatrix(self): def asMatrix(self):
qbarhalf = 0.5*(self.w**2 - self.x**2 - self.y**2 - self.z**2) # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
qbarhalf = 0.5*(self.q**2 - np.dot(self.p,self.p))
return 2.0*np.array( return 2.0*np.array(
[[ qbarhalf + self.x**2 , self.x*self.y - self.w*self.z, self.x*self.z + self.w*self.y], [[ qbarhalf + self.p[0]**2 ,
[ self.x*self.y + self.w*self.z, qbarhalf + self.y**2 , self.y*self.z - self.w*self.x], self.p[0]*self.p[1] -P* self.q*self.p[2],
[ self.x*self.z - self.w*self.y, self.y*self.z + self.w*self.x, qbarhalf + self.z**2 ], self.p[0]*self.p[2] +P* self.q*self.p[1] ],
[ self.p[0]*self.p[1] +P* self.q*self.p[2],
qbarhalf + self.p[1]**2 ,
self.p[1]*self.p[2] -P* self.q*self.p[0] ],
[ self.p[0]*self.p[2] -P* self.q*self.p[1],
self.p[1]*self.p[2] +P* self.q*self.p[0],
qbarhalf + self.p[2]**2 ],
]) ])
def asAngleAxis(self, def asAngleAxis(self,
degrees = False): degrees = False):
if self.w > 1: if self.q > 1.:
self.normalize() self.normalize()
s = math.sqrt(1. - self.w**2) s = math.sqrt(1. - self.q**2)
x = 2*self.w**2 - 1. x = 2*self.q**2 - 1.
y = 2*self.w * s y = 2*self.q * s
angle = math.atan2(y,x) angle = math.atan2(y,x)
if angle < 0.0: if angle < 0.0:
@ -338,26 +271,28 @@ class Quaternion:
s *= -1. s *= -1.
return (np.degrees(angle) if degrees else angle, return (np.degrees(angle) if degrees else angle,
np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else [self.x / s, self.y / s, self.z / s])) np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s))
def asRodrigues(self): def asRodrigues(self):
return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w return np.inf*np.ones(3) if self.q == 0.0 else self.p/self.q
def asEulers(self, def asEulers(self,
degrees = False): degrees = False):
"""Orientation as Bunge-Euler angles.""" """Orientation as Bunge-Euler angles."""
q03 = self.w**2+self.z**2 # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
q12 = self.x**2+self.y**2 P = -1.0
q03 = self.q**2 + self.p[2]**2
q12 = self.p[0]**2 + self.p[1]**2
chi = np.sqrt(q03*q12) chi = np.sqrt(q03*q12)
if abs(chi) < 1e-10 and abs(q12) < 1e-10: if abs(chi) < 1e-10 and abs(q12) < 1e-10:
eulers = np.array([math.atan2(-2*self.w*self.z,self.w**2-self.z**2),0,0]) eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0])
elif abs(chi) < 1e-10 and abs(q03) < 1e-10: elif abs(chi) < 1e-10 and abs(q03) < 1e-10:
eulers = np.array([math.atan2( 2*self.x*self.y,self.x**2-self.y**2),np.pi,0]) eulers = np.array([math.atan2( 2 *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0])
else: else:
eulers = np.array([math.atan2((self.x*self.z-self.w*self.y)/chi,(-self.w*self.x-self.y*self.z)/chi), eulers = np.array([math.atan2((self.p[0]*self.p[2]-P*self.q*self.p[1])/chi,(-P*self.q*self.p[0]-self.p[1]*self.p[2])/chi),
math.atan2(2*chi,q03-q12), math.atan2(2*chi,q03-q12),
math.atan2((self.w*self.y+self.x*self.z)/chi,( self.y*self.z-self.w*self.x)/chi), math.atan2((P*self.q*self.p[1]+self.p[0]*self.p[2])/chi,( self.p[1]*self.p[2]-P*self.q*self.p[0])/chi),
]) ])
return np.degrees(eulers) if degrees else eulers return np.degrees(eulers) if degrees else eulers
@ -371,25 +306,26 @@ class Quaternion:
@classmethod @classmethod
def fromRandom(cls,randomSeed = None): def fromRandom(cls,randomSeed = None):
import binascii
if randomSeed is None: if randomSeed is None:
randomSeed = int(os.urandom(4).encode('hex'), 16) randomSeed = int(binascii.hexlify(os.urandom(4)),16)
np.random.seed(randomSeed) np.random.seed(randomSeed)
r = np.random.random(3) r = np.random.random(3)
w = math.cos(2.0*math.pi*r[0])*math.sqrt(r[2]) w = math.cos(2.0*math.pi*r[0])*math.sqrt(r[2])
x = math.sin(2.0*math.pi*r[1])*math.sqrt(1.0-r[2]) x = math.sin(2.0*math.pi*r[1])*math.sqrt(1.0-r[2])
y = math.cos(2.0*math.pi*r[1])*math.sqrt(1.0-r[2]) y = math.cos(2.0*math.pi*r[1])*math.sqrt(1.0-r[2])
z = math.sin(2.0*math.pi*r[0])*math.sqrt(r[2]) z = math.sin(2.0*math.pi*r[0])*math.sqrt(r[2])
return cls([w,x,y,z]) return cls(quat=[w,x,y,z])
@classmethod @classmethod
def fromRodrigues(cls, rodrigues): def fromRodrigues(cls, rodrigues):
if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues) if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues)
halfangle = math.atan(np.linalg.norm(rodrigues)) norm = np.linalg.norm(rodrigues)
halfangle = math.atan(norm)
s = math.sin(halfangle)
c = math.cos(halfangle) c = math.cos(halfangle)
w = c return cls(q=c,p=s*rodrigues/norm)
x,y,z = rodrigues/c
return cls([w,x,y,z])
@classmethod @classmethod
@ -397,22 +333,19 @@ class Quaternion:
angle, angle,
axis, axis,
degrees = False): degrees = False):
if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype='d') if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype=float)
axis = axis.astype(float)/np.linalg.norm(axis) axis = axis.astype(float)/np.linalg.norm(axis)
angle = np.radians(angle) if degrees else angle angle = np.radians(angle) if degrees else angle
s = math.sin(0.5 * angle) s = math.sin(0.5 * angle)
w = math.cos(0.5 * angle) c = math.cos(0.5 * angle)
x = axis[0] * s return cls(q=c,p=axis*s)
y = axis[1] * s
z = axis[2] * s
return cls([w,x,y,z])
@classmethod @classmethod
def fromEulers(cls, def fromEulers(cls,
eulers, eulers,
degrees = False): degrees = False):
if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d') if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype=float)
eulers = np.radians(eulers) if degrees else eulers eulers = np.radians(eulers) if degrees else eulers
sigma = 0.5*(eulers[0]+eulers[2]) sigma = 0.5*(eulers[0]+eulers[2])
@ -420,11 +353,13 @@ class Quaternion:
c = np.cos(0.5*eulers[1]) c = np.cos(0.5*eulers[1])
s = np.sin(0.5*eulers[1]) s = np.sin(0.5*eulers[1])
w = c * np.cos(sigma) # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
x = -s * np.cos(delta) P = -1.0
y = -s * np.sin(delta) w = c * np.cos(sigma)
z = -c * np.sin(sigma) x = -P * s * np.cos(delta)
return cls([w,x,y,z]) y = -P * s * np.sin(delta)
z = -P * c * np.sin(sigma)
return cls(quat=[w,x,y,z])
# Modified Method to calculate Quaternion from Orientation Matrix, # Modified Method to calculate Quaternion from Orientation Matrix,
@ -435,16 +370,18 @@ class Quaternion:
if m.shape != (3,3) and np.prod(m.shape) == 9: if m.shape != (3,3) and np.prod(m.shape) == 9:
m = m.reshape(3,3) m = m.reshape(3,3)
w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2]) # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
x = 0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2]) P = -1.0
y = 0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2]) w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2])
z = 0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2]) x = P*0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
y = P*0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2])
z = P*0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2])
x *= -1 if m[2,1] < m[1,2] else 1 x *= -1 if m[2,1] < m[1,2] else 1
y *= -1 if m[0,2] < m[2,0] else 1 y *= -1 if m[0,2] < m[2,0] else 1
z *= -1 if m[1,0] < m[0,1] else 1 z *= -1 if m[1,0] < m[0,1] else 1
return cls( np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2)) return cls(quat=np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2))
@classmethod @classmethod
@ -458,36 +395,30 @@ class Quaternion:
assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion) assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion)
Q = cls() Q = cls()
costheta = q1.w * q2.w + q1.x * q2.x + q1.y * q2.y + q1.z * q2.z costheta = q1.q*q2.q + np.dot(q1.p,q2.p)
if costheta < 0.: if costheta < 0.:
costheta = -costheta costheta = -costheta
q1 = q1.conjugated() q1 = q1.conjugated()
elif costheta > 1: elif costheta > 1.:
costheta = 1 costheta = 1.
theta = math.acos(costheta) theta = math.acos(costheta)
if abs(theta) < 0.01: if abs(theta) < 0.01:
Q.w = q2.w Q.q = q2.q
Q.x = q2.x Q.p = q2.p
Q.y = q2.y
Q.z = q2.z
return Q return Q
sintheta = math.sqrt(1.0 - costheta * costheta) sintheta = math.sqrt(1.0 - costheta * costheta)
if abs(sintheta) < 0.01: if abs(sintheta) < 0.01:
Q.w = (q1.w + q2.w) * 0.5 Q.q = (q1.q + q2.q) * 0.5
Q.x = (q1.x + q2.x) * 0.5 Q.p = (q1.p + q2.p) * 0.5
Q.y = (q1.y + q2.y) * 0.5
Q.z = (q1.z + q2.z) * 0.5
return Q return Q
ratio1 = math.sin((1 - t) * theta) / sintheta ratio1 = math.sin((1.0 - t) * theta) / sintheta
ratio2 = math.sin(t * theta) / sintheta ratio2 = math.sin( t * theta) / sintheta
Q.w = q1.w * ratio1 + q2.w * ratio2 Q.q = q1.q * ratio1 + q2.q * ratio2
Q.x = q1.x * ratio1 + q2.x * ratio2 Q.p = q1.p * ratio1 + q2.p * ratio2
Q.y = q1.y * ratio1 + q2.y * ratio2
Q.z = q1.z * ratio1 + q2.z * ratio2
return Q return Q
@ -513,7 +444,7 @@ class Symmetry:
def __repr__(self): def __repr__(self):
"""Readbable string""" """Readbable string"""
return '%s' % (self.lattice) return '{}'.format(self.lattice)
def __eq__(self, other): def __eq__(self, other):
@ -526,7 +457,7 @@ class Symmetry:
def __cmp__(self,other): def __cmp__(self,other):
"""Linear ordering""" """Linear ordering"""
myOrder = Symmetry.lattices.index(self.lattice) myOrder = Symmetry.lattices.index(self.lattice)
otherOrder = Symmetry.lattices.index(other.lattice) otherOrder = Symmetry.lattices.index(other.lattice)
return (myOrder > otherOrder) - (myOrder < otherOrder) return (myOrder > otherOrder) - (myOrder < otherOrder)
@ -722,7 +653,7 @@ class Symmetry:
else: else:
return True return True
v = np.array(vector,dtype = float) v = np.array(vector,dtype=float)
if proper: # check both improper ... if proper: # check both improper ...
theComponents = np.dot(basis['improper'],v) theComponents = np.dot(basis['improper'],v)
inSST = np.all(theComponents >= 0.0) inSST = np.all(theComponents >= 0.0)
@ -737,10 +668,10 @@ class Symmetry:
if color: # have to return color array if color: # have to return color array
if inSST: if inSST:
rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps
rgb = np.minimum(np.ones(3,'d'),rgb) # limit to maximum intensity rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity
rgb /= max(rgb) # normalize to (HS)V = 1 rgb /= max(rgb) # normalize to (HS)V = 1
else: else:
rgb = np.zeros(3,'d') rgb = np.zeros(3,dtype=float)
return (inSST,rgb) return (inSST,rgb)
else: else:
return inSST return inSST
@ -780,8 +711,9 @@ class Orientation:
self.quaternion = Quaternion.fromRodrigues(Rodrigues) self.quaternion = Quaternion.fromRodrigues(Rodrigues)
elif isinstance(quaternion, Quaternion): # based on given quaternion elif isinstance(quaternion, Quaternion): # based on given quaternion
self.quaternion = quaternion.homomorphed() self.quaternion = quaternion.homomorphed()
elif isinstance(quaternion, np.ndarray) and quaternion.shape == (4,): # based on given quaternion-like array elif (isinstance(quaternion, np.ndarray) and quaternion.shape == (4,)) or \
self.quaternion = Quaternion(quaternion).homomorphed() (isinstance(quaternion, list) and len(quaternion) == 4 ): # based on given quaternion-like array
self.quaternion = Quaternion(quat=quaternion).homomorphed()
self.symmetry = Symmetry(symmetry) self.symmetry = Symmetry(symmetry)
@ -794,10 +726,12 @@ class Orientation:
def __repr__(self): def __repr__(self):
"""Value as all implemented representations""" """Value as all implemented representations"""
return 'Symmetry: %s\n' % (self.symmetry) + \ return '\n'.join([
'Quaternion: %s\n' % (self.quaternion) + \ 'Symmetry: {}'.format(self.symmetry),
'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \ 'Quaternion: {}'.format(self.quaternion),
'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers(degrees=True))) ) 'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ),
'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ),
])
def asQuaternion(self): def asQuaternion(self):
return self.quaternion.asList() return self.quaternion.asList()
@ -927,7 +861,7 @@ class Orientation:
M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa
eig, vec = np.linalg.eig(M/N) eig, vec = np.linalg.eig(M/N)
return Orientation(quaternion = Quaternion(quatArray = np.real(vec.T[eig.argmax()])), return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])),
symmetry = reference.symmetry.lattice) symmetry = reference.symmetry.lattice)

4
processing/post/addGrainID.py
View File

@ -200,9 +200,9 @@ for name in filenames:
if gID != -1 and gID not in alreadyChecked: # indexed point belonging to a grain not yet tested? if gID != -1 and gID not in alreadyChecked: # indexed point belonging to a grain not yet tested?
alreadyChecked[gID] = True # remember not to check again alreadyChecked[gID] = True # remember not to check again
disorientation = o.disorientation(orientations[gID],SST = False)[0] # compare against other orientation disorientation = o.disorientation(orientations[gID],SST = False)[0] # compare against other orientation
if disorientation.quaternion.w > cos_disorientation: # within threshold ... if disorientation.quaternion.q > cos_disorientation: # within threshold ...
candidates.append(gID) # remember as potential candidate candidates.append(gID) # remember as potential candidate
if disorientation.quaternion.w >= bestDisorientation.w: # ... and better than current best? if disorientation.quaternion.q >= bestDisorientation.q: # ... and better than current best?
matched = True matched = True
matchedID = gID # remember that grain matchedID = gID # remember that grain
bestDisorientation = disorientation.quaternion bestDisorientation = disorientation.quaternion

8
processing/pre/geom_addPrimitive.py
View File

@ -64,11 +64,11 @@ if options.dimension is None:
parser.error('no dimension specified.') parser.error('no dimension specified.')
if options.angleaxis is not None: if options.angleaxis is not None:
options.angleaxis = list(map(float,options.angleaxis)) options.angleaxis = list(map(float,options.angleaxis))
rotation = damask.Quaternion().fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0], rotation = damask.Quaternion.fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0],
options.angleaxis[1:4]) options.angleaxis[1:4])
elif options.quaternion is not None: elif options.quaternion is not None:
options.quaternion = map(float,options.quaternion) options.quaternion = list(map(float,options.quaternion))
rotation = damask.Quaternion(options.quaternion) rotation = damask.Quaternion(quat=options.quaternion)
else: else:
rotation = damask.Quaternion() rotation = damask.Quaternion()

12
processing/pre/geom_fromTable.py
View File

@ -43,7 +43,7 @@ parser.add_option('-e', '--eulers',
parser.add_option('-d', '--degrees', parser.add_option('-d', '--degrees',
dest = 'degrees', dest = 'degrees',
action = 'store_true', action = 'store_true',
help = 'angles are given in degrees [%default]') help = 'all angles are in degrees')
parser.add_option('-m', '--matrix', parser.add_option('-m', '--matrix',
dest = 'matrix', dest = 'matrix',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
@ -71,7 +71,7 @@ parser.add_option('--axes',
parser.add_option('-s', '--symmetry', parser.add_option('-s', '--symmetry',
dest = 'symmetry', dest = 'symmetry',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'crystal symmetry %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:]))) help = 'crystal symmetry of each phase %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
parser.add_option('--homogenization', parser.add_option('--homogenization',
dest = 'homogenization', dest = 'homogenization',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
@ -234,7 +234,7 @@ for name in filenames:
o = damask.Orientation(Eulers = myData[colOri:colOri+3]*toRadians, o = damask.Orientation(Eulers = myData[colOri:colOri+3]*toRadians,
symmetry = mySym) symmetry = mySym)
elif inputtype == 'matrix': elif inputtype == 'matrix':
o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3).transpose(), o = damask.Orientation(matrix = myData[colOri:colOri+9].reshape(3,3),
symmetry = mySym) symmetry = mySym)
elif inputtype == 'frame': elif inputtype == 'frame':
o = damask.Orientation(matrix = np.hstack((myData[colOri[0]:colOri[0]+3], o = damask.Orientation(matrix = np.hstack((myData[colOri[0]:colOri[0]+3],
@ -246,7 +246,7 @@ for name in filenames:
o = damask.Orientation(quaternion = myData[colOri:colOri+4], o = damask.Orientation(quaternion = myData[colOri:colOri+4],
symmetry = mySym) symmetry = mySym)
cos_disorientations = -np.ones(1,dtype='f') # largest possible disorientation cos_disorientations = -np.ones(1,dtype=float) # largest possible disorientation
closest_grain = -1 # invalid neighbor closest_grain = -1 # invalid neighbor
if options.tolerance > 0.0: # only try to compress orientations if asked to if options.tolerance > 0.0: # only try to compress orientations if asked to
@ -258,7 +258,7 @@ for name in filenames:
if len(grains) > 0: # check immediate neighborhood first if len(grains) > 0: # check immediate neighborhood first
cos_disorientations = np.array([o.disorientation(orientations[grainID], cos_disorientations = np.array([o.disorientation(orientations[grainID],
SST = False)[0].quaternion.w \ SST = False)[0].quaternion.q \
for grainID in grains]) # store disorientation per grainID for grainID in grains]) # store disorientation per grainID
closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself
match = 'local' match = 'local'
@ -269,7 +269,7 @@ for name in filenames:
if len(grains) > 0: if len(grains) > 0:
cos_disorientations = np.array([o.disorientation(orientations[grainID], cos_disorientations = np.array([o.disorientation(orientations[grainID],
SST = False)[0].quaternion.w \ SST = False)[0].quaternion.q \
for grainID in grains]) # store disorientation per grainID for grainID in grains]) # store disorientation per grainID
closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself
match = 'global' match = 'global'

2
processing/pre/hybridIA_linODFsampling.py
View File

@ -244,7 +244,7 @@ for name in filenames:
continue continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
randomSeed = int(os.urandom(4).encode('hex'), 16) if options.randomSeed is None else options.randomSeed # random seed per file for second phase randomSeed = int(os.urandom(4).hex(), 16) if options.randomSeed is None else options.randomSeed # random seed per file for second phase
random.seed(randomSeed) random.seed(randomSeed)
# ------------------------------------------ read header and data --------------------------------- # ------------------------------------------ read header and data ---------------------------------

10
src/IO.f90
View File

@ -191,7 +191,9 @@ recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
l,i, & l,i, &
myStat myStat
if (merge(cnt,0_pInt,present(cnt))>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName)) if (present(cnt)) then
if (cnt>10_pInt) call IO_error(106_pInt,ext_msg=trim(fileName))
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! read data as stream ! read data as stream
@ -684,7 +686,11 @@ function IO_stringValue(string,chunkPos,myChunk,silent)
logical :: warn logical :: warn
warn = merge(silent,.false.,present(silent)) if (present(silent)) then
warn = silent
else
warn = .false.
endif
IO_stringValue = '' IO_stringValue = ''
valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then

12
src/config.f90
View File

@ -513,8 +513,12 @@ character(len=65536) function getString(this,key,defaultVal,raw)
type(tPartitionedStringList), pointer :: item type(tPartitionedStringList), pointer :: item
logical :: found, & logical :: found, &
whole whole
if (present(raw)) then
whole = raw
else
whole = .false.
endif
whole = merge(raw,.false.,present(raw)) ! whole string or white space splitting
found = present(defaultVal) found = present(defaultVal)
if (found) then if (found) then
getString = trim(defaultVal) getString = trim(defaultVal)
@ -661,7 +665,11 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
cumulative cumulative
cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')') cumulative = (key(1:1) == '(' .and. key(len_trim(key):len_trim(key)) == ')')
whole = merge(raw,.false.,present(raw)) if (present(raw)) then
whole = raw
else
whole = .false.
endif
found = .false. found = .false.
item => this item => this

2
src/kinematics_thermal_expansion.f90
View File

@ -225,7 +225,7 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc,
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient
) / & ) / &
(1.0_pReal \ (1.0_pReal &
+ lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. & + lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**1 / 1. &
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. & + lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**2 / 2. &
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. & + lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**3 / 3. &

11
src/math.f90
View File

@ -302,7 +302,7 @@ subroutine math_check
endif endif
end subroutine math_check end subroutine math_check
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Quicksort algorithm for two-dimensional integer arrays !> @brief Quicksort algorithm for two-dimensional integer arrays
@ -2625,12 +2625,9 @@ real(pReal) pure function math_clip(a, left, right)
real(pReal), intent(in) :: a real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: left, right real(pReal), intent(in), optional :: left, right
math_clip = a
math_clip = min ( & if (present(left)) math_clip = max(left,math_clip)
max (merge(left, -huge(a), present(left)), a), & if (present(right)) math_clip = min(right,math_clip)
merge(right, huge(a), present(right)) &
)
if (present(left) .and. present(right)) & if (present(left) .and. present(right)) &
math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right) math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right)

6
src/spectral_damage.f90
View File

@ -120,8 +120,8 @@ subroutine spectral_damage_init()
trim(snes_type) == 'vinewtonssls') then trim(snes_type) == 'vinewtonssls') then
call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
call VecSet(lBound,0.0,ierr); CHKERRQ(ierr) call VecSet(lBound,0.0_pReal,ierr); CHKERRQ(ierr)
call VecSet(uBound,1.0,ierr); CHKERRQ(ierr) call VecSet(uBound,1.0_pReal,ierr); CHKERRQ(ierr)
call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc. call SNESVISetVariableBounds(damage_snes,lBound,uBound,ierr) !< variable bounds for variational inequalities like contact mechanics, damage etc.
call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr) call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
@ -134,7 +134,7 @@ subroutine spectral_damage_init()
xend = xstart + xend - 1 xend = xstart + xend - 1
yend = ystart + yend - 1 yend = ystart + yend - 1
zend = zstart + zend - 1 zend = zstart + zend - 1
call VecSet(solution,1.0,ierr); CHKERRQ(ierr) call VecSet(solution,1.0_pReal,ierr); CHKERRQ(ierr)
allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_current(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal) allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)