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 22-NewStyle_disloUCLA-2

This commit is contained in:
Martin Diehl 2018-12-21 06:46:35 +01:00
parent fe1183e010 9bdd4d1d68
commit d4c7e8f33b
41 changed files with 605 additions and 956 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
- release
Post_General:
stage: postprocessing
script: PostProcessing/test.py
except:
- master
- release
#Post_General:
# stage: postprocessing
# script: PostProcessing/test.py
# except:
# - master
# - release
Post_GeometryReconstruction:
stage: postprocessing
@ -364,12 +364,12 @@ Phenopowerlaw_singleSlip:
- master
- release
TextureComponents:
stage: spectral
script: TextureComponents/test.py
except:
- master
- release
#TextureComponents:
# stage: spectral
# script: TextureComponents/test.py
# except:
# - master
# - release
###################################################################################################

5
DAMASK_prerequisites.sh
View File

@ -79,7 +79,7 @@ ls $PETSC_DIR/lib
firstLevel "Python"
DEFAULT_PYTHON=python3
for executable in python python2 python3 python2.7; do
for executable in python python3; do
getDetails $executable '--version'
done
secondLevel "Details on $DEFAULT_PYTHON:"
@ -119,6 +119,9 @@ for executable in mpirun mpiexec; do
getDetails $executable '--version'
done
firstLevel "CMake"
getDetails cmake --version
firstLevel "Abaqus"
cd installation/mods_Abaqus # to have the right environment file
for executable in abaqus abq2017 abq2018; do

2
PRIVATE

@ -1 +1 @@
Subproject commit e9f93abaecafbfbf11072ae70bca213a7201ed38
Subproject commit b9a52a85cd65cc27a8e863302bd984abdcad1455

2
VERSION
View File

@ -1 +1 @@
v2.0.2-1100-g65ff2157
v2.0.2-1187-gcd8ee450

7
env/DAMASK.csh vendored
View File

@ -7,6 +7,11 @@ set DAMASK_ROOT=`python -c "import os,sys; print(os.path.realpath(os.path.expand
source $DAMASK_ROOT/CONFIG
# add BRANCH if DAMASK_ROOT is a git repository
cd $DAMASK_ROOT >/dev/null
set BRANCH = `git branch 2>/dev/null| grep -E '^\* ')`
cd - >/dev/null
# if DAMASK_BIN is present
if ( $?DAMASK_BIN) then
set path = ($DAMASK_BIN $path)
@ -41,7 +46,7 @@ if ( $?prompt ) then
echo https://damask.mpie.de
echo
echo Using environment with ...
echo "DAMASK $DAMASK_ROOT"
echo "DAMASK $DAMASK_ROOT $BRANCH"
echo "Spectral Solver $SOLVER"
echo "Post Processing $PROCESSING"
if ( $?PETSC_DIR) then

7
env/DAMASK.sh vendored
View File

@ -30,6 +30,9 @@ set() {
source $DAMASK_ROOT/CONFIG
unset -f set
# add BRANCH if DAMASK_ROOT is a git repository
cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd - >/dev/null
# add DAMASK_BIN if present
[ "x$DAMASK_BIN" != "x" ] && PATH=$DAMASK_BIN:$PATH
@ -59,7 +62,7 @@ if [ ! -z "$PS1" ]; then
echo https://damask.mpie.de
echo
echo Using environment with ...
echo "DAMASK $DAMASK_ROOT"
echo "DAMASK $DAMASK_ROOT $BRANCH"
echo "Spectral Solver $SOLVER"
echo "Post Processing $PROCESSING"
if [ "x$PETSC_DIR" != "x" ]; then
@ -94,7 +97,7 @@ fi
export DAMASK_NUM_THREADS
export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH
for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN; do
for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN BRANCH; do
unset "${var}"
done
for var in DAMASK MSC; do

19
env/DAMASK.zsh vendored
View File

@ -21,16 +21,19 @@ set() {
source $DAMASK_ROOT/CONFIG
unset -f set
# add BRANCH if DAMASK_ROOT is a git repository
cd $DAMASK_ROOT >/dev/null; BRANCH=$(git branch 2>/dev/null| grep -E '^\* '); cd - >/dev/null
# add DAMASK_BIN if present
[ "x$DAMASK_BIN != x" ] && PATH=$DAMASK_BIN:$PATH
[[ "x$DAMASK_BIN" != "x" ]] && PATH=$DAMASK_BIN:$PATH
SOLVER=$(which DAMASK_spectral || true 2>/dev/null)
[ "x$SOLVER" = "x" ] && SOLVER=$(blink 'Not found!')
[[ "x$SOLVER" == "x" ]] && SOLVER=$(blink 'Not found!')
PROCESSING=$(which postResults || true 2>/dev/null)
[ "x$PROCESSING" = "x" ] && PROCESSING=$(blink 'Not found!')
[[ "x$PROCESSING" == "x" ]] && PROCESSING=$(blink 'Not found!')
[ "x$DAMASK_NUM_THREADS" = "x" ] && DAMASK_NUM_THREADS=1
[[ "x$DAMASK_NUM_THREADS" == "x" ]] && DAMASK_NUM_THREADS=1
# currently, there is no information that unlimited causes problems
# still, http://software.intel.com/en-us/forums/topic/501500 suggest to fix it
@ -50,16 +53,16 @@ if [ ! -z "$PS1" ]; then
echo https://damask.mpie.de
echo
echo "Using environment with ..."
echo "DAMASK $DAMASK_ROOT"
echo "DAMASK $DAMASK_ROOT $BRANCH"
echo "Spectral Solver $SOLVER"
echo "Post Processing $PROCESSING"
if [ "x$PETSC_DIR" != "x" ]; then
if [ "x$PETSC_DIR" != "x" ]; then
echo -n "PETSc location "
[ -d $PETSC_DIR ] && echo $PETSC_DIR || blink $PETSC_DIR
[[ $(canonicalPath "$PETSC_DIR") == $PETSC_DIR ]] \
|| echo " ~~> "$(canonicalPath "$PETSC_DIR")
fi
[[ "x$PETSC_ARCH" == "x" ]] \
[[ "x$PETSC_ARCH" == "x" ]] \
|| echo "PETSc architecture $PETSC_ARCH"
echo -n "MSC.Marc/Mentat "
[ -d $MSC_ROOT ] && echo $MSC_ROOT || blink $MSC_ROOT
@ -87,7 +90,7 @@ fi
export DAMASK_NUM_THREADS
export PYTHONPATH=$DAMASK_ROOT/lib:$PYTHONPATH
for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN; do
for var in BASE STAT SOLVER PROCESSING FREE DAMASK_BIN BRANCH; do
unset "${var}"
done
for var in DAMASK MSC; do

1
examples/ConfigFiles/Phase_Phenopowerlaw_Aluminum.config
View File

@ -6,7 +6,6 @@ plasticity phenopowerlaw
(output) shearrate_slip
(output) resolvedstress_slip
(output) accumulated_shear_slip
(output) totalshear
lattice_structure fcc
Nslip 12 # per family

1
examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Ferrite.config
View File

@ -19,4 +19,3 @@ tausat_slip 222.e6 412.7e6 # per family, optimization long
h0_slipslip 1000.0e6
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
w0_slip 2.0
(output) totalshear

1
examples/ConfigFiles/Phase_Phenopowerlaw_BCC-Martensite.config
View File

@ -19,4 +19,3 @@ tausat_slip 872.9e6 971.2e6 # per family
h0_slipslip 563.0e9
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
a_slip 2.0
(output) totalshear

2
examples/ConfigFiles/Phase_Phenopowerlaw_Gold.config
View File

@ -14,11 +14,9 @@ plasticity phenopowerlaw
(output) resistance_slip
(output) shearrate_slip
(output) resolvedstress_slip
(output) totalshear
(output) resistance_twin
(output) shearrate_twin
(output) resolvedstress_twin
(output) totalvolfrac_twin
lattice_structure fcc
Nslip 12 # per family

2
examples/ConfigFiles/Phase_Phenopowerlaw_Magnesium.config
View File

@ -9,11 +9,9 @@ elasticity hooke
(output) resistance_slip
(output) shearrate_slip
(output) resolvedstress_slip
(output) totalshear
(output) resistance_twin
(output) shearrate_twin
(output) resolvedstress_twin
(output) totalvolfrac_twin
lattice_structure hex
covera_ratio 1.62350 # from Tromans 2011, Elastic Anisotropy of HCP Metal Crystals and Polycrystals

2
examples/ConfigFiles/Phase_Phenopowerlaw_cpTi-alpha.config
View File

@ -5,11 +5,9 @@ elasticity hooke
# (output) resistance_slip
# (output) shearrate_slip
# (output) resolvedstress_slip
# (output) totalshear
# (output) resistance_twin
# (output) shearrate_twin
# (output) resolvedstress_twin
# (output) totalvolfrac_twin
lattice_structure hex
covera_ratio 1.587

2
examples/ConfigFiles/Phase_Phenopowerlaw_multiField.config
View File

@ -6,12 +6,10 @@ plasticity phenopowerlaw
(output) shearrate_slip
(output) resolvedstress_slip
(output) accumulated_shear_slip
(output) totalshear
(output) resistance_twin
(output) shearrate_twin
(output) resolvedstress_twin
(output) accumulated_shear_twin
(output) totalvolfrac_twin
lattice_structure fcc
Nslip 12 # per family

2
lib/damask/__init__.py
View File

@ -13,7 +13,7 @@ from .asciitable import ASCIItable # noqa
from .config import Material # noqa
from .colormaps import Colormap, Color # noqa
from .orientation import Quaternion, Rodrigues, Symmetry, Orientation # noqa
from .orientation import Quaternion, Symmetry, Orientation # noqa
#from .block import Block # only one class
from .result import Result # noqa

4
lib/damask/asciitable.py
View File

@ -493,8 +493,8 @@ class ASCIItable():
(d if str(c) != str(labels[present[i]]) else
1)))
use = np.array(columns) if len(columns) > 0 else None
self.tags = list(np.array(self.tags)[use]) # update labels with valid subset
self.tags = list(np.array(self.__IO__['tags'])[use]) # update labels with valid subset
self.data = np.loadtxt(self.__IO__['in'],usecols=use,ndmin=2)
# self.data = np.genfromtxt(self.__IO__['in'],dtype=None,names=self.tags,usecols=use)

450
lib/damask/orientation.py
View File

@ -7,24 +7,6 @@
import math,os
import numpy as np
# ******************************************************************************************
class Rodrigues:
def __init__(self, vector = np.zeros(3)):
self.vector = vector
def asQuaternion(self):
norm = np.linalg.norm(self.vector)
halfAngle = np.arctan(norm)
return Quaternion(np.cos(halfAngle),np.sin(halfAngle)*self.vector/norm)
def asAngleAxis(self):
norm = np.linalg.norm(self.vector)
halfAngle = np.arctan(norm)
return (2.0*halfAngle,self.vector/norm)
# ******************************************************************************************
class Quaternion:
u"""
@ -33,7 +15,7 @@ class Quaternion:
All methods and naming conventions based on Rowenhorst_etal2015
Convention 1: coordinate frames are right-handed
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 4: Euler angle triplets are implemented using the Bunge convention,
with the angular ranges as [0, 2π],[0, π],[0, 2π]
@ -48,253 +30,167 @@ class Quaternion:
"""
def __init__(self,
quatArray = [1.0,0.0,0.0,0.0]):
"""Initializes to identity if not given"""
self.w, \
self.x, \
self.y, \
self.z = quatArray
quat = None,
q = 1.0,
p = np.zeros(3,dtype=float)):
"""Initializes to identity unless specified"""
self.q = quat[0] if quat is not None else q
self.p = np.array(quat[1:4]) if quat is not None else p
self.homomorph()
def __iter__(self):
"""Components"""
return iter([self.w,self.x,self.y,self.z])
return iter(self.asList())
def __copy__(self):
"""Create copy"""
Q = Quaternion([self.w,self.x,self.y,self.z])
return Q
"""Copy"""
return self.__class__(q=self.q,p=self.p.copy())
copy = __copy__
def __repr__(self):
"""Readbable string"""
return 'Quaternion(real=%+.6f, imag=<%+.6f, %+.6f, %+.6f>)' % \
(self.w, self.x, self.y, self.z)
"""Readable string"""
return 'Quaternion(real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p)
def __pow__(self, exponent):
"""Power"""
omega = math.acos(self.w)
vRescale = math.sin(exponent*omega)/math.sin(omega)
Q = Quaternion()
Q.w = math.cos(exponent*omega)
Q.x = self.x * vRescale
Q.y = self.y * vRescale
Q.z = self.z * vRescale
return Q
omega = math.acos(self.q)
return self.__class__(q= math.cos(exponent*omega),
p=self.p * math.sin(exponent*omega)/math.sin(omega))
def __ipow__(self, exponent):
"""In-place power"""
omega = math.acos(self.w)
vRescale = math.sin(exponent*omega)/math.sin(omega)
self.w = np.cos(exponent*omega)
self.x *= vRescale
self.y *= vRescale
self.z *= vRescale
omega = math.acos(self.q)
self.q = math.cos(exponent*omega)
self.p *= math.sin(exponent*omega)/math.sin(omega)
return self
def __mul__(self, other):
"""Multiplication"""
try: # quaternion
Aw = self.w
Ax = self.x
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
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
try: # quaternion
return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p))
except: pass
try: # vector (perform active rotation, i.e. q*v*q.conjugated)
w = self.w
x = self.x
y = self.y
z = self.z
Vx = other[0]
Vy = other[1]
Vz = other[2]
try: # vector (perform passive rotation)
( x, y, z) = self.p
(Vx,Vy,Vz) = other[0:3]
A = self.q*self.q - np.dot(self.p,self.p)
B = 2.0 * (x*Vx + y*Vy + z*Vz)
C = 2.0 * P*self.q
return np.array([\
w * w * Vx + 2 * y * w * Vz - 2 * z * w * Vy + \
x * x * Vx + 2 * y * x * Vy + 2 * z * x * Vz - \
z * z * Vx - y * 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 ])
return np.array([
A*Vx + B*x + C*(y*Vz - z*Vy),
A*Vy + B*y + C*(z*Vx - x*Vz),
A*Vz + B*z + C*(x*Vy - y*Vx),
])
except: pass
try: # scalar
Q = self.copy()
Q.w *= other
Q.x *= other
Q.y *= other
Q.z *= other
return Q
try: # scalar
return self.__class__(q=self.q*other,
p=self.p*other)
except:
return self.copy()
def __imul__(self, other):
"""In-place multiplication"""
try: # Quaternion
Aw = self.w
Ax = self.x
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
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
try: # Quaternion
self.q = self.q*other.q - np.dot(self.p,other.p)
self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
except: pass
return self
def __div__(self, other):
"""Division"""
if isinstance(other, (int,float)):
w = self.w / other
x = self.x / other
y = self.y / other
z = self.z / other
return self.__class__([w,x,y,z])
return self.__class__(q=self.q / other,
p=self.p / other)
else:
return NotImplemented
def __idiv__(self, other):
"""In-place division"""
if isinstance(other, (int,float)):
self.w /= other
self.x /= other
self.y /= other
self.z /= other
self.q /= other
self.p /= other
return self
def __add__(self, other):
"""Addition"""
if isinstance(other, Quaternion):
w = self.w + other.w
x = self.x + other.x
y = self.y + other.y
z = self.z + other.z
return self.__class__([w,x,y,z])
return self.__class__(q=self.q + other.q,
p=self.p + other.p)
else:
return NotImplemented
def __iadd__(self, other):
"""In-place addition"""
if isinstance(other, Quaternion):
self.w += other.w
self.x += other.x
self.y += other.y
self.z += other.z
self.q += other.q
self.p += other.p
return self
def __sub__(self, other):
"""Subtraction"""
if isinstance(other, Quaternion):
Q = self.copy()
Q.w -= other.w
Q.x -= other.x
Q.y -= other.y
Q.z -= other.z
return Q
return self.__class__(q=self.q - other.q,
p=self.p - other.p)
else:
return self.copy()
return NotImplemented
def __isub__(self, other):
"""In-place subtraction"""
if isinstance(other, Quaternion):
self.w -= other.w
self.x -= other.x
self.y -= other.y
self.z -= other.z
self.q -= other.q
self.p -= other.p
return self
def __neg__(self):
"""Additive inverse"""
self.w = -self.w
self.x = -self.x
self.y = -self.y
self.z = -self.z
self.q = -self.q
self.p = -self.p
return self
def __abs__(self):
"""Norm"""
return math.sqrt(self.w ** 2 + \
self.x ** 2 + \
self.y ** 2 + \
self.z ** 2)
return math.sqrt(self.q ** 2 + np.dot(self.p,self.p))
magnitude = __abs__
def __eq__(self,other):
"""Equal at e-8 precision"""
return (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) \
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)
"""Equal (sufficiently close) to each other"""
return np.isclose(( self-other).magnitude(),0.0) \
or np.isclose((-self-other).magnitude(),0.0)
def __ne__(self,other):
"""Not equal at e-8 precision"""
return not self.__eq__(self,other)
"""Not equal (sufficiently close) to each other"""
return not self.__eq__(other)
def __cmp__(self,other):
"""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):
return self.w ** 2 + \
self.x ** 2 + \
self.y ** 2 + \
self.z ** 2
def identity(self):
self.w = 1.
self.x = 0.
self.y = 0.
self.z = 0.
return self
return self.q ** 2 + np.dot(self.p,self.p)
def normalize(self):
d = self.magnitude()
if d > 0.0:
self /= d
self.q /= d
self.p /= d
return self
def conjugate(self):
self.x = -self.x
self.y = -self.y
self.z = -self.z
return self
def inverse(self):
d = self.magnitude()
if d > 0.0:
self.conjugate()
self /= d
self.p = -self.p
return self
def homomorph(self):
if self.w < 0.0:
self.w = -self.w
self.x = -self.x
self.y = -self.y
self.z = -self.z
if self.q < 0.0:
self.q = -self.q
self.p = -self.p
return self
def normalized(self):
@ -303,63 +199,73 @@ class Quaternion:
def conjugated(self):
return self.copy().conjugate()
def inversed(self):
return self.copy().inverse()
def homomorphed(self):
return self.copy().homomorph()
def asList(self):
return [i for i in self]
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 [self.q]+list(self.p)
def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.)
return np.outer(self.asList(),self.asList())
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(
[[ qbarhalf + self.x**2 , self.x*self.y - self.w*self.z, self.x*self.z + self.w*self.y],
[ self.x*self.y + self.w*self.z, qbarhalf + self.y**2 , self.y*self.z - self.w*self.x],
[ self.x*self.z - self.w*self.y, self.y*self.z + self.w*self.x, qbarhalf + self.z**2 ],
[[ qbarhalf + self.p[0]**2 ,
self.p[0]*self.p[1] -P* self.q*self.p[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,
degrees = False):
if self.w > 1:
self.normalize()
degrees = False,
flat = False):
s = math.sqrt(1. - self.w**2)
x = 2*self.w**2 - 1.
y = 2*self.w * s
angle = 2.0*math.acos(self.q)
angle = math.atan2(y,x)
if angle < 0.0:
angle *= -1.
s *= -1.
if np.isclose(angle,0.0):
angle = 0.0
axis = np.array([0.0,0.0,1.0])
elif np.isclose(self.q,0.0):
angle = math.pi
axis = self.p
else:
axis = np.sign(self.q)*self.p/np.linalg.norm(self.p)
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]))
angle = np.degrees(angle) if degrees else angle
return np.hstack((angle,axis)) if flat else (angle,axis)
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 np.isclose(self.q,0.0) else self.p/self.q
def asEulers(self,
degrees = False):
"""Orientation as Bunge-Euler angles."""
q03 = self.w**2+self.z**2
q12 = self.x**2+self.y**2
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
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)
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])
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])
if np.isclose(chi,0.0) and np.isclose(q12,0.0):
eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0])
elif np.isclose(chi,0.0) and np.isclose(q03,0.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:
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((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),
])
eulers %= 2.0*math.pi # enforce positive angles
return np.degrees(eulers) if degrees else eulers
@ -371,25 +277,28 @@ class Quaternion:
@classmethod
def fromRandom(cls,randomSeed = None):
import binascii
if randomSeed is None:
randomSeed = int(os.urandom(4).encode('hex'), 16)
randomSeed = int(binascii.hexlify(os.urandom(4)),16)
np.random.seed(randomSeed)
r = np.random.random(3)
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])
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])
return cls([w,x,y,z])
A = math.sqrt(max(0.0,r[2]))
B = math.sqrt(max(0.0,1.0-r[2]))
w = math.cos(2.0*math.pi*r[0])*A
x = math.sin(2.0*math.pi*r[1])*B
y = math.cos(2.0*math.pi*r[1])*B
z = math.sin(2.0*math.pi*r[0])*A
return cls(quat=[w,x,y,z])
@classmethod
def fromRodrigues(cls, 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)
w = c
x,y,z = rodrigues/c
return cls([w,x,y,z])
return cls(q=c,p=s*rodrigues/norm)
@classmethod
@ -397,22 +306,19 @@ class Quaternion:
angle,
axis,
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)
angle = np.radians(angle) if degrees else angle
s = math.sin(0.5 * angle)
w = math.cos(0.5 * angle)
x = axis[0] * s
y = axis[1] * s
z = axis[2] * s
return cls([w,x,y,z])
c = math.cos(0.5 * angle)
return cls(q=c,p=axis*s)
@classmethod
def fromEulers(cls,
eulers,
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
sigma = 0.5*(eulers[0]+eulers[2])
@ -420,11 +326,13 @@ class Quaternion:
c = np.cos(0.5*eulers[1])
s = np.sin(0.5*eulers[1])
w = c * np.cos(sigma)
x = -s * np.cos(delta)
y = -s * np.sin(delta)
z = -c * np.sin(sigma)
return cls([w,x,y,z])
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
w = c * np.cos(sigma)
x = -P * s * np.cos(delta)
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,
@ -435,16 +343,18 @@ class Quaternion:
if m.shape != (3,3) and np.prod(m.shape) == 9:
m = m.reshape(3,3)
w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2])
x = 0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
y = 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])
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
w = 0.5*math.sqrt(max(0.0,1.0+m[0,0]+m[1,1]+m[2,2]))
x = P*0.5*math.sqrt(max(0.0,1.0+m[0,0]-m[1,1]-m[2,2]))
y = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]+m[1,1]-m[2,2]))
z = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]-m[1,1]+m[2,2]))
x *= -1 if m[2,1] < m[1,2] else 1
y *= -1 if m[0,2] < m[2,0] 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
@ -458,36 +368,30 @@ class Quaternion:
assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion)
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.:
costheta = -costheta
q1 = q1.conjugated()
elif costheta > 1:
costheta = 1
elif costheta > 1.:
costheta = 1.
theta = math.acos(costheta)
if abs(theta) < 0.01:
Q.w = q2.w
Q.x = q2.x
Q.y = q2.y
Q.z = q2.z
Q.q = q2.q
Q.p = q2.p
return Q
sintheta = math.sqrt(1.0 - costheta * costheta)
if abs(sintheta) < 0.01:
Q.w = (q1.w + q2.w) * 0.5
Q.x = (q1.x + q2.x) * 0.5
Q.y = (q1.y + q2.y) * 0.5
Q.z = (q1.z + q2.z) * 0.5
Q.q = (q1.q + q2.q) * 0.5
Q.p = (q1.p + q2.p) * 0.5
return Q
ratio1 = math.sin((1 - t) * theta) / sintheta
ratio2 = math.sin(t * theta) / sintheta
ratio1 = math.sin((1.0 - t) * theta) / sintheta
ratio2 = math.sin( t * theta) / sintheta
Q.w = q1.w * ratio1 + q2.w * ratio2
Q.x = q1.x * ratio1 + q2.x * ratio2
Q.y = q1.y * ratio1 + q2.y * ratio2
Q.z = q1.z * ratio1 + q2.z * ratio2
Q.q = q1.q * ratio1 + q2.q * ratio2
Q.p = q1.p * ratio1 + q2.p * ratio2
return Q
@ -512,21 +416,21 @@ class Symmetry:
def __repr__(self):
"""Readbable string"""
return '%s' % (self.lattice)
"""Readable string"""
return '{}'.format(self.lattice)
def __eq__(self, other):
"""Equal"""
"""Equal to other"""
return self.lattice == other.lattice
def __neq__(self, other):
"""Not equal"""
"""Not equal to other"""
return not self.__eq__(other)
def __cmp__(self,other):
"""Linear ordering"""
myOrder = Symmetry.lattices.index(self.lattice)
myOrder = Symmetry.lattices.index(self.lattice)
otherOrder = Symmetry.lattices.index(other.lattice)
return (myOrder > otherOrder) - (myOrder < otherOrder)
@ -598,7 +502,7 @@ class Symmetry:
]
return list(map(Quaternion,
np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))]))
np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))]))
def equivalentQuaternions(self,
@ -610,7 +514,7 @@ class Symmetry:
def inFZ(self,R):
"""Check whether given Rodrigues vector falls into fundamental zone of own symmetry."""
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentially passed quaternion
if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentally passed quaternion
# fundamental zone in Rodrigues space is point symmetric around origin
R = abs(R)
if self.lattice == 'cubic':
@ -722,7 +626,7 @@ class Symmetry:
else:
return True
v = np.array(vector,dtype = float)
v = np.array(vector,dtype=float)
if proper: # check both improper ...
theComponents = np.dot(basis['improper'],v)
inSST = np.all(theComponents >= 0.0)
@ -737,10 +641,10 @@ class Symmetry:
if color: # have to return color array
if inSST:
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
else:
rgb = np.zeros(3,'d')
rgb = np.zeros(3,dtype=float)
return (inSST,rgb)
else:
return inSST
@ -780,8 +684,9 @@ class Orientation:
self.quaternion = Quaternion.fromRodrigues(Rodrigues)
elif isinstance(quaternion, Quaternion): # based on given quaternion
self.quaternion = quaternion.homomorphed()
elif isinstance(quaternion, np.ndarray) and quaternion.shape == (4,): # based on given quaternion-like array
self.quaternion = Quaternion(quaternion).homomorphed()
elif (isinstance(quaternion, np.ndarray) and quaternion.shape == (4,)) or \
(isinstance(quaternion, list) and len(quaternion) == 4 ): # based on given quaternion-like array
self.quaternion = Quaternion(quat=quaternion).homomorphed()
self.symmetry = Symmetry(symmetry)
@ -794,10 +699,12 @@ class Orientation:
def __repr__(self):
"""Value as all implemented representations"""
return 'Symmetry: %s\n' % (self.symmetry) + \
'Quaternion: %s\n' % (self.quaternion) + \
'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \
'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers(degrees=True))) )
return '\n'.join([
'Symmetry: {}'.format(self.symmetry),
'Quaternion: {}'.format(self.quaternion),
'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):
return self.quaternion.asList()
@ -813,8 +720,9 @@ class Orientation:
rodrigues = property(asRodrigues)
def asAngleAxis(self,
degrees = False):
return self.quaternion.asAngleAxis(degrees)
degrees = False,
flat = False):
return self.quaternion.asAngleAxis(degrees,flat)
angleAxis = property(asAngleAxis)
def asMatrix(self):
@ -927,7 +835,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
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)

24
lib/damask/util.py
View File

@ -132,6 +132,30 @@ class extendableOption(Option):
else:
Option.take_action(self, action, dest, opt, value, values, parser)
# Print iterations progress
# from https://gist.github.com/aubricus/f91fb55dc6ba5557fbab06119420dd6a
def progressBar(iteration, total, prefix='', suffix='', decimals=1, bar_length=100):
"""
Call in a loop to create terminal progress bar
@params:
iteration - Required : current iteration (Int)
total - Required : total iterations (Int)
prefix - Optional : prefix string (Str)
suffix - Optional : suffix string (Str)
decimals - Optional : positive number of decimals in percent complete (Int)
bar_length - Optional : character length of bar (Int)
"""
str_format = "{0:." + str(decimals) + "f}"
percents = str_format.format(100 * (iteration / float(total)))
filled_length = int(round(bar_length * iteration / float(total)))
bar = '█' * filled_length + '-' * (bar_length - filled_length)
sys.stderr.write('\r%s |%s| %s%s %s' % (prefix, bar, percents, '%', suffix)),
if iteration == total: sys.stderr.write('\n\n')
sys.stderr.flush()
# -----------------------------
class backgroundMessage(threading.Thread):
"""Reporting with animation to indicate progress"""

89
processing/post/addAPS34IDEstrainCoords.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
@ -19,55 +19,50 @@ Transform X,Y,Z,F APS BeamLine 34 coordinates to x,y,z APS strain coordinates.
""", version = scriptID)
parser.add_option('-f','--frame', dest='frame', nargs=4, type='string', metavar='string string string string',
help='APS X,Y,Z coords, and depth F')
parser.set_defaults(frame = None)
parser.add_option('-f',
'--frame',
dest='frame',
metavar='string',
help='APS X,Y,Z coords')
parser.add_option('--depth',
dest='depth',
metavar='string',
help='depth')
(options,filenames) = parser.parse_args()
if options.frame is None:
parser.error('no data column specified...')
parser.error('frame not specified')
if options.depth is None:
parser.error('depth not specified')
# --- loop over input files ------------------------------------------------------------------------
datainfo = {'len':3,
'label':[]
}
datainfo['label'] += options.frame
# --- loop over input files -------------------------------------------------------------------------
if filenames == []:
filenames = ['STDIN']
if filenames == []: filenames = [None]
for name in filenames:
if name == 'STDIN':
file = {'name':'STDIN', 'input':sys.stdin, 'output':sys.stdout, 'croak':sys.stderr}
file['croak'].write('\033[1m'+scriptName+'\033[0m\n')
else:
if not os.path.exists(name): continue
file = {'name':name, 'input':open(name), 'output':open(name+'_tmp','w'), 'croak':sys.stderr}
file['croak'].write('\033[1m'+scriptName+'\033[0m: '+file['name']+'\n')
try: table = damask.ASCIItable(name = name,
buffered = False)
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks -----------------------------------------
errors = []
if table.label_dimension(options.frame) != 3:
errors.append('input {} does not have dimension 3.'.format(options.frame))
if table.label_dimension(options.depth) != 1:
errors.append('input {} does not have dimension 1.'.format(options.depth))
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
table = damask.ASCIItable(file['input'],file['output'],buffered=False) # make unbuffered ASCII_table
table.head_read() # read ASCII header info
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
# --------------- figure out columns to process ---------------------------------------------------
active = []
column = {}
columnMissing = False
for label in datainfo['label']:
key = label
if key in table.labels(raw = True):
active.append(label)
column[label] = table.labels.index(key) # remember columns of requested data
else:
file['croak'].write('column %s not found...\n'%label)
columnMissing = True
if columnMissing: continue
# ------------------------------------------ assemble header ---------------------------------------
table.labels_append(['%i_coord'%(i+1) for i in range(3)]) # extend ASCII header with new labels
table.head_write()
@ -77,21 +72,15 @@ for name in filenames:
RotMat2TSL=np.array([[1., 0., 0.],
[0., np.cos(theta), np.sin(theta)], # Orientation to account for -135 deg
[0., -np.sin(theta), np.cos(theta)]]) # rotation for TSL convention
vec = np.zeros(4)
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
for i,label in enumerate(active):
vec[i] = table.data[column[label]]
coord = list(map(float,table.data[table.label_index(options.frame):table.label_index(options.frame)+3]))
depth = float(table.data[table.label_index(options.depth)])
table.data_append(np.dot(RotMat2TSL,np.array([-vec[0], -vec[1],-vec[2]+vec[3]])))
table.data_append(np.dot(RotMat2TSL,np.array([-coord[0],-coord[1],-coord[2]+depth])))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output result -----------------------------------------
outputAlive and table.output_flush() # just in case of buffered ASCII table
# ------------------------------------------ output finalization -----------------------------------
table.input_close() # close input ASCII table (works for stdin)
table.output_close() # close output ASCII table (works for stdout)
if file['name'] != 'STDIN':
os.rename(file['name']+'_tmp',file['name']) # overwrite old one with tmp new
table.close() # close ASCII tables

121
processing/post/addGrainID.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,time,copy
import os,sys,copy
import numpy as np
import damask
from optparse import OptionParser
@ -29,49 +29,28 @@ parser.add_option('-d',
parser.add_option('-s',
'--symmetry',
dest = 'symmetry',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'crystal symmetry [%default]')
parser.add_option('-e',
'--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'label of Euler angles')
parser.add_option('--degrees',
dest = 'degrees',
action = 'store_true',
help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m',
'--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'label of orientation matrix')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'label of crystal frame a vector')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'label of crystal frame b vector')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'label of crystal frame c vector')
parser.add_option('-q',
'--quaternion',
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
type = 'string', metavar = 'string',
help = 'label of quaternion')
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.add_option('-p',
'--pos', '--position',
dest = 'pos',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('--quiet',
dest='verbose',
action = 'store_false',
help = 'hide status bar (useful when piping to file)')
parser.set_defaults(disorientation = 5,
verbose = True,
quaternion = 'orientation',
symmetry = 'cubic',
pos = 'pos',
degrees = False,
)
(options, filenames) = parser.parse_args()
@ -79,22 +58,6 @@ parser.set_defaults(disorientation = 5,
if options.radius is None:
parser.error('no radius specified.')
input = [options.eulers is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
]
if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
cos_disorientation = np.cos(np.radians(options.disorientation/2.)) # cos of half the disorientation angle
# --- loop over input files -------------------------------------------------------------------------
@ -118,9 +81,9 @@ for name in filenames:
if not 3 >= table.label_dimension(options.pos) >= 1:
errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
if not np.all(table.label_dimension(label) == dim):
errors.append('input "{}" does not have dimension {}.'.format(label,dim))
else: column = table.label_index(label)
if not np.all(table.label_dimension(options.quaternion) == 4):
errors.append('input "{}" does not have dimension 4.'.format(options.quaternion))
else: column = table.label_index(options.quaternion)
if remarks != []: damask.util.croak(remarks)
if errors != []:
@ -131,34 +94,18 @@ for name in filenames:
# ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append('grainID_{}@{:g}'.format('+'.join(label)
if isinstance(label, (list,tuple))
else label,
options.disorientation)) # report orientation source and disorientation
table.labels_append('grainID_{}@{:g}'.format(options.quaternion,options.disorientation)) # report orientation source and disorientation
table.head_write()
# ------------------------------------------ process data ------------------------------------------
# ------------------------------------------ build KD tree -----------------------------------------
# --- start background messaging
bg = damask.util.backgroundMessage()
bg.start()
bg.set_message('reading positions...')
table.data_readArray(options.pos) # read position vectors
grainID = -np.ones(len(table.data),dtype=int)
start = tick = time.clock()
bg.set_message('building KD tree...')
Npoints = table.data.shape[0]
kdtree = spatial.KDTree(copy.deepcopy(table.data))
# ------------------------------------------ assign grain IDs --------------------------------------
tick = time.clock()
orientations = [] # quaternions found for grain
memberCounts = [] # number of voxels in grain
p = 0 # point counter
@ -169,26 +116,11 @@ for name in filenames:
table.data_rewind()
while table.data_read(): # read next data line of ASCII table
if p > 0 and p % 1000 == 0:
if options.verbose and Npoints > 100 and p%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=p,total=Npoints)
time_delta = (time.clock()-tick) * (len(grainID) - p) / p
bg.set_message('(%02i:%02i:%02i) processing point %i of %i (grain count %i)...'\
%(time_delta//3600,time_delta%3600//60,time_delta%60,p,len(grainID),np.count_nonzero(memberCounts)))
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(map(float,table.data[column:column+3]))*toRadians,
symmetry = options.symmetry).reduced()
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(map(float,table.data[column:column+9])).reshape(3,3).transpose(),
symmetry = options.symmetry).reduced()
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3])).reshape(3,3),
symmetry = options.symmetry).reduced()
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(map(float,table.data[column:column+4])),
symmetry = options.symmetry).reduced()
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
matched = False
alreadyChecked = {}
@ -200,9 +132,9 @@ for name in filenames:
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
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
if disorientation.quaternion.w >= bestDisorientation.w: # ... and better than current best?
if disorientation.quaternion.q >= bestDisorientation.q: # ... and better than current best?
matched = True
matchedID = gID # remember that grain
bestDisorientation = disorientation.quaternion
@ -233,13 +165,12 @@ for name in filenames:
outputAlive = True
p = 0
damask.util.progressBar(iteration=1,total=1)
while outputAlive and table.data_read(): # read next data line of ASCII table
table.data_append(1+packingMap[grainID[p]]) # add (condensed) grain ID
outputAlive = table.data_write() # output processed line
p += 1
bg.set_message('done after {} seconds'.format(time.clock()-start))
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables

79
processing/post/addIPFcolor.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import os,sys
import numpy as np
from optparse import OptionParser
import damask
@ -18,66 +18,29 @@ Add RGB color value corresponding to TSL-OIM scheme for inverse pole figures.
""", version = scriptID)
parser.add_option('-p', '--pole',
parser.add_option('-p',
'--pole',
dest = 'pole',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'lab frame direction for inverse pole figure [%default]')
parser.add_option('-s', '--symmetry',
parser.add_option('-s',
'--symmetry',
dest = 'symmetry',
type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string',
help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
parser.add_option('-e', '--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'Euler angles label')
parser.add_option('-d', '--degrees',
dest = 'degrees',
action = 'store_true',
help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m', '--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'orientation matrix label')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'crystal frame c vector label')
parser.add_option('-q', '--quaternion',
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
type = 'string', metavar = 'string',
help = 'quaternion label')
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.set_defaults(pole = (0.0,0.0,1.0),
quaternion = 'orientation',
symmetry = damask.Symmetry.lattices[-1],
degrees = False,
)
(options, filenames) = parser.parse_args()
input = [options.eulers is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
]
if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
pole = np.array(options.pole)
pole /= np.linalg.norm(pole)
@ -98,12 +61,12 @@ for name in filenames:
# ------------------------------------------ sanity checks ----------------------------------------
if not np.all(table.label_dimension(label) == dim):
damask.util.croak('input {} does not have dimension {}.'.format(label,dim))
if not table.label_dimension(options.quaternion) == 4:
damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion))
table.close(dismiss = True) # close ASCIItable and remove empty file
continue
column = table.label_index(label)
column = table.label_index(options.quaternion)
# ------------------------------------------ assemble header ---------------------------------------
@ -115,20 +78,8 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,
symmetry = options.symmetry).reduced()
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),
symmetry = options.symmetry).reduced()
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3),
symmetry = options.symmetry).reduced()
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
table.data_append(o.IPFcolor(pole))
outputAlive = table.data_write() # output processed line

2
processing/post/addInfo.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os

2
processing/post/addMises.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math

143
processing/post/addOrientations.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import os,sys
import numpy as np
from optparse import OptionParser
import damask
@ -9,6 +9,31 @@ import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# convention conformity checks
# --------------------------------------------------------------------
def check_Eulers(eulers):
if np.any(eulers < 0.0) or np.any(eulers > 2.0*np.pi) or eulers[1] > np.pi: # Euler angles within valid range?
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eulers))
return eulers
def check_quaternion(q):
if q[0] < 0.0: # positive first quaternion component?
raise ValueError('quaternion has negative first component.\n{}'.format(q[0]))
if not np.isclose(np.linalg.norm(q), 1.0): # unit quaternion?
raise ValueError('quaternion is not of unit length.\n{} {} {} {}'.format(*q))
return q
def check_matrix(M):
if not np.isclose(np.linalg.det(M),1.0): # proper rotation?
raise ValueError('matrix is not a proper rotation.\n{}'.format(M))
if not np.isclose(np.dot(M[0],M[1]), 0.0) \
or not np.isclose(np.dot(M[1],M[2]), 0.0) \
or not np.isclose(np.dot(M[2],M[0]), 0.0): # all orthogonal?
raise ValueError('matrix is not orthogonal.\n{}'.format(M))
return M
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
@ -21,58 +46,64 @@ Additional (globally fixed) rotations of the lab frame and/or crystal frame can
""", version = scriptID)
outputChoices = ['quaternion','rodrigues','eulers']
parser.add_option('-o', '--output',
outputChoices = {
'quaternion': ['quat',4],
'rodrigues': ['rodr',3],
'eulers': ['eulr',3],
'matrix': ['mtrx',9],
'angleaxis': ['aaxs',4],
}
parser.add_option('-o',
'--output',
dest = 'output',
action = 'extend', metavar = '<string LIST>',
help = 'output orientation formats {{{}}}'.format(', '.join(outputChoices)))
parser.add_option('-s', '--symmetry',
dest = 'symmetry',
type = 'choice', choices = damask.Symmetry.lattices[1:], metavar='string',
help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
parser.add_option('-d', '--degrees',
parser.add_option('-d',
'--degrees',
dest = 'degrees',
action = 'store_true',
help = 'angles are given in degrees [%default]')
parser.add_option('-R', '--labrotation',
help = 'all angles in degrees')
parser.add_option('-R',
'--labrotation',
dest='labrotation',
type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
help = 'angle and axis of additional lab frame rotation')
parser.add_option('-r', '--crystalrotation',
parser.add_option('-r',
'--crystalrotation',
dest='crystalrotation',
type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
help = 'angle and axis of additional crystal frame rotation')
parser.add_option( '--eulers',
parser.add_option('--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'Euler angles label')
parser.add_option( '--rodrigues',
parser.add_option('--rodrigues',
dest = 'rodrigues',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'Rodrigues vector label')
parser.add_option( '--matrix',
parser.add_option('--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'orientation matrix label')
parser.add_option( '--quaternion',
parser.add_option('--quaternion',
dest = 'quaternion',
type = 'string', metavar = 'string',
metavar = 'string',
help = 'quaternion label')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'crystal frame c vector label')
parser.add_option('-x',
dest = 'x',
metavar = 'string',
help = 'label of lab x vector (expressed in crystal coords)')
parser.add_option('-y',
dest = 'y',
metavar = 'string',
help = 'label of lab y vector (expressed in crystal coords)')
parser.add_option('-z',
dest = 'z',
metavar = 'string',
help = 'label of lab z vector (expressed in crystal coords)')
parser.set_defaults(output = [],
symmetry = damask.Symmetry.lattices[-1],
labrotation = (0.,1.,1.,1.), # no rotation about 1,1,1
crystalrotation = (0.,1.,1.,1.), # no rotation about 1,1,1
degrees = False,
@ -86,9 +117,9 @@ if options.output == [] or (not set(options.output).issubset(set(outputChoices))
input = [options.eulers is not None,
options.rodrigues is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.x is not None and \
options.y is not None and \
options.z is not None,
options.matrix is not None,
options.quaternion is not None,
]
@ -97,13 +128,14 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
(options.rodrigues,3,'rodrigues'),
([options.a,options.b,options.c],[3,3,3],'frame'),
([options.x,options.y,options.z],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
r = damask.Quaternion().fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation
R = damask.Quaternion().fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation
toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
r = damask.Quaternion.fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation
R = damask.Quaternion.fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation
# --- loop over input files ------------------------------------------------------------------------
@ -137,32 +169,31 @@ for name in filenames:
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for output in options.output:
if output == 'quaternion': table.labels_append(['{}_{}_{}({})'.format(i+1,'quat',options.symmetry,label) for i in range(4)])
elif output == 'rodrigues': table.labels_append(['{}_{}_{}({})'.format(i+1,'rodr',options.symmetry,label) for i in range(3)])
elif output == 'eulers': table.labels_append(['{}_{}_{}({})'.format(i+1,'eulr',options.symmetry,label) for i in range(3)])
if output in outputChoices:
table.labels_append(['{}_{}({})'.format(i+1,outputChoices[output][0],label) \
for i in range(outputChoices[output][1])])
table.head_write()
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,
symmetry = options.symmetry).reduced()
if inputtype == 'eulers':
o = damask.Orientation(Eulers = check_Eulers(np.array(list(map(float,table.data[column:column+3])))*toRadians))
elif inputtype == 'rodrigues':
o = damask.Orientation(Rodrigues= np.array(list(map(float,table.data[column:column+3]))),
symmetry = options.symmetry).reduced()
o = damask.Orientation(Rodrigues = np.array(list(map(float,table.data[column:column+3]))))
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),
symmetry = options.symmetry).reduced()
o = damask.Orientation(matrix = check_matrix(np.array(list(map(float,table.data[column:column+9]))).reshape(3,3)))
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3),
symmetry = options.symmetry).reduced()
M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3).T
o = damask.Orientation(matrix = check_matrix(M/np.linalg.norm(M,axis=0)))
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
o = damask.Orientation(quaternion = check_quaternion(np.array(list(map(float,table.data[column:column+4])))))
o.quaternion = r*o.quaternion*R # apply additional lab and crystal frame rotations
@ -170,6 +201,8 @@ for name in filenames:
if output == 'quaternion': table.data_append(o.asQuaternion())
elif output == 'rodrigues': table.data_append(o.asRodrigues())
elif output == 'eulers': table.data_append(o.asEulers(degrees=options.degrees))
elif output == 'matrix': table.data_append(o.asMatrix())
elif output == 'angleaxis': table.data_append(o.asAngleAxis(degrees=options.degrees,flat=True))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------

86
processing/post/addPole.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
import os,sys
import numpy as np
from optparse import OptionParser
import damask
@ -14,70 +14,32 @@ scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description = """
Add x,y coordinates of stereographic projection of given direction (pole) in crystal frame.
Add coordinates of stereographic projection of given direction (pole) in crystal frame.
""", version = scriptID)
parser.add_option('-p', '--pole',
parser.add_option('-p',
'--pole',
dest = 'pole',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'crystal frame direction for pole figure [%default]')
parser.add_option('--polar',
dest = 'polar',
action = 'store_true',
help = 'output polar coordinates r,phi [%default]')
parser.add_option('-e', '--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'Euler angles label')
parser.add_option('-d', '--degrees',
dest = 'degrees',
action = 'store_true',
help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m', '--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'orientation matrix label')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'crystal frame c vector label')
parser.add_option('-q', '--quaternion',
help = 'output polar coordinates (r,φ) instead of Cartesian coordinates (x,y)')
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
type = 'string', metavar = 'string',
help = 'quaternion label')
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.set_defaults(pole = (1.0,0.0,0.0),
degrees = False,
quaternion = 'orientation',
polar = False,
)
(options, filenames) = parser.parse_args()
input = [options.eulers is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
]
if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
pole = np.array(options.pole)
pole /= np.linalg.norm(pole)
@ -98,18 +60,13 @@ for name in filenames:
# ------------------------------------------ sanity checks ----------------------------------------
errors = []
remarks = []
if not np.all(table.label_dimension(label) == dim): errors.append('input {} does not have dimension {}.'.format(label,dim))
else: column = table.label_index(label)
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
if not table.label_dimension(options.quaternion) == 4:
damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion))
table.close(dismiss = True) # close ASCIItable and remove empty file
continue
column = table.label_index(options.quaternion)
# ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
@ -119,16 +76,7 @@ for name in filenames:
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians)
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose())
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3))
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))))
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))))
rotatedPole = o.quaternion*pole # rotate pole according to crystal orientation
(x,y) = rotatedPole[0:2]/(1.+abs(pole[2])) # stereographic projection

83
processing/post/addSchmidfactors.py
View File

@ -109,64 +109,42 @@ Add columns listing Schmid factors (and optional trace vector of selected system
""", version = scriptID)
latticeChoices = ('fcc','bcc','hex')
parser.add_option('-l','--lattice',
parser.add_option('-l',
'--lattice',
dest = 'lattice', type = 'choice', choices = latticeChoices, metavar='string',
help = 'type of lattice structure [%default] {}'.format(latticeChoices))
parser.add_option('--covera',
dest = 'CoverA', type = 'float', metavar = 'float',
help = 'C over A ratio for hexagonal systems')
parser.add_option('-f', '--force',
parser.add_option('-f',
'--force',
dest = 'force',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'force direction in lab frame [%default]')
parser.add_option('-n', '--normal',
parser.add_option('-n',
'--normal',
dest = 'normal',
type = 'float', nargs = 3, metavar = 'float float float',
help = 'stress plane normal in lab frame [%default]')
parser.add_option('-e', '--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'Euler angles label')
parser.add_option('-d', '--degrees',
dest = 'degrees',
action = 'store_true',
help = 'Euler angles are given in degrees [%default]')
parser.add_option('-m', '--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'orientation matrix label')
parser.add_option('-a',
dest = 'a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest = 'b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar = 'string',
help = 'crystal frame c vector label')
parser.add_option('-q', '--quaternion',
help = 'stress plane normal in lab frame, per default perpendicular to the force')
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
type = 'string', metavar = 'string',
help = 'quaternion label')
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.set_defaults(force = (0.0,0.0,1.0),
quaternion='orientation',
normal = None,
lattice = latticeChoices[0],
CoverA = math.sqrt(8./3.),
degrees = False,
)
(options, filenames) = parser.parse_args()
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
force = np.array(options.force)
force /= np.linalg.norm(force)
if options.normal:
if options.normal is not None:
normal = np.array(options.normal)
normal /= np.linalg.norm(normal)
if abs(np.dot(force,normal)) > 1e-3:
@ -174,22 +152,6 @@ if options.normal:
else:
normal = force
input = [options.eulers is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
]
if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested
slip_direction = np.zeros((len(slipSystems[options.lattice]),3),'f')
slip_normal = np.zeros_like(slip_direction)
@ -227,13 +189,12 @@ for name in filenames:
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
if not np.all(table.label_dimension(label) == dim):
damask.util.croak('input {} does not have dimension {}.'.format(label,dim))
if not table.label_dimension(options.quaternion) == 4:
damask.util.croak('input {} does not have dimension 4.'.format(options.quaternion))
table.close(dismiss = True) # close ASCIItable and remove empty file
continue
column = table.label_index(label)
column = table.label_index(options.quaternion)
# ------------------------------------------ assemble header ---------------------------------------
@ -251,17 +212,7 @@ for name in filenames:
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers':
o = damask.Orientation(Eulers = np.array(list(map(float,table.data[column:column+3])))*toRadians,)
elif inputtype == 'matrix':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column:column+9]))).reshape(3,3).transpose(),)
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3),)
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),)
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))))
table.data_append( np.abs( np.sum(slip_direction * (o.quaternion * force) ,axis=1) \
* np.sum(slip_normal * (o.quaternion * normal),axis=1)))

5
processing/post/binXY.py
View File

@ -1,4 +1,4 @@
#!/usr/bin/env python2.7
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys
@ -118,10 +118,9 @@ for name in filenames:
minmax[c] = np.log(minmax[c]) # change minmax to log, too
delta = minmax[:,1]-minmax[:,0]
(grid,xedges,yedges) = np.histogram2d(table.data[:,0],table.data[:,1],
bins=options.bins,
range=minmax,
range=minmax[:2],
weights=None if options.weight is None else table.data[:,2])
if options.normCol:

133
processing/post/postResults.py
View File

@ -121,12 +121,8 @@ class MPIEspectral_result: # mimic py_post result object
self._logscales = self._keyedPackedArray('logscales',count=self.N_loadcases,type='i')
self.size = self._keyedPackedArray('size:',count=3,type='d')
if self.size == [None,None,None]: # no 'size' found, try legacy alias 'dimension'
self.size = self._keyedPackedArray('dimension',count=3,type='d')
self.grid = self._keyedPackedArray('grid:',count=3,type='i')
if self.grid == [None,None,None]: # no 'grid' found, try legacy alias 'resolution'
self.grid = self._keyedPackedArray('resolution',count=3,type='i')
self.N_nodes = (self.grid[0]+1)*(self.grid[1]+1)*(self.grid[2]+1)
self.N_elements = self.grid[0] * self.grid[1] * self.grid[2]
@ -142,13 +138,8 @@ class MPIEspectral_result: # mimic py_post result object
# parameters for file handling depending on output format
if options.legacy:
self.tagLen=8
self.fourByteLimit = 2**31 -1 -8
else:
self.tagLen=0
self.tagLen=0
self.expectedFileSize = self.dataOffset+self.N_increments*(self.tagLen+self.N_elements*self.N_element_scalars*8)
if options.legacy: self.expectedFileSize+=self.expectedFileSize//self.fourByteLimit*8 # add extra 8 bytes for additional headers at 4 GB limits
if self.expectedFileSize != self.filesize:
print('\n**\n* Unexpected file size. Incomplete simulation or file corrupted!\n**')
@ -280,42 +271,16 @@ class MPIEspectral_result: # mimic py_post result object
return self.N_element_scalars
def element_scalar(self,e,idx):
if not options.legacy:
incStart = self.dataOffset \
+ self.position*8*self.N_elements*self.N_element_scalars
where = (e*self.N_element_scalars + idx)*8
try:
self.file.seek(incStart+where)
value = struct.unpack('d',self.file.read(8))[0]
except:
print('seeking {}'.format(incStart+where))
print('e {} idx {}'.format(e,idx))
sys.exit(1)
else:
self.fourByteLimit = 2**31 -1 -8
# header & footer + extra header and footer for 4 byte int range (Fortran)
# values
incStart = self.dataOffset \
+ self.position*8*( 1 + self.N_elements*self.N_element_scalars*8//self.fourByteLimit \
+ self.N_elements*self.N_element_scalars)
where = (e*self.N_element_scalars + idx)*8
try:
if where%self.fourByteLimit + 8 >= self.fourByteLimit: # danger of reading into fortran record footer at 4 byte limit
data=''
for i in range(8):
self.file.seek(incStart+where+(where//self.fourByteLimit)*8+4)
data += self.file.read(1)
where += 1
value = struct.unpack('d',data)[0]
else:
self.file.seek(incStart+where+(where//self.fourByteLimit)*8+4)
value = struct.unpack('d',self.file.read(8))[0]
except:
print('seeking {}'.format(incStart+where+(where//self.fourByteLimit)*8+4))
print('e {} idx {}'.format(e,idx))
sys.exit(1)
incStart = self.dataOffset \
+ self.position*8*self.N_elements*self.N_element_scalars
where = (e*self.N_element_scalars + idx)*8
try:
self.file.seek(incStart+where)
value = struct.unpack('d',self.file.read(8))[0]
except:
print('seeking {}'.format(incStart+where))
print('e {} idx {}'.format(e,idx))
sys.exit(1)
return [elemental_scalar(node,value) for node in self.element(e).items]
@ -645,8 +610,6 @@ of already processed data points for evaluation.
parser.add_option('-i','--info', action='store_true', dest='info',
help='list contents of resultfile')
parser.add_option('-l','--legacy', action='store_true', dest='legacy',
help='data format of spectral solver is in legacy format (no MPI out)')
parser.add_option('-n','--nodal', action='store_true', dest='nodal',
help='data is extrapolated to nodal value')
parser.add_option( '--prefix', dest='prefix',
@ -673,10 +636,7 @@ parser.add_option('-p','--type', dest='filetype',
help = 'type of result file [auto]')
parser.add_option('-q','--quiet', dest='verbose',
action = 'store_false',
help = 'suppress verbose output')
parser.add_option('--verbose', dest='verbose',
action = 'store_true',
help = 'enable verbose output')
help = 'hide status bar (useful when piping to file)')
group_material = OptionGroup(parser,'Material identifier')
@ -718,9 +678,8 @@ parser.add_option_group(group_general)
parser.add_option_group(group_special)
parser.set_defaults(info = False,
verbose = False,
legacy = False,
nodal = False,
verbose = True,
prefix = '',
suffix = '',
dir = 'postProc',
@ -747,6 +706,8 @@ if files == []:
parser.print_help()
parser.error('no file specified...')
damask.util.report(scriptName,files[0])
if not os.path.exists(files[0]):
parser.print_help()
parser.error('invalid file "%s" specified...'%files[0])
@ -803,12 +764,6 @@ if not options.constitutiveResult: options.constitutiveResult = []
options.sort.reverse()
options.sep.reverse()
# --- start background messaging
if options.verbose:
bg = damask.util.backgroundMessage()
bg.start()
# --- parse .output and .t16 files
if os.path.splitext(files[0])[1] == '':
@ -825,18 +780,13 @@ me = {
'Constitutive': options.phase,
}
if options.verbose: bg.set_message('parsing .output files...')
for what in me:
outputFormat[what] = ParseOutputFormat(filename, what, me[what])
if '_id' not in outputFormat[what]['specials']:
print("\nsection '{}' not found in <{}>".format(me[what], what))
print('\n'.join(map(lambda x:' [%s]'%x, outputFormat[what]['specials']['brothers'])))
if options.verbose: bg.set_message('opening result file...')
p = OpenPostfile(filename+extension,options.filetype,options.nodal)
if options.verbose: bg.set_message('parsing result file...')
stat = ParsePostfile(p, filename, outputFormat)
if options.filetype == 'marc':
stat['NumberOfIncrements'] -= 1 # t16 contains one "virtual" increment (at 0)
@ -879,8 +829,10 @@ if options.info:
# --- build connectivity maps
elementsOfNode = {}
for e in range(stat['NumberOfElements']):
if options.verbose and e%1000 == 0: bg.set_message('connect elem %i...'%e)
Nelems = stat['NumberOfElements']
for e in range(Nelems):
if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=e,total=Nelems,prefix='1/3: connecting elements')
for n in map(p.node_sequence,p.element(e).items):
if n not in elementsOfNode:
elementsOfNode[n] = [p.element_id(e)]
@ -899,10 +851,13 @@ index = {}
groups = []
groupCount = 0
memberCount = 0
damask.util.progressBar(iteration=1,total=1,prefix='1/3: connecting elements')
if options.nodalScalar:
for n in range(stat['NumberOfNodes']):
if options.verbose and n%1000 == 0: bg.set_message('scan node %i...'%n)
Npoints = stat['NumberOfNodes']
for n in range(Npoints):
if options.verbose and Npoints > 100 and e%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=n,total=Npoints,prefix='2/3: scanning nodes ')
myNodeID = p.node_id(n)
myNodeCoordinates = [p.node(n).x, p.node(n).y, p.node(n).z]
myElemID = 0
@ -911,32 +866,35 @@ if options.nodalScalar:
# generate an expression that is only true for the locations specified by options.filter
filter = substituteLocation(options.filter, [myElemID,myNodeID,myIpID,myGrainID], myNodeCoordinates)
if filter != '' and not eval(filter): # for all filter expressions that are not true:...
continue # ... ignore this data point and continue with next
if filter != '' and not eval(filter): # for all filter expressions that are not true:...
continue # ... ignore this data point and continue with next
# --- group data locations
# generate a unique key for a group of separated data based on the separation criterium for the location
grp = substituteLocation('#'.join(options.sep), [myElemID,myNodeID,myIpID,myGrainID], myNodeCoordinates)
if grp not in index: # create a new group if not yet present
if grp not in index: # create a new group if not yet present
index[grp] = groupCount
groups.append([[0,0,0,0,0.0,0.0,0.0]]) # initialize with avg location
groups.append([[0,0,0,0,0.0,0.0,0.0]]) # initialize with avg location
groupCount += 1
groups[index[grp]][0][:4] = mapIncremental('','unique',
len(groups[index[grp]])-1,
groups[index[grp]][0][:4],
[myElemID,myNodeID,myIpID,myGrainID]) # keep only if unique average location
[myElemID,myNodeID,myIpID,myGrainID]) # keep only if unique average location
groups[index[grp]][0][4:] = mapIncremental('','avg',
len(groups[index[grp]])-1,
groups[index[grp]][0][4:],
myNodeCoordinates) # incrementally update average location
groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member
myNodeCoordinates) # incrementally update average location
groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,0]) # append a new list defining each group member
memberCount += 1
damask.util.progressBar(iteration=1,total=1,prefix='2/3: scanning nodes ')
else:
for e in range(stat['NumberOfElements']):
if options.verbose and e%1000 == 0: bg.set_message('scan elem %i...'%e)
Nelems = stat['NumberOfElements']
for e in range(Nelems):
if options.verbose and Nelems > 100 and e%(Nelems//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=e,total=Nelems,prefix='2/3: scanning elements ')
myElemID = p.element_id(e)
myIpCoordinates = ipCoords(p.element(e).type, list(map(lambda node: [node.x, node.y, node.z],
list(map(p.node, map(p.node_sequence, p.element(e).items))))))
@ -976,6 +934,7 @@ else:
myIpCoordinates[n]) # incrementally update average location
groups[index[grp]].append([myElemID,myNodeID,myIpID,myGrainID,n]) # append a new list defining each group member
memberCount += 1
damask.util.progressBar(iteration=1,total=1,prefix='2/3: scanning elements ')
# --------------------------- sort groups --------------------------------
@ -1002,7 +961,6 @@ if 'none' not in map(str.lower, options.sort):
theKeys.append('x[0][%i]'%where[criterium])
sortKeys = eval('lambda x:(%s)'%(','.join(theKeys)))
if options.verbose: bg.set_message('sorting groups...')
groups.sort(key = sortKeys) # in-place sorting to save mem
@ -1021,8 +979,6 @@ standard = ['inc'] + \
# --------------------------- loop over positions --------------------------------
if options.verbose: bg.set_message('getting map between positions and increments...')
incAtPosition = {}
positionOfInc = {}
@ -1048,8 +1004,8 @@ increments = [incAtPosition[x] for x in locations] # build list of increments to
time_start = time.time()
Nincs = len([i for i in locations])
for incCount,position in enumerate(locations): # walk through locations
p.moveto(position+offset_pos) # wind to correct position
# --------------------------- file management --------------------------------
@ -1075,16 +1031,14 @@ for incCount,position in enumerate(locations): # walk through locations
# --------------------------- read and map data per group --------------------------------
member = 0
for group in groups:
Ngroups = len(groups)
for j,group in enumerate(groups):
f = incCount*Ngroups + j
if options.verbose and (Ngroups*Nincs) > 100 and f%((Ngroups*Nincs)//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=f,total=Ngroups*Nincs,prefix='3/3: processing points ')
N = 0 # group member counter
for (e,n,i,g,n_local) in group[1:]: # loop over group members
member += 1
if member%1000 == 0:
time_delta = ((len(locations)*memberCount)/float(member+incCount*memberCount)-1.0)*(time.time()-time_start)
if options.verbose: bg.set_message('(%02i:%02i:%02i) processing point %i of %i from increment %i (position %i)...'
%(time_delta//3600,time_delta%3600//60,time_delta%60,member,memberCount,increments[incCount],position))
newby = [] # current member's data
if options.nodalScalar:
@ -1172,6 +1126,7 @@ for incCount,position in enumerate(locations): # walk through locations
group[0] + \
mappedResult)
)) + '\n')
damask.util.progressBar(iteration=1,total=1,prefix='3/3: processing points ')
if fileOpen:
file.close()

65
processing/post/rotateData.py
View File

@ -18,19 +18,15 @@ Rotate vector and/or tensor column data by given angle around given axis.
""", version = scriptID)
parser.add_option('-v','--vector',
dest = 'vector',
parser.add_option('-d', '--data',
dest = 'data',
action = 'extend', metavar = '<string LIST>',
help = 'column heading of vector(s) to rotate')
parser.add_option('-t','--tensor',
dest = 'tensor',
action = 'extend', metavar = '<string LIST>',
help = 'column heading of tensor(s) to rotate')
help = 'vector/tensor value(s) label(s)')
parser.add_option('-r', '--rotation',
dest = 'rotation',
type = 'float', nargs = 4, metavar = ' '.join(['float']*4),
help = 'angle and axis to rotate data [%default]')
parser.add_option('-d', '--degrees',
parser.add_option('--degrees',
dest = 'degrees',
action = 'store_true',
help = 'angles are given in degrees [%default]')
@ -41,7 +37,7 @@ parser.set_defaults(rotation = (0.,1.,1.,1.),
(options,filenames) = parser.parse_args()
if options.vector is None and options.tensor is None:
if options.data is None:
parser.error('no data column specified.')
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians
@ -59,27 +55,24 @@ for name in filenames:
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header ------------------------------------------
# --- interpret header ----------------------------------------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks ----------------------------------------
items = {
'tensor': {'dim': 9, 'shape': [3,3], 'labels':options.tensor, 'active':[], 'column': []},
'vector': {'dim': 3, 'shape': [3], 'labels':options.vector, 'active':[], 'column': []},
}
errors = []
remarks = []
column = {}
for type, data in items.items():
for what in data['labels']:
dim = table.label_dimension(what)
if dim != data['dim']: remarks.append('column {} is not a {}.'.format(what,type))
else:
items[type]['active'].append(what)
items[type]['column'].append(table.label_index(what))
active = {'vector':[],'tensor':[]}
for i,dim in enumerate(table.label_dimension(options.data)):
label = options.data[i]
if dim == -1:
remarks.append('"{}" not found...'.format(label))
elif dim == 3:
remarks.append('adding vector "{}"...'.format(label))
active['vector'].append(label)
elif dim == 9:
remarks.append('adding tensor "{}"...'.format(label))
active['tensor'].append(label)
if remarks != []: damask.util.croak(remarks)
if errors != []:
@ -95,20 +88,14 @@ for name in filenames:
# ------------------------------------------ process data ------------------------------------------
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
datatype = 'vector'
for column in items[datatype]['column']: # loop over all requested labels
table.data[column:column+items[datatype]['dim']] = \
q * np.array(list(map(float,table.data[column:column+items[datatype]['dim']])))
datatype = 'tensor'
for column in items[datatype]['column']: # loop over all requested labels
table.data[column:column+items[datatype]['dim']] = \
np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+items[datatype]['dim']]))).\
reshape(items[datatype]['shape']),R.transpose())).reshape(items[datatype]['dim'])
for v in active['vector']:
column = table.label_index(v)
table.data[column:column+3] = q * np.array(list(map(float,table.data[column:column+3])))
for t in active['tensor']:
column = table.label_index(t)
table.data[column:column+9] = \
np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+9]))).reshape((3,3)),
R.transpose())).reshape((9))
outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization -----------------------------------

0
processing/post/showTable.py → processing/post/viewTable.py
View File

8
processing/pre/geom_addPrimitive.py
View File

@ -64,11 +64,11 @@ if options.dimension is None:
parser.error('no dimension specified.')
if options.angleaxis is not None:
options.angleaxis = list(map(float,options.angleaxis))
rotation = damask.Quaternion().fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0],
options.angleaxis[1:4])
rotation = damask.Quaternion.fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0],
options.angleaxis[1:4])
elif options.quaternion is not None:
options.quaternion = map(float,options.quaternion)
rotation = damask.Quaternion(options.quaternion)
options.quaternion = list(map(float,options.quaternion))
rotation = damask.Quaternion(quat=options.quaternion)
else:
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',
dest = 'degrees',
action = 'store_true',
help = 'angles are given in degrees [%default]')
help = 'all angles are in degrees')
parser.add_option('-m', '--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
@ -71,7 +71,7 @@ parser.add_option('--axes',
parser.add_option('-s', '--symmetry',
dest = 'symmetry',
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',
dest = 'homogenization',
type = 'int', metavar = 'int',
@ -234,7 +234,7 @@ for name in filenames:
o = damask.Orientation(Eulers = myData[colOri:colOri+3]*toRadians,
symmetry = mySym)
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)
elif inputtype == 'frame':
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],
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
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
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
closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself
match = 'local'
@ -269,7 +269,7 @@ for name in filenames:
if len(grains) > 0:
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
closest_grain = np.argmax(cos_disorientations) # grain among grains with closest orientation to myself
match = 'global'

5
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -323,12 +323,13 @@ for name in filenames:
]
if hasEulers:
config_header += ['<texture>']
theAxes = [] if options.axes is None else ['axes\t{} {} {}'.format(*options.axes)]
for ID in grainIDs:
eulerID = np.nonzero(grains == ID)[0][0] # find first occurrence of this grain id
config_header += ['[Grain{}]'.format(str(ID).zfill(formatwidth)),
'(gauss)\tphi1 {:g}\tPhi {:g}\tphi2 {:g}\tscatter 0.0\tfraction 1.0'.format(*eulers[eulerID])
]
if options.axes is not None: config_header.append('axes\t{} {} {}'.format(*options.axes))
] + theAxes
config_header += ['<!skip>']
table.labels_clear()
table.info_clear()

2
processing/pre/geom_rotate.py
View File

@ -60,8 +60,6 @@ eulers = np.array(damask.orientation.Orientation(
degrees = options.degrees,
).asEulers(degrees=True))
damask.util.croak('{} {} {}'.format(*eulers))
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]

2
processing/pre/hybridIA_linODFsampling.py
View File

@ -244,7 +244,7 @@ for name in filenames:
continue
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)
# ------------------------------------------ read header and data ---------------------------------

10
src/IO.f90
View File

@ -191,7 +191,9 @@ recursive function IO_recursiveRead(fileName,cnt) result(fileContent)
l,i, &
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
@ -684,7 +686,11 @@ function IO_stringValue(string,chunkPos,myChunk,silent)
logical :: warn
warn = merge(silent,.false.,present(silent))
if (present(silent)) then
warn = silent
else
warn = .false.
endif
IO_stringValue = ''
valuePresent: if (myChunk > chunkPos(1) .or. myChunk < 1_pInt) then

18
src/config.f90
View File

@ -142,7 +142,7 @@ subroutine config_init()
case (trim(material_partPhase))
call parseFile(phase_name,config_phase,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Phase parsed'; flush(6)
case (trim(material_partMicrostructure))
call parseFile(microstructure_name,config_microstructure,line,fileContent(i+1:))
@ -150,7 +150,7 @@ subroutine config_init()
case (trim(material_partCrystallite))
call parseFile(crystallite_name,config_crystallite,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Crystallite parsed'; flush(6)
case (trim(material_partHomogenization))
call parseFile(homogenization_name,config_homogenization,line,fileContent(i+1:))
@ -158,7 +158,7 @@ subroutine config_init()
case (trim(material_partTexture))
call parseFile(texture_name,config_texture,line,fileContent(i+1:))
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6)
if (iand(myDebug,debug_levelBasic) /= 0_pInt) write(6,'(a)') ' Texture parsed'; flush(6)
end select
@ -513,8 +513,12 @@ character(len=65536) function getString(this,key,defaultVal,raw)
type(tPartitionedStringList), pointer :: item
logical :: found, &
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)
if (found) then
getString = trim(defaultVal)
@ -661,7 +665,11 @@ function getStrings(this,key,defaultVal,requiredShape,raw)
cumulative
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.
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,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,2,phase)*(T - TRef)**2 / 2. &
+ 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
end subroutine math_check
!--------------------------------------------------------------------------------------------------
!> @brief Quicksort algorithm for two-dimensional integer arrays
@ -2625,12 +2625,9 @@ real(pReal) pure elemental function math_clip(a, left, right)
real(pReal), intent(in) :: a
real(pReal), intent(in), optional :: left, right
math_clip = min ( &
max (merge(left, -huge(a), present(left)), a), &
merge(right, huge(a), present(right)) &
)
math_clip = a
if (present(left)) math_clip = max(left,math_clip)
if (present(right)) math_clip = min(right,math_clip)
if (present(left) .and. present(right)) &
math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right)

120
src/plastic_phenopowerlaw.f90
View File

@ -24,12 +24,10 @@ module plastic_phenopowerlaw
accumulatedshear_slip_ID, &
shearrate_slip_ID, &
resolvedstress_slip_ID, &
totalshear_ID, &
resistance_twin_ID, &
accumulatedshear_twin_ID, &
shearrate_twin_ID, &
resolvedstress_twin_ID, &
totalvolfrac_twin_ID
resolvedstress_twin_ID
end enum
type, private :: tParameters
@ -55,7 +53,7 @@ module plastic_phenopowerlaw
xi_twin_0, & !< initial critical shear stress for twin
xi_slip_sat, & !< maximum critical shear stress for slip
nonSchmidCoeff, &
H_int, & !< per family hardening activity (optional) !ToDo: Better name!
H_int, & !< per family hardening activity (optional)
gamma_twin_char !< characteristic shear for twins
real(pReal), allocatable, dimension(:,:) :: &
interaction_SlipSlip, & !< slip resistance from slip activity
@ -80,9 +78,6 @@ module plastic_phenopowerlaw
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance)
type, private :: tPhenopowerlawState
real(pReal), pointer, dimension(:) :: &
sumGamma, & ! ToDo: why not make a dependent state?
sumF ! ToDo: why not make a dependent state?
real(pReal), pointer, dimension(:,:) :: &
xi_slip, &
xi_twin, &
@ -153,12 +148,6 @@ subroutine plastic_phenopowerlaw_init
real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::]
character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::]
type(tParameters) :: &
prm
type(tPhenopowerlawState) :: &
stt, &
dot
integer(kind(undefined_ID)) :: &
outputID !< ID of each post result output
@ -166,7 +155,7 @@ subroutine plastic_phenopowerlaw_init
structure = '',&
extmsg = ''
character(len=65536), dimension(:), allocatable :: &
outputs
outputs
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
@ -245,8 +234,8 @@ subroutine plastic_phenopowerlaw_init
! sanity checks
if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_slip '
if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//'a_slip ' ! ToDo: negative values ok?
if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//'n_slip ' ! ToDo: negative values ok?
if (dEq0(prm%a_slip)) extmsg = trim(extmsg)//'a_slip ' ! ToDo: negative values ok?
if (dEq0(prm%n_slip)) extmsg = trim(extmsg)//'n_slip ' ! ToDo: negative values ok?
if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//'xi_slip_0 '
if (any(prm%xi_slip_sat < prm%xi_slip_0)) extmsg = trim(extmsg)//'xi_slip_sat '
else slipActive
@ -279,10 +268,11 @@ subroutine plastic_phenopowerlaw_init
! sanity checks
if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_twin '
if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//'n_twin ' ! ToDo: negative values ok?
if (dEq0(prm%n_twin)) extmsg = trim(extmsg)//'n_twin ' ! ToDo: negative values ok?
else twinActive
allocate(prm%interaction_TwinTwin(0,0))
allocate(prm%xi_twin_0(0))
allocate(prm%gamma_twin_char(0))
endif twinActive
!--------------------------------------------------------------------------------------------------
@ -295,8 +285,8 @@ subroutine plastic_phenopowerlaw_init
config_phase(p)%getFloats('interaction_twinslip'), &
structure(1:3))
else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0
allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0
prm%h0_TwinSlip = 0.0_pReal
endif slipAndTwinActive
@ -338,12 +328,6 @@ subroutine plastic_phenopowerlaw_init
outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0_pInt)
outputSize = prm%totalNtwin
case ('totalshear')
outputID = merge(totalshear_ID,undefined_ID,prm%totalNslip>0_pInt)
outputSize = 1_pInt
case ('totalvolfrac_twin')
outputID = merge(totalvolfrac_twin_ID,undefined_ID,prm%totalNtwin>0_pInt)
outputSize = 1_pInt
end select
if (outputID /= undefined_ID) then
@ -358,8 +342,7 @@ subroutine plastic_phenopowerlaw_init
! allocate state arrays
NipcMyPhase = count(material_phase == p) ! number of IPCs containing my phase
sizeState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip &
+ size(['tau_twin ','gamma_twin']) * prm%TotalNtwin &
+ size(['sum(gamma)','sum(f) ']) ! ToDo: only needed if either twin or slip active!
+ size(['tau_twin ','gamma_twin']) * prm%TotalNtwin
sizeDotState = sizeState
call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, &
@ -383,18 +366,6 @@ subroutine plastic_phenopowerlaw_init
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance
startIndex = endIndex + 1_pInt
endIndex = endIndex + 1_pInt
stt%sumGamma => plasticState(p)%state (startIndex,:)
dot%sumGamma => plasticState(p)%dotState(startIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear
startIndex = endIndex + 1_pInt
endIndex = endIndex + 1_pInt
stt%sumF=>plasticState(p)%state (startIndex,:)
dot%sumF=>plasticState(p)%dotState(startIndex,:)
plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolTwinFrac
startIndex = endIndex + 1_pInt
endIndex = endIndex + prm%totalNslip
stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:)
@ -421,6 +392,8 @@ end subroutine plastic_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!> @details asumme that deformation by dislocation glide affects twinned and untwinned volume
! equally (Taylor assumption). Twinning happens only in untwinned volume (
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
@ -443,18 +416,15 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
gdot_slip_pos,gdot_slip_neg
real(pReal), dimension(param(instance)%totalNtwin) :: &
gdot_twin,dgdot_dtautwin
type(tParameters) :: prm
type(tPhenopowerlawState) :: stt
associate(prm => param(instance), stt => state(instance))
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
associate(prm => param(instance), stt => state(instance))
call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
slipSystems: do i = 1_pInt, prm%totalNslip
Lp = Lp + (1.0_pReal-stt%sumF(of))*(gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) &
@ -468,9 +438,9 @@ subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
+ dgdot_dtautwin(i)*prm%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i)
enddo twinSystems
end associate
end subroutine plastic_phenopowerlaw_LpAndItsTangent
@ -490,29 +460,28 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
i
real(pReal) :: &
c_SlipSlip,c_TwinSlip,c_TwinTwin, &
xi_slip_sat_offset
xi_slip_sat_offset,&
sumGamma,sumF
real(pReal), dimension(param(instance)%totalNslip) :: &
left_SlipSlip,right_SlipSlip, &
gdot_slip_pos,gdot_slip_neg
type(tParameters) :: prm
type(tPhenopowerlawState) :: dot,stt
associate(prm => param(instance), stt => state(instance), dot => dotState(instance))
dot%whole(:,of) = 0.0_pReal
sumGamma = sum(stt%gamma_slip(:,of))
sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)
!--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*stt%sumF(of)** prm%twinB)
c_TwinSlip = prm%h0_TwinSlip * stt%sumGamma(of)**prm%twinE
c_TwinTwin = prm%h0_TwinTwin * stt%sumF(of)**prm%twinD
c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*sumF** prm%twinB)
c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%twinE
c_TwinTwin = prm%h0_TwinTwin * sumF**prm%twinD
!--------------------------------------------------------------------------------------------------
! calculate left and right vectors
left_SlipSlip = 1.0_pReal + prm%H_int
xi_slip_sat_offset = prm%spr*sqrt(stt%sumF(of))
xi_slip_sat_offset = prm%spr*sqrt(sumF)
right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip &
* sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset))
@ -520,17 +489,13 @@ subroutine plastic_phenopowerlaw_dotState(Mp,instance,of)
! shear rates
call kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg)
dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg)
dot%sumGamma(of) = sum(dot%gamma_slip(:,of))
call kinetics_twin(prm,stt,of,Mp,dot%gamma_twin(:,of))
if (prm%totalNtwin > 0_pInt) dot%sumF(of) = merge(sum(dot%gamma_twin(:,of)/prm%gamma_twin_char), &
0.0_pReal, &
stt%sumF(of) < 0.98_pReal)
!--------------------------------------------------------------------------------------------------
! hardening
hardeningSlip: do i = 1_pInt, prm%totalNslip
dot%xi_slip(i,of) = dot_product(prm%interaction_SlipSlip(i,:),right_SlipSlip*dot%gamma_slip(:,of)) &
* c_SlipSlip * left_SlipSlip(i) &
* c_SlipSlip * left_SlipSlip(i) &
+ dot_product(prm%interaction_SlipTwin(i,:),dot%gamma_twin(:,of))
enddo hardeningSlip
@ -546,8 +511,9 @@ end subroutine plastic_phenopowerlaw_dotState
!--------------------------------------------------------------------------------------------------
!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress
!> @details: Shear rates are calculated only optionally. NOTE: Against the common convention, the
!> result (i.e. intent(out)) variables are the last to have the optional arguments at the end
!> @details Shear rates are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_slip(prm,stt,of,Mp,gdot_slip_pos,gdot_slip_neg, &
dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
@ -619,9 +585,11 @@ end subroutine kinetics_slip
!--------------------------------------------------------------------------------------------------
!> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress
!> @details: Shear rates are calculated only optionally. NOTE: Against the common convention, the
!> result (i.e. intent(out)) variables are the last to have the optional arguments at the end
!> @brief calculates shear rates on twin systems and derivatives with respect to resolved stress.
! twinning is assumed to take place only in untwinned volume.
!> @details Derivates are calculated only optionally.
! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
! have the optional arguments at the end.
!--------------------------------------------------------------------------------------------------
pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin)
use prec, only: &
@ -652,7 +620,8 @@ pure subroutine kinetics_twin(prm,stt,of,Mp,gdot_twin,dgdot_dtau_twin)
enddo
where(tau_twin > 0.0_pReal)
gdot_twin = (1.0_pReal-stt%sumF(of))*prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction
* prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin
else where
gdot_twin = 0.0_pReal
end where
@ -690,13 +659,10 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
real(pReal), dimension(param(instance)%totalNslip) :: &
gdot_slip_pos,gdot_slip_neg
type(tParameters) :: prm
type(tPhenopowerlawState) :: stt
associate( prm => param(instance), stt => state(instance))
postResults = 0.0_pReal
c = 0_pInt
associate( prm => param(instance), stt => state(instance))
outputsLoop: do o = 1_pInt,size(prm%outputID)
select case(prm%outputID(o))
@ -732,15 +698,9 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults)
enddo
c = c + prm%totalNtwin
case (totalshear_ID)
postResults(c+1_pInt) = stt%sumGamma(of)
c = c + 1_pInt
case (totalvolfrac_twin_ID)
postResults(c+1_pInt) = stt%sumF(of)
c = c + 1_pInt
end select
enddo outputsLoop
end associate
end function plastic_phenopowerlaw_postResults

6
src/spectral_damage.f90
View File

@ -120,8 +120,8 @@ subroutine spectral_damage_init()
trim(snes_type) == 'vinewtonssls') then
call DMGetGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMGetGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
call VecSet(lBound,0.0,ierr); CHKERRQ(ierr)
call VecSet(uBound,1.0,ierr); CHKERRQ(ierr)
call VecSet(lBound,0.0_pReal,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 DMRestoreGlobalVector(damage_grid,lBound,ierr); CHKERRQ(ierr)
call DMRestoreGlobalVector(damage_grid,uBound,ierr); CHKERRQ(ierr)
@ -134,7 +134,7 @@ subroutine spectral_damage_init()
xend = xstart + xend - 1
yend = ystart + yend - 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_lastInc(grid(1),grid(2),grid3), source=1.0_pReal)
allocate(damage_stagInc(grid(1),grid(2),grid3), source=1.0_pReal)