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 NewStyleKinematicHardening-2

This commit is contained in:
Martin Diehl 2018-12-14 06:00:28 +01:00
parent 20671b8ed3 c038dc9f93
commit 85f1368480
28 changed files with 1751 additions and 5777 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

52
.gitlab-ci.yml
View File

@ -151,13 +151,20 @@ Spectral_geometryPacking:
- release
###################################################################################################
Post_General:
Post_AverageDown:
stage: postprocessing
script: PostProcessing/test.py
script: averageDown/test.py
except:
- master
- release
#Post_General:
# stage: postprocessing
# script: PostProcessing/test.py
# except:
# - master
# - release
Post_GeometryReconstruction:
stage: postprocessing
script: Spectral_geometryReconstruction/test.py
@ -193,6 +200,13 @@ Post_ParaviewRelated:
- master
- release
Post_OrientationConversion:
stage: postprocessing
script: OrientationConversion/test.py
except:
- master
- release
###################################################################################################
Compile_Spectral_Intel:
stage: compilePETScIntel
@ -243,6 +257,13 @@ Compile_Intel_Prepare:
- release
###################################################################################################
Thermal:
stage: spectral
script: Thermal/test.py
except:
- master
- release
Spectral_PackedGeometry:
stage: spectral
script: Spectral_PackedGeometry/test.py
@ -343,19 +364,12 @@ Phenopowerlaw_singleSlip:
- master
- release
HybridIA:
stage: spectral
script: HybridIA/test.py
except:
- master
- release
TextureComponents:
stage: spectral
script: TextureComponents/test.py
except:
- master
- release
#TextureComponents:
# stage: spectral
# script: TextureComponents/test.py
# except:
# - master
# - release
###################################################################################################
@ -423,14 +437,6 @@ SpectralExample:
only:
- development
AbaqusExample:
stage: example
script:
- module load $IntelAbaqus $Abaqus
- Abaqus_example/test.py
only:
- development
###################################################################################################
SpectralRuntime:
stage: performance

2
CMakeLists.txt
View File

@ -47,7 +47,7 @@ linker:
# CMake will execute each target in the ${petsc_config_makefile}
# to acquire corresponding PETSc Variables.
find_program (MAKE_EXECUTABLE NAMES make gmake)
find_program (MAKE_EXECUTABLE NAMES gmake make)
# Find the PETSc includes directory settings
execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "includes"
RESULT_VARIABLE PETSC_INCLUDES_RETURN

3
DAMASK_prerequisites.sh
View File

@ -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 dd70ef078e0b4422b4bf81030e896a7bd098e841
Subproject commit e9f93abaecafbfbf11072ae70bca213a7201ed38

2
VERSION
View File

@ -1 +1 @@
v2.0.2-955-g7663d20b
v2.0.2-1124-g0ff9e9fa

3
env/DAMASK.csh vendored
View File

@ -44,7 +44,6 @@ if ( $?prompt ) then
echo "DAMASK $DAMASK_ROOT"
echo "Spectral Solver $SOLVER"
echo "Post Processing $PROCESSING"
echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
if ( $?PETSC_DIR) then
echo "PETSc location $PETSC_DIR"
endif
@ -52,8 +51,10 @@ if ( $?prompt ) then
echo "MSC.Marc/Mentat $MSC_ROOT"
endif
echo
echo "Multithreading DAMASK_NUM_THREADS=$DAMASK_NUM_THREADS"
echo `limit datasize`
echo `limit stacksize`
echo
endif
setenv DAMASK_NUM_THREADS $DAMASK_NUM_THREADS

1
env/DAMASK.sh vendored
View File

@ -88,6 +88,7 @@ if [ ! -z "$PS1" ]; then
size=$(( 1024*$(ulimit -s) )); \
print('{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \
['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0]))")
echo
fi
export DAMASK_NUM_THREADS

1
env/DAMASK.zsh vendored
View File

@ -81,6 +81,7 @@ if [ ! -z "$PS1" ]; then
size=$(( 1024*$(ulimit -s) )); \
print('{:.4g} {}'.format(size / (1 << ((int(math.log(size,2) / 10) if size else 0) * 10)), \
['bytes','KiB','MiB','GiB','TiB','EiB','ZiB'][int(math.log(size,2) / 10) if size else 0]))")
echo
fi
export DAMASK_NUM_THREADS

2832
examples/AbaqusStandard/Job_sx-px.inp
View File

File diff suppressed because it is too large Load Diff

BIN
examples/AbaqusStandard/SX_PX_compression.cae
View File

Binary file not shown.

23
examples/AbaqusStandard/SX_PX_compression.jnl
View File

@ -1,23 +0,0 @@
# Save by m.diehl on 2017_12_06-18.38.26; build 2017 2016_09_27-23.54.59 126836
from abaqus import *
upgradeMdb(
'/nethome/storage/raid4/m.diehl/DAMASK/examples/AbaqusStandard/SX_PX_compression-6.9-1.cae'
,
'/nethome/storage/raid4/m.diehl/DAMASK/examples/AbaqusStandard/SX_PX_compression.cae')
# Save by m.diehl on 2017_12_06-18.38.26; build 2017 2016_09_27-23.54.59 126836
from part import *
from material import *
from section import *
from assembly import *
from step import *
from interaction import *
from load import *
from mesh import *
from optimization import *
from job import *
from sketch import *
from visualization import *
from connectorBehavior import *
mdb.jobs['Job_sx-px'].setValues(description='compression', userSubroutine=
'$HOME/DAMASK/src/DAMASK_abaqus_std.f')
# Save by m.diehl on 2017_12_06-18.39.44; build 2017 2016_09_27-23.54.59 126836

93
examples/AbaqusStandard/material.config
View File

@ -1,93 +0,0 @@
#-------------------#
<homogenization>
#-------------------#
[dummy]
mech none
[poly]
mech isostrain
Nconstituents 10
#-------------------#
<microstructure>
#-------------------#
[Aluminum_001]
crystallite 1
(constituent) phase 1 texture 1 fraction 1.0
[Aluminum_10]
crystallite 1
(constituent) phase 1 texture 2 fraction 1.0
#-------------------#
<crystallite>
#-------------------#
[orientation]
(output) eulerangles
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4)
#-------------------#
<phase>
#-------------------#
[Aluminum_phenopowerlaw]
# slip only
elasticity hooke
plasticity phenopowerlaw
(output) resistance_slip
(output) shearrate_slip
(output) resolvedstress_slip
(output) totalshear
(output) resistance_twin
(output) shearrate_twin
(output) resolvedstress_twin
(output) totalvolfrac
lattice_structure fcc
Nslip 12
Ntwin 0
c11 106.75e9
c12 60.41e9
c44 28.34e9
gdot0_slip 0.001
n_slip 20
tau0_slip 31e6 # per family
tausat_slip 63e6 # per family
a_slip 2.25
gdot0_twin 0.001
n_twin 20
tau0_twin 31e6 # per family
s_pr 0 # push-up factor for slip saturation due to twinning
twin_b 0
twin_c 0
twin_d 0
twin_e 0
h0_slipslip 75e6
h0_sliptwin 0
h0_twinslip 0
h0_twintwin 0
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
atol_resistance 1
#-------------------#
<texture>
#-------------------#
[001]
(gauss) phi1 0.000 Phi 45.000 phi2 0.000 scatter 0.000 fraction 1.000
[random]

1
examples/AbaqusStandard/numerics.config
View File

@ -1 +0,0 @@
fixed_seed 1697667030

358
lib/damask/orientation.py
View File

@ -33,7 +33,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,206 +48,138 @@ 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"""
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
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
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
return self.__class__(q=self.q*other,
p=self.p*other)
except:
return self.copy()
def __imul__(self, other):
"""In-place multiplication"""
# Rowenhorst_etal2015 MSMSE: value of P is selected as -1
P = -1.0
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
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)
return (self-other).magnitude() < 1e-8 or (-self-other).magnitude() < 1e-8
def __ne__(self,other):
"""Not equal at e-8 precision"""
@ -255,31 +187,26 @@ class Quaternion:
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
return self.q ** 2 + np.dot(self.p,self.p)
def identity(self):
self.w = 1.
self.x = 0.
self.y = 0.
self.z = 0.
self.q = 1.
self.p = np.zeros(3,dtype=float)
return self
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
self.p = -self.p
return self
def inverse(self):
@ -290,11 +217,9 @@ class Quaternion:
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):
@ -310,27 +235,35 @@ class Quaternion:
return self.copy().homomorph()
def asList(self):
return [i for i in self]
return [self.q]+list(self.p)
def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.)
return np.outer([i for i in self],[i for i in self])
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:
if self.q > 1.:
self.normalize()
s = math.sqrt(1. - self.w**2)
x = 2*self.w**2 - 1.
y = 2*self.w * s
s = math.sqrt(1. - self.q**2)
x = 2*self.q**2 - 1.
y = 2*self.q * s
angle = math.atan2(y,x)
if angle < 0.0:
@ -338,26 +271,28 @@ class Quaternion:
s *= -1.
return (np.degrees(angle) if degrees else angle,
np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else [self.x / s, self.y / s, self.z / s]))
np.array([1.0, 0.0, 0.0] if np.abs(angle) < 1e-6 else self.p / s))
def asRodrigues(self):
return np.inf*np.ones(3) if self.w == 0.0 else np.array([self.x, self.y, self.z])/self.w
return np.inf*np.ones(3) if self.q == 0.0 else self.p/self.q
def asEulers(self,
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])
eulers = np.array([math.atan2(-2*P*self.q*self.p[2],self.q**2-self.p[2]**2),0,0])
elif abs(chi) < 1e-10 and abs(q03) < 1e-10:
eulers = np.array([math.atan2( 2*self.x*self.y,self.x**2-self.y**2),np.pi,0])
eulers = np.array([math.atan2( 2 *self.p[0]*self.p[1],self.p[0]**2-self.p[1]**2),np.pi,0])
else:
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),
])
return np.degrees(eulers) if degrees else eulers
@ -371,25 +306,26 @@ 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])
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 +333,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 +353,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 +370,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(1.+m[0,0]+m[1,1]+m[2,2])
x = P*0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
y = P*0.5*math.sqrt(1.-m[0,0]+m[1,1]-m[2,2])
z = P*0.5*math.sqrt(1.-m[0,0]-m[1,1]+m[2,2])
x *= -1 if m[2,1] < m[1,2] else 1
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 +395,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
@ -513,7 +444,7 @@ class Symmetry:
def __repr__(self):
"""Readbable string"""
return '%s' % (self.lattice)
return '{}'.format(self.lattice)
def __eq__(self, other):
@ -526,7 +457,7 @@ class Symmetry:
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)
@ -722,7 +653,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 +668,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 +711,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 +726,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()
@ -927,7 +861,7 @@ class Orientation:
M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa
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)

4
processing/post/addGrainID.py
View File

@ -200,9 +200,9 @@ for name in filenames:
if gID != -1 and gID not in alreadyChecked: # indexed point belonging to a grain not yet tested?
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

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'

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 ---------------------------------

2
processing/pre/seeds_fromRandom.py
View File

@ -127,7 +127,7 @@ options.fraction = np.array(options.fraction)
options.grid = np.array(options.grid)
gridSize = options.grid.prod()
if options.randomSeed is None: options.randomSeed = int(os.urandom(4).encode('hex'), 16)
if options.randomSeed is None: options.randomSeed = int(os.urandom(4).hex(), 16)
np.random.seed(options.randomSeed) # init random generators
random.seed(options.randomSeed)

230
src/IO.f90
View File

@ -33,7 +33,6 @@ module IO
IO_write_jobIntFile, &
IO_read_realFile, &
IO_read_intFile, &
IO_hybridIA, &
IO_isBlank, &
IO_getTag, &
IO_stringPos, &
@ -192,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
@ -583,223 +584,6 @@ logical function IO_abaqus_hasNoPart(fileUnit)
620 end function IO_abaqus_hasNoPart
#endif
!--------------------------------------------------------------------------------------------------
!> @brief hybrid IA sampling of ODFfile
!--------------------------------------------------------------------------------------------------
function IO_hybridIA(Nast,ODFfileName)
use prec, only: &
tol_math_check
implicit none
integer(pInt), intent(in) :: Nast !< number of samples?
real(pReal), dimension(3,Nast) :: IO_hybridIA
character(len=*), intent(in) :: ODFfileName !< name of ODF file including total path
!--------------------------------------------------------------------------------------------------
! math module is not available
real(pReal), parameter :: PI = 3.141592653589793_pReal
real(pReal), parameter :: INRAD = PI/180.0_pReal
integer(pInt) :: i,j,bin,NnonZero,Nset,Nreps,reps,phi1,Phi,phi2
integer(pInt), allocatable, dimension(:) :: chunkPos
integer(pInt), dimension(3) :: steps !< number of steps in phi1, Phi, and phi2 direction
integer(pInt), dimension(4) :: columns !< columns in linearODF file where eulerangles and density are located
integer(pInt), dimension(:), allocatable :: binSet
real(pReal) :: center,sum_dV_V,prob,dg_0,C,lowerC,upperC,rnd
real(pReal), dimension(2,3) :: limits !< starting and end values for eulerangles
real(pReal), dimension(3) :: deltas, & !< angular step size in phi1, Phi, and phi2 direction
eulers !< euler angles when reading from file
real(pReal), dimension(:,:,:), allocatable :: dV_V
character(len=65536) :: line, keyword
integer(pInt) :: headerLength
integer(pInt), parameter :: FILEUNIT = 999_pInt
IO_hybridIA = 0.0_pReal ! initialize return value for case of error
write(6,'(/,a,/)',advance='no') ' Using linear ODF file: '//trim(ODFfileName)
write(6,'(/,a)') ' Eisenlohr et al., Computational Materials Science, 42(4):670–678, 2008'
write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2007.09.015'
!--------------------------------------------------------------------------------------------------
! parse header of ODF file
call IO_open_file(FILEUNIT,ODFfileName)
headerLength = 0_pInt
line=IO_read(FILEUNIT)
chunkPos = IO_stringPos(line)
keyword = IO_lc(IO_StringValue(line,chunkPos,2_pInt,.true.))
if (keyword(1:4) == 'head') then
headerLength = IO_intValue(line,chunkPos,1_pInt) + 1_pInt
else
call IO_error(error_ID=156_pInt, ext_msg='no header found')
endif
!--------------------------------------------------------------------------------------------------
! figure out columns containing data
do i = 1_pInt, headerLength-1_pInt
line=IO_read(FILEUNIT)
enddo
columns = 0_pInt
chunkPos = IO_stringPos(line)
do i = 1_pInt, chunkPos(1)
select case ( IO_lc(IO_StringValue(line,chunkPos,i,.true.)) )
case ('phi1')
columns(1) = i
case ('phi')
columns(2) = i
case ('phi2')
columns(3) = i
case ('intensity')
columns(4) = i
end select
enddo
if (any(columns<1)) call IO_error(error_ID = 156_pInt, ext_msg='could not find expected header')
!--------------------------------------------------------------------------------------------------
! determine limits, number of steps and step size
limits(1,1:3) = 721.0_pReal
limits(2,1:3) = -1.0_pReal
steps = 0_pInt
line=IO_read(FILEUNIT)
do while (trim(line) /= IO_EOF)
chunkPos = IO_stringPos(line)
eulers=[IO_floatValue(line,chunkPos,columns(1)),&
IO_floatValue(line,chunkPos,columns(2)),&
IO_floatValue(line,chunkPos,columns(3))]
steps = steps + merge(1,0,eulers>limits(2,1:3))
limits(1,1:3) = min(limits(1,1:3),eulers)
limits(2,1:3) = max(limits(2,1:3),eulers)
line=IO_read(FILEUNIT)
enddo
deltas = (limits(2,1:3)-limits(1,1:3))/real(steps-1_pInt,pReal)
write(6,'(/,a,/,3(2x,f12.4,1x))',advance='no') ' Starting angles / ° = ',limits(1,1:3)
write(6,'(/,a,/,3(2x,f12.4,1x))',advance='no') ' Ending angles / ° = ',limits(2,1:3)
write(6,'(/,a,/,3(2x,f12.4,1x))',advance='no') ' Angular steps / ° = ',deltas
if (all(abs(limits(1,1:3)) < tol_math_check)) then
write(6,'(/,a,/)',advance='no') ' assuming vertex centered data'
center = 0.0_pReal ! no need to shift
if (any(mod(int(limits(2,1:3),pInt),90)==0)) &
call IO_error(error_ID = 156_pInt, ext_msg='linear ODF data repeated at right boundary')
else
write(6,'(/,a,/)',advance='no') ' assuming cell centered data'
center = 0.5_pReal ! shift data by half of a bin
endif
limits = limits*INRAD
deltas = deltas*INRAD
!--------------------------------------------------------------------------------------------------
! read in data
allocate(dV_V(steps(3),steps(2),steps(1)),source=0.0_pReal)
sum_dV_V = 0.0_pReal
dg_0 = deltas(1)*deltas(3)*2.0_pReal*sin(deltas(2)/2.0_pReal)
NnonZero = 0_pInt
call IO_checkAndRewind(FILEUNIT) ! forward
do i = 1_pInt, headerLength
line=IO_read(FILEUNIT)
enddo
do phi1=1_pInt,steps(1); do Phi=1_pInt,steps(2); do phi2=1_pInt,steps(3)
line=IO_read(FILEUNIT)
chunkPos = IO_stringPos(line)
eulers=[IO_floatValue(line,chunkPos,columns(1)),& ! read in again for consistency check only
IO_floatValue(line,chunkPos,columns(2)),&
IO_floatValue(line,chunkPos,columns(3))]*INRAD
if (any(abs((real([phi1,phi,phi2],pReal) -1.0_pReal + center)*deltas-eulers)>tol_math_check)) & ! check if data is in expected order (phi2 fast) and correct for Fortran starting at 1
call IO_error(error_ID = 156_pInt, ext_msg='linear ODF data not in expected order')
prob = IO_floatValue(line,chunkPos,columns(4))
if (prob > 0.0_pReal) then
NnonZero = NnonZero+1_pInt
sum_dV_V = sum_dV_V+prob
else
prob = 0.0_pReal
endif
dV_V(phi2,Phi,phi1) = prob*dg_0*sin((real(Phi-1_pInt,pReal)+center)*deltas(2))
enddo; enddo; enddo
close(FILEUNIT)
dV_V = dV_V/sum_dV_V ! normalize to 1
!--------------------------------------------------------------------------------------------------
! now fix bounds
Nset = max(Nast,NnonZero) ! if less than non-zero voxel count requested, sample at least that much
lowerC = 0.0_pReal
upperC = real(Nset, pReal)
do while (hybridIA_reps(dV_V,steps,upperC) < Nset)
lowerC = upperC
upperC = upperC*2.0_pReal
enddo
!--------------------------------------------------------------------------------------------------
! binary search for best C
do
C = (upperC+lowerC)/2.0_pReal
Nreps = hybridIA_reps(dV_V,steps,C)
if (abs(upperC-lowerC) < upperC*1.0e-14_pReal) then
C = upperC
Nreps = hybridIA_reps(dV_V,steps,C)
exit
elseif (Nreps < Nset) then
lowerC = C
elseif (Nreps > Nset) then
upperC = C
else
exit
endif
enddo
allocate(binSet(Nreps))
bin = 0_pInt ! bin counter
i = 1_pInt ! set counter
do phi1=1_pInt,steps(1); do Phi=1_pInt,steps(2) ;do phi2=1_pInt,steps(3)
reps = nint(C*dV_V(phi2,Phi,phi1), pInt)
binSet(i:i+reps-1) = bin
bin = bin+1_pInt ! advance bin
i = i+reps ! advance set
enddo; enddo; enddo
do i=1_pInt,Nast
if (i < Nast) then
call random_number(rnd)
j = nint(rnd*real(Nreps-i,pReal)+real(i,pReal)+0.5_pReal,pInt)
else
j = i
endif
bin = binSet(j)
IO_hybridIA(1,i) = deltas(1)*(real(mod(bin/(steps(3)*steps(2)),steps(1)),pReal)+center) ! phi1
IO_hybridIA(2,i) = deltas(2)*(real(mod(bin/ steps(3) ,steps(2)),pReal)+center) ! Phi
IO_hybridIA(3,i) = deltas(3)*(real(mod(bin ,steps(3)),pReal)+center) ! phi2
binSet(j) = binSet(i)
enddo
contains
!--------------------------------------------------------------------------------------------------
!> @brief counts hybrid IA repetitions
!--------------------------------------------------------------------------------------------------
integer(pInt) pure function hybridIA_reps(dV_V,steps,C)
implicit none
integer(pInt), intent(in), dimension(3) :: steps !< number of bins in Euler space
real(pReal), intent(in), dimension(steps(3),steps(2),steps(1)) :: dV_V !< needs description
real(pReal), intent(in) :: C !< needs description
integer(pInt) :: phi1,Phi,phi2
hybridIA_reps = 0_pInt
do phi1=1_pInt,steps(1); do Phi =1_pInt,steps(2); do phi2=1_pInt,steps(3)
hybridIA_reps = hybridIA_reps+nint(C*dV_V(phi2,Phi,phi1), pInt)
enddo; enddo; enddo
end function hybridIA_reps
end function IO_hybridIA
!--------------------------------------------------------------------------------------------------
!> @brief identifies strings without content
@ -902,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
@ -1758,7 +1546,6 @@ integer(pInt) function IO_verifyIntValue (string,validChars,myName)
validChars, & !< valid characters in string
myName !< name of caller function (for debugging)
integer(pInt) :: readStatus, invalidWhere
!character(len=len(trim(string))) :: trimmed does not work with ifort 14.0.1
IO_verifyIntValue = 0_pInt
@ -1788,7 +1575,6 @@ real(pReal) function IO_verifyFloatValue (string,validChars,myName)
myName !< name of caller function (for debugging)
integer(pInt) :: readStatus, invalidWhere
!character(len=len(trim(string))) :: trimmed does not work with ifort 14.0.1
IO_verifyFloatValue = 0.0_pReal

12
src/config.f90
View File

@ -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

22
src/constitutive.f90
View File

@ -513,7 +513,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
of = phasememberAt(ipc,ip,el)
@ -526,9 +526,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
temperature(ho)%p(tme),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
@ -925,8 +925,9 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
call plastic_kinehardening_dotState(Mp,instance,of)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), &
ipc,ip,el)
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_dislotwin_dotState (Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_dotState (math_Mandel33to6(Mp),temperature(ho)%p(tme), &
@ -1143,22 +1144,29 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
case (PLASTICITY_ISOTROPIC_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_isotropic_postResults(S6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_postResults(startPos:endPos) = &
plastic_kinehardening_postResults(Mp,instance,of)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_postResults(startPos:endPos) = &
plastic_dislotwin_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
plastic_dislotwin_postResults(Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_disloucla_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_nonlocal_postResults (S6,FeArray,ip,el)

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. &

73
src/lattice.f90
View File

@ -72,7 +72,7 @@ module lattice
!--------------------------------------------------------------------------------------------------
! face centered cubic
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &
LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for fcc
LATTICE_fcc_NslipSystem = int([12, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0],pInt) !< # of slip systems per family for fcc
integer(pInt), dimension(LATTICE_maxNtwinFamily), parameter, public :: &
LATTICE_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< # of twin systems per family for fcc
@ -104,11 +104,19 @@ module lattice
1, 1, 0, 1,-1,-1, & ! A6
0, 1, 1, -1, 1,-1, & ! D1
1, 0,-1, -1, 1,-1, & ! D4
-1,-1, 0, -1, 1,-1 & ! D6
-1,-1, 0, -1, 1,-1, & ! D6
! Slip system <110>{110}
1, 1, 0, 1,-1, 0, &
1,-1, 0, 1, 1, 0, &
1, 0, 1, 1, 0,-1, &
1, 0,-1, 1, 0, 1, &
0, 1, 1, 0, 1,-1, &
0, 1,-1, 0, 1, 1 &
],pReal),shape(LATTICE_FCC_SYSTEMSLIP)) !< Slip system <110>{111} directions. Sorted according to Eisenlohr & Hantcherli
character(len=*), dimension(1), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = &
['<0 1 -1>{1 1 1}']
character(len=*), dimension(2), parameter, public :: LATTICE_FCC_SLIPFAMILY_NAME = &
['<0 1 -1>{1 1 1}', &
'<0 1 -1>{0 1 1}']
real(pReal), dimension(3+3,LATTICE_fcc_Ntwin), parameter, private :: &
LATTICE_fcc_systemTwin = reshape(real( [&
@ -166,25 +174,38 @@ module lattice
integer(pInt), dimension(LATTICE_fcc_Nslip,lattice_fcc_Nslip), parameter, public :: &
LATTICE_fcc_interactionSlipSlip = reshape(int( [&
1,2,2,4,6,5,3,5,5,4,5,6, & ! ---> slip
2,1,2,6,4,5,5,4,6,5,3,5, & ! |
2,2,1,5,5,3,5,6,4,6,5,4, & ! |
4,6,5,1,2,2,4,5,6,3,5,5, & ! v slip
6,4,5,2,1,2,5,3,5,5,4,6, &
5,5,3,2,2,1,6,5,4,5,6,4, &
3,5,5,4,5,6,1,2,2,4,6,5, &
5,4,6,5,3,5,2,1,2,6,4,5, &
5,6,4,6,5,4,2,2,1,5,5,3, &
4,5,6,3,5,5,4,6,5,1,2,2, &
5,3,5,5,4,6,6,4,5,2,1,2, &
6,5,4,5,6,4,5,5,3,2,2,1 &
],pInt),shape(LATTICE_FCC_INTERACTIONSLIPSLIP),order=[2,1]) !< Slip--slip interaction types for fcc
1, 2, 2, 4, 6, 5, 3, 5, 5, 4, 5, 6, 9,10, 9,10,11,12, & ! ---> slip
2, 1, 2, 6, 4, 5, 5, 4, 6, 5, 3, 5, 9,10,11,12, 9,10, & ! |
2, 2, 1, 5, 5, 3, 5, 6, 4, 6, 5, 4, 11,12, 9,10, 9,10, & ! |
4, 6, 5, 1, 2, 2, 4, 5, 6, 3, 5, 5, 9,10,10, 9,12,11, & ! v slip
6, 4, 5, 2, 1, 2, 5, 3, 5, 5, 4, 6, 9,10,12,11,10, 9, &
5, 5, 3, 2, 2, 1, 6, 5, 4, 5, 6, 4, 11,12,10, 9,10, 9, &
3, 5, 5, 4, 5, 6, 1, 2, 2, 4, 6, 5, 10, 9,10, 9,11,12, &
5, 4, 6, 5, 3, 5, 2, 1, 2, 6, 4, 5, 10, 9,12,11, 9,10, &
5, 6, 4, 6, 5, 4, 2, 2, 1, 5, 5, 3, 12,11,10, 9, 9,10, &
4, 5, 6, 3, 5, 5, 4, 6, 5, 1, 2, 2, 10, 9, 9,10,12,11, &
5, 3, 5, 5, 4, 6, 6, 4, 5, 2, 1, 2, 10, 9,11,12,10, 9, &
6, 5, 4, 5, 6, 4, 5, 5, 3, 2, 2, 1, 12,11, 9,10,10, 9, &
9, 9,11, 9, 9,11,10,10,12,10,10,12, 1, 7, 8, 8, 8, 8, &
10,10,12,10,10,12, 9, 9,11, 9, 9,11, 7, 1, 8, 8, 8, 8, &
9,11, 9,10,12,10,10,12,10, 9,11, 9, 8, 8, 1, 7, 8, 8, &
10,12,10, 9,11, 9, 9,11, 9,10,12,10, 8, 8, 7, 1, 8, 8, &
11, 9, 9,12,10,10,11, 9, 9,12,10,10, 8, 8, 8, 8, 1, 7, &
12,10,10,11, 9, 9,12,10,10,11, 9, 9, 8, 8, 8, 8, 7, 1 &
],pInt),[LATTICE_fcc_Nslip,LATTICE_fcc_Nslip],order=[2,1]) !< Slip--slip interaction types for fcc
!< 1: self interaction
!< 2: coplanar interaction
!< 3: collinear interaction
!< 4: Hirth locks
!< 5: glissile junctions
!< 6: Lomer locks
!< 7: crossing (similar to Hirth locks in <110>{111} for two {110} planes)
!< 8: similar to Lomer locks in <110>{111} for two {110} planes
!< 9: similar to Lomer locks in <110>{111} btw one {110} and one {111} plane
!<10: similar to glissile junctions in <110>{111} btw one {110} and one {111} plane
!<11: crossing btw one {110} and one {111} plane
!<12: collinear btw one {110} and one {111} plane
integer(pInt), dimension(LATTICE_fcc_Nslip,LATTICE_fcc_Ntwin), parameter, public :: &
LATTICE_fcc_interactionSlipTwin = reshape(int( [&
1,1,1,3,3,3,2,2,2,3,3,3, & ! ---> twin
@ -198,7 +219,14 @@ module lattice
3,3,3,3,3,3,1,1,1,2,2,2, &
3,3,3,2,2,2,3,3,3,1,1,1, &
2,2,2,3,3,3,3,3,3,1,1,1, &
3,3,3,3,3,3,2,2,2,1,1,1 &
3,3,3,3,3,3,2,2,2,1,1,1, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4 &
],pInt),shape(LATTICE_FCC_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip--twin interaction types for fcc
!< 1: coplanar interaction
!< 2: screw trace between slip system and twin habit plane (easy cross slip)
@ -235,7 +263,14 @@ module lattice
3,3,3,3,3,3,1,1,1,2,2,2, &
3,3,3,2,2,2,3,3,3,1,1,1, &
2,2,2,3,3,3,3,3,3,1,1,1, &
3,3,3,3,3,3,2,2,2,1,1,1 &
3,3,3,3,3,3,2,2,2,1,1,1, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4 &
],pInt),shape(LATTICE_FCCTOHEX_INTERACTIONSLIPTRANS),order=[2,1]) !< Slip--trans interaction types for fcc
integer(pInt), dimension(LATTICE_fcc_Ntrans,LATTICE_fcc_Nslip), parameter, public :: &

173
src/material.f90
View File

@ -169,18 +169,13 @@ module material
homogenization_maxNgrains !< max number of grains in any USED homogenization
integer(pInt), dimension(:), allocatable, public, protected :: &
material_homogenizationAt, & !< homogenization ID of each element (copy of mesh_homogenizationAt)
phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
phase_Noutput, & !< number of '(output)' items per phase
phase_elasticityInstance, & !< instance of particular elasticity of each phase
phase_plasticityInstance !< instance of particular plasticity of each phase
integer(pInt), dimension(:), allocatable, public, protected :: &
crystallite_Noutput !< number of '(output)' items per crystallite setting
integer(pInt), dimension(:), allocatable, public, protected :: &
phase_plasticityInstance, & !< instance of particular plasticity of each phase
crystallite_Noutput, & !< number of '(output)' items per crystallite setting
homogenization_Ngrains, & !< number of grains in each homogenization
homogenization_Noutput, & !< number of '(output)' items per homogenization
homogenization_typeInstance, & !< instance of particular type of each homogenization
@ -189,7 +184,7 @@ module material
vacancyflux_typeInstance, & !< instance of particular type of each vacancy flux
porosity_typeInstance, & !< instance of particular type of each porosity model
hydrogenflux_typeInstance, & !< instance of particular type of each hydrogen flux
microstructure_crystallite !< crystallite setting ID of each microstructure
microstructure_crystallite !< crystallite setting ID of each microstructure ! DEPRECATED !!!!
real(pReal), dimension(:), allocatable, public, protected :: &
thermal_initialT, & !< initial temperature per each homogenization
@ -198,12 +193,27 @@ module material
porosity_initialPhi, & !< initial posority per each homogenization
hydrogenflux_initialCh !< initial hydrogen concentration per each homogenization
! NEW MAPPINGS
integer(pInt), dimension(:), allocatable, public, protected :: &
material_homogenizationAt, & !< homogenization ID of each element (copy of mesh_homogenizationAt)
material_homogenizationMemberAt, & !< position of the element within its homogenization instance
material_aggregateAt, & !< aggregate ID of each element FUTURE USE FOR OUTPUT
material_aggregatMemberAt !< position of the element within its aggregate instance FUTURE USE FOR OUTPUT
integer(pInt), dimension(:,:), allocatable, public, protected :: &
material_phaseAt, & !< phase ID of each element
material_phaseMemberAt, & !< position of the element within its phase instance
material_crystalliteAt, & !< crystallite ID of each element CURRENTLY NOT PER CONSTITUTENT
material_crystalliteMemberAt !< position of the element within its crystallite instance CURRENTLY NOT PER CONSTITUTENT
! END NEW MAPPINGS
! DEPRECATED: use material_phaseAt
integer(pInt), dimension(:,:,:), allocatable, public :: &
material_phase !< phase (index) of each grain,IP,element
! BEGIN DEPRECATED: use material_homogenizationAt
! DEPRECATED: use material_homogenizationAt
integer(pInt), dimension(:,:), allocatable, public :: &
material_homog !< homogenization (index) of each IP,element
! END DEPRECATED
type(tPlasticState), allocatable, dimension(:), public :: &
plasticState
type(tSourceState), allocatable, dimension(:), public :: &
@ -227,10 +237,6 @@ module material
microstructure_elemhomo, & !< flag to indicate homogeneous microstructure distribution over element's IPs
phase_localPlasticity !< flags phases with local constitutive law
character(len=65536), dimension(:), allocatable, private :: &
texture_ODFfile !< name of each ODF file
integer(pInt), private :: &
microstructure_maxNconstituents, & !< max number of constituents in any phase
texture_maxNgauss, & !< max number of Gauss components in any texture
@ -258,11 +264,13 @@ module material
logical, dimension(:), allocatable, private :: &
homogenization_active
! BEGIN DEPRECATED
integer(pInt), dimension(:,:,:), allocatable, public :: phaseAt !< phase ID of every material point (ipc,ip,el)
integer(pInt), dimension(:,:,:), allocatable, public :: phasememberAt !< memberID of given phase at every material point (ipc,ip,el)
integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingCrystallite
integer(pInt), dimension(:,:,:), allocatable, public, target :: mappingHomogenization !< mapping from material points to offset in heterogenous state/field
integer(pInt), dimension(:,:), allocatable, public, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
integer(pInt), dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
! END DEPRECATED
type(tHomogMapping), allocatable, dimension(:), public :: &
thermalMapping, & !< mapping for thermal state/fields
@ -368,7 +376,6 @@ subroutine material_init()
use mesh, only: &
mesh_homogenizationAt, &
mesh_NipsPerElem, &
mesh_maxNips, &
mesh_NcpElems, &
FE_geomtype
@ -378,11 +385,10 @@ subroutine material_init()
integer(pInt) :: &
g, & !< grain number
i, & !< integration point number
e, & !< element number
phase
integer(pInt), dimension(:), allocatable :: ConstitutivePosition
integer(pInt), dimension(:), allocatable :: CrystallitePosition
integer(pInt), dimension(:), allocatable :: HomogenizationPosition
e !< element number
integer(pInt), dimension(:), allocatable :: &
CounterPhase, &
CounterHomogenization
myDebug = debug_level(debug_material)
@ -473,30 +479,34 @@ subroutine material_init()
call material_populateGrains
allocate(phaseAt ( homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
allocate(phasememberAt ( homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0_pInt)
allocate(mappingHomogenization (2, mesh_maxNips,mesh_NcpElems),source=0_pInt)
allocate(mappingCrystallite (2,homogenization_maxNgrains, mesh_NcpElems),source=0_pInt)
allocate(mappingHomogenizationConst( mesh_maxNips,mesh_NcpElems),source=1_pInt)
! BEGIN DEPRECATED
allocate(phaseAt ( homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
allocate(phasememberAt ( homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
allocate(mappingHomogenization (2, mesh_nIPsPerElem,mesh_NcpElems),source=0_pInt)
allocate(mappingHomogenizationConst( mesh_nIPsPerElem,mesh_NcpElems),source=1_pInt)
! END DEPRECATED
allocate(ConstitutivePosition (size(config_phase)), source=0_pInt)
allocate(HomogenizationPosition(size(config_homogenization)),source=0_pInt)
allocate(CrystallitePosition (size(config_phase)), source=0_pInt)
allocate(material_homogenizationAt,source=mesh_homogenizationAt)
allocate(CounterPhase (size(config_phase)), source=0_pInt)
allocate(CounterHomogenization(size(config_homogenization)),source=0_pInt)
ElemLoop:do e = 1_pInt,mesh_NcpElems
! BEGIN DEPRECATED
do e = 1_pInt,mesh_NcpElems
myHomog = mesh_homogenizationAt(e)
IPloop:do i = 1_pInt, mesh_NipsPerElem
HomogenizationPosition(myHomog) = HomogenizationPosition(myHomog) + 1_pInt
mappingHomogenization(1:2,i,e) = [HomogenizationPosition(myHomog),myHomog]
GrainLoop:do g = 1_pInt,homogenization_Ngrains(myHomog)
phase = material_phase(g,i,e)
ConstitutivePosition(phase) = ConstitutivePosition(phase)+1_pInt ! not distinguishing between instances of same phase
phaseAt(g,i,e) = phase
phasememberAt(g,i,e) = ConstitutivePosition(phase)
enddo GrainLoop
enddo IPloop
enddo ElemLoop
do i = 1_pInt, mesh_NipsPerElem
CounterHomogenization(myHomog) = CounterHomogenization(myHomog) + 1_pInt
mappingHomogenization(1:2,i,e) = [CounterHomogenization(myHomog),myHomog]
do g = 1_pInt,homogenization_Ngrains(myHomog)
myPhase = material_phase(g,i,e)
CounterPhase(myPhase) = CounterPhase(myPhase)+1_pInt ! not distinguishing between instances of same phase
phaseAt(g,i,e) = myPhase
phasememberAt(g,i,e) = CounterPhase(myPhase)
enddo
enddo
enddo
! END DEPRECATED
! REMOVE !!!!!
! hack needed to initialize field values used during constitutive and crystallite initializations
do myHomog = 1,size(config_homogenization)
thermalMapping (myHomog)%p => mappingHomogenizationConst
@ -937,9 +947,7 @@ subroutine material_parseTexture
integer(pInt) :: section, gauss, fiber, j, t, i
character(len=65536), dimension(:), allocatable :: strings ! Values for given key in material config
integer(pInt), dimension(:), allocatable :: chunkPos
character(len=65536) :: tag
allocate(texture_ODFfile(size(config_texture))); texture_ODFfile=''
allocate(texture_symmetry(size(config_texture)), source=1_pInt)
allocate(texture_Ngauss(size(config_texture)), source=0_pInt)
allocate(texture_Nfiber(size(config_texture)), source=0_pInt)
@ -985,9 +993,6 @@ subroutine material_parseTexture
if(dNeq(math_det33(texture_transformation(1:3,1:3,t)),1.0_pReal)) call IO_error(157_pInt,t)
endif
tag=''
texture_ODFfile(t) = config_texture(t)%getString('hybridia',defaultVal=tag)
if (config_texture(t)%keyExists('symmetry')) then
select case (config_texture(t)%getString('symmetry'))
case('orthotropic')
@ -1122,7 +1127,7 @@ end subroutine material_allocatePlasticState
!--------------------------------------------------------------------------------------------------
!> @brief populates the grains
!> @details populates the grains by identifying active microstructure/homogenization pairs,
!! calculates the volume of the grains and deals with texture components and hybridIA
!! calculates the volume of the grains and deals with texture components
!--------------------------------------------------------------------------------------------------
subroutine material_populateGrains
use prec, only: &
@ -1141,7 +1146,6 @@ subroutine material_populateGrains
mesh_elemType, &
mesh_homogenizationAt, &
mesh_microstructureAt, &
mesh_maxNips, &
mesh_NcpElems, &
mesh_ipVolume, &
FE_geomtype
@ -1152,8 +1156,7 @@ subroutine material_populateGrains
homogenization_name, &
microstructure_name
use IO, only: &
IO_error, &
IO_hybridIA
IO_error
use debug, only: &
debug_level, &
debug_material, &
@ -1181,12 +1184,11 @@ subroutine material_populateGrains
myDebug = debug_level(debug_material)
allocate(material_volume(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0.0_pReal)
allocate(material_phase(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt)
allocate(material_homog(mesh_maxNips,mesh_NcpElems), source=0_pInt)
allocate(material_homogenizationAt,source=mesh_homogenizationAt)
allocate(material_texture(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), source=0_pInt)
allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
allocate(material_volume(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems), source=0.0_pReal)
allocate(material_phase(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems), source=0_pInt)
allocate(material_homog(mesh_nIPsPerElem,mesh_NcpElems), source=0_pInt)
allocate(material_texture(homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems), source=0_pInt)
allocate(material_EulerAngles(3,homogenization_maxNgrains,mesh_nIPsPerElem,mesh_NcpElems),source=0.0_pReal)
allocate(Ngrains(size(config_homogenization),size(config_microstructure)), source=0_pInt)
allocate(Nelems (size(config_homogenization),size(config_microstructure)), source=0_pInt)
@ -1330,39 +1332,31 @@ subroutine material_populateGrains
real(texture_symmetry(textureID),pReal),pInt) ! max number of unique orientations (excl. symmetry)
!--------------------------------------------------------------------------------------------------
! ...has texture components
if (texture_ODFfile(textureID) == '') then
gauss: do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components
do g = 1_pInt,int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
texture_Gauss( 4,t,textureID))
enddo
constituentGrain = &
constituentGrain + int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID)) ! advance counter for grains of current constituent
enddo gauss
! has texture components
gauss: do t = 1_pInt,texture_Ngauss(textureID) ! loop over Gauss components
do g = 1_pInt,int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID),pInt) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleGaussOri(texture_Gauss(1:3,t,textureID),&
texture_Gauss( 4,t,textureID))
enddo
constituentGrain = &
constituentGrain + int(real(myNorientations,pReal)*texture_Gauss(5,t,textureID)) ! advance counter for grains of current constituent
enddo gauss
fiber: do t = 1_pInt,texture_Nfiber(textureID) ! loop over fiber components
do g = 1_pInt,int(real(myNorientations,pReal)*texture_Fiber(6,t,textureID),pInt) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
texture_Fiber(3:4,t,textureID),&
texture_Fiber( 5,t,textureID))
enddo
constituentGrain = &
constituentGrain + int(real(myNorientations,pReal)*texture_fiber(6,t,textureID),pInt) ! advance counter for grains of current constituent
enddo fiber
fiber: do t = 1_pInt,texture_Nfiber(textureID) ! loop over fiber components
do g = 1_pInt,int(real(myNorientations,pReal)*texture_Fiber(6,t,textureID),pInt) ! loop over required grain count
orientationOfGrain(:,grain+constituentGrain+g) = &
math_sampleFiberOri(texture_Fiber(1:2,t,textureID),&
texture_Fiber(3:4,t,textureID),&
texture_Fiber( 5,t,textureID))
enddo
constituentGrain = &
constituentGrain + int(real(myNorientations,pReal)*texture_fiber(6,t,textureID),pInt) ! advance counter for grains of current constituent
enddo fiber
random: do constituentGrain = constituentGrain+1_pInt,myNorientations ! fill remainder with random
orientationOfGrain(:,grain+constituentGrain) = math_sampleRandomOri()
enddo random
!--------------------------------------------------------------------------------------------------
! ...has hybrid IA
else
orientationOfGrain(1:3,grain+1_pInt:grain+myNorientations) = &
IO_hybridIA(myNorientations,texture_ODFfile(textureID))
if (all(dEq(orientationOfGrain(1:3,grain+1_pInt),-1.0_pReal))) call IO_error(156_pInt)
endif
random: do constituentGrain = constituentGrain+1_pInt,myNorientations ! fill remainder with random
orientationOfGrain(:,grain+constituentGrain) = math_sampleRandomOri()
enddo random
!--------------------------------------------------------------------------------------------------
! ...texture transformation
@ -1479,12 +1473,7 @@ subroutine material_populateGrains
enddo microstructureLoop
enddo homogenizationLoop
deallocate(volumeOfGrain)
deallocate(phaseOfGrain)
deallocate(textureOfGrain)
deallocate(orientationOfGrain)
deallocate(texture_transformation)
deallocate(Nelems)
deallocate(elemsOfHomogMicro)
call config_deallocate('material.config/microstructure')

9
src/math.f90
View File

@ -2625,12 +2625,9 @@ real(pReal) pure 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)

3522
src/plastic_dislotwin.f90
View File

File diff suppressed because it is too large Load Diff

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)