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

Merge remote-tracking branch 'origin/development' into 19-NewStylePhenopowerlaw

This commit is contained in:
Martin Diehl 2018-11-25 13:39:34 +01:00
parent 1808b37357 5290a9ac83
commit 8ccec40e35
16 changed files with 184 additions and 3417 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

36
.gitlab-ci.yml
View File

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

2
CMakeLists.txt
View File

@ -47,7 +47,7 @@ linker:
# CMake will execute each target in the ${petsc_config_makefile} # CMake will execute each target in the ${petsc_config_makefile}
# to acquire corresponding PETSc Variables. # 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 # Find the PETSc includes directory settings
execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "includes" execute_process (COMMAND ${MAKE_EXECUTABLE} --no-print-directory -f ${petsc_config_makefile} "includes"
RESULT_VARIABLE PETSC_INCLUDES_RETURN RESULT_VARIABLE PETSC_INCLUDES_RETURN

2
PRIVATE

@ -1 +1 @@
Subproject commit ee5a63d34abbac295207354fddf30e6d7cc258cd Subproject commit 7f22579bf986795e307c471834eeff0806d0ac84

2
VERSION
View File

@ -1 +1 @@
v2.0.2-933-ge9bf0824 v2.0.2-976-gde4873e9

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

177
lib/damask/orientation.py
View File

@ -27,15 +27,22 @@ class Rodrigues:
# ****************************************************************************************** # ******************************************************************************************
class Quaternion: class Quaternion:
""" u"""
Orientation represented as unit quaternion. Orientation represented as unit quaternion.
All methods and naming conventions based on http://www.euclideanspace.com/maths/algebra/realNormedAlgebra/quaternions. 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
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π]
Convention 5: the rotation angle ω is limited to the interval [0, π]
w is the real part, (x, y, z) are the imaginary parts. w is the real part, (x, y, z) are the imaginary parts.
Representation of rotation is in ACTIVE form!
(Derived directly or through angleAxis, Euler angles, or active matrix) Vector "a" (defined in coordinate system "A") is passively rotated
Vector "a" (defined in coordinate system "A") is actively rotated to new coordinates "b". resulting in new coordinates "b" when expressed in system "B".
b = Q * a b = Q * a
b = np.dot(Q.asMatrix(),a) b = np.dot(Q.asMatrix(),a)
""" """
@ -309,10 +316,12 @@ class Quaternion:
return np.outer([i for i in self],[i for i in self]) return np.outer([i for i in self],[i for i in self])
def asMatrix(self): def asMatrix(self):
return np.array( qbarhalf = 0.5*(self.w**2 - self.x**2 - self.y**2 - self.z**2)
[[1.0-2.0*(self.y*self.y+self.z*self.z), 2.0*(self.x*self.y-self.z*self.w), 2.0*(self.x*self.z+self.y*self.w)], return 2.0*np.array(
[ 2.0*(self.x*self.y+self.z*self.w), 1.0-2.0*(self.x*self.x+self.z*self.z), 2.0*(self.y*self.z-self.x*self.w)], [[ qbarhalf + self.x**2 , self.x*self.y - self.w*self.z, self.x*self.z + self.w*self.y],
[ 2.0*(self.x*self.z-self.y*self.w), 2.0*(self.x*self.w+self.y*self.z), 1.0-2.0*(self.x*self.x+self.y*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 ],
])
def asAngleAxis(self, def asAngleAxis(self,
degrees = False): degrees = False):
@ -335,52 +344,23 @@ class Quaternion:
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.w == 0.0 else np.array([self.x, self.y, self.z])/self.w
def asEulers(self, def asEulers(self,
type = "bunge", degrees = False):
degrees = False, """Orientation as Bunge-Euler angles."""
standardRange = False): q03 = self.w**2+self.z**2
""" q12 = self.x**2+self.y**2
Orientation as Bunge-Euler angles. chi = np.sqrt(q03*q12)
Conversion of ACTIVE rotation to Euler angles taken from: if abs(chi) < 1e-10 and abs(q12) < 1e-10:
Melcher, A.; Unser, A.; Reichhardt, M.; Nestler, B.; Poetschke, M.; Selzer, M. eulers = np.array([math.atan2(-2*self.w*self.z,self.w**2-self.z**2),0,0])
Conversion of EBSD data by a quaternion based algorithm to be used for grain structure simulations elif abs(chi) < 1e-10 and abs(q03) < 1e-10:
Technische Mechanik 30 (2010) pp 401--413. eulers = np.array([math.atan2( 2*self.x*self.y,self.x**2-self.y**2),np.pi,0])
"""
angles = [0.0,0.0,0.0]
if type.lower() == 'bunge' or type.lower() == 'zxz':
if abs(self.x) < 1e-4 and abs(self.y) < 1e-4:
x = self.w**2 - self.z**2
y = 2.*self.w*self.z
angles[0] = math.atan2(y,x)
elif abs(self.w) < 1e-4 and abs(self.z) < 1e-4:
x = self.x**2 - self.y**2
y = 2.*self.x*self.y
angles[0] = math.atan2(y,x)
angles[1] = math.pi
else: else:
chi = math.sqrt((self.w**2 + self.z**2)*(self.x**2 + self.y**2)) eulers = np.array([math.atan2((self.x*self.z-self.w*self.y)/chi,(-self.w*self.x-self.y*self.z)/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),
])
x = (self.w * self.x - self.y * self.z)/2./chi return np.degrees(eulers) if degrees else eulers
y = (self.w * self.y + self.x * self.z)/2./chi
angles[0] = math.atan2(y,x)
x = self.w**2 + self.z**2 - (self.x**2 + self.y**2)
y = 2.*chi
angles[1] = math.atan2(y,x)
x = (self.w * self.x + self.y * self.z)/2./chi
y = (self.z * self.x - self.y * self.w)/2./chi
angles[2] = math.atan2(y,x)
if standardRange:
angles[0] %= 2*math.pi
if angles[1] < 0.0:
angles[1] += math.pi
angles[2] *= -1.0
angles[2] %= 2*math.pi
return np.degrees(angles) if degrees else angles
# # Static constructors # # Static constructors
@ -408,7 +388,7 @@ class Quaternion:
halfangle = math.atan(np.linalg.norm(rodrigues)) halfangle = math.atan(np.linalg.norm(rodrigues))
c = math.cos(halfangle) c = math.cos(halfangle)
w = c w = c
x,y,z = c*rodrigues x,y,z = rodrigues/c
return cls([w,x,y,z]) return cls([w,x,y,z])
@ -431,24 +411,19 @@ class Quaternion:
@classmethod @classmethod
def fromEulers(cls, def fromEulers(cls,
eulers, eulers,
type = 'Bunge',
degrees = False): degrees = False):
if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d') if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype='d')
eulers = np.radians(eulers) if degrees else eulers eulers = np.radians(eulers) if degrees else eulers
c = np.cos(0.5 * eulers) sigma = 0.5*(eulers[0]+eulers[2])
s = np.sin(0.5 * eulers) delta = 0.5*(eulers[0]-eulers[2])
c = np.cos(0.5*eulers[1])
s = np.sin(0.5*eulers[1])
if type.lower() == 'bunge' or type.lower() == 'zxz': w = c * np.cos(sigma)
w = c[0] * c[1] * c[2] - s[0] * c[1] * s[2] x = -s * np.cos(delta)
x = c[0] * s[1] * c[2] + s[0] * s[1] * s[2] y = -s * np.sin(delta)
y = - c[0] * s[1] * s[2] + s[0] * s[1] * c[2] z = -c * np.sin(sigma)
z = c[0] * c[1] * s[2] + s[0] * c[1] * c[2]
else:
w = c[0] * c[1] * c[2] - s[0] * s[1] * s[2]
x = s[0] * s[1] * c[2] + c[0] * c[1] * s[2]
y = s[0] * c[1] * c[2] + c[0] * s[1] * s[2]
z = c[0] * s[1] * c[2] - s[0] * c[1] * s[2]
return cls([w,x,y,z]) return cls([w,x,y,z])
@ -460,49 +435,16 @@ class Quaternion:
if m.shape != (3,3) and np.prod(m.shape) == 9: if m.shape != (3,3) and np.prod(m.shape) == 9:
m = m.reshape(3,3) m = m.reshape(3,3)
tr = np.trace(m) w = 0.5*math.sqrt(1.+m[0,0]+m[1,1]+m[2,2])
if tr > 1e-8: x = 0.5*math.sqrt(1.+m[0,0]-m[1,1]-m[2,2])
s = math.sqrt(tr + 1.0)*2.0 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])
return cls( x *= -1 if m[2,1] < m[1,2] else 1
[ s*0.25, y *= -1 if m[0,2] < m[2,0] else 1
(m[2,1] - m[1,2])/s, z *= -1 if m[1,0] < m[0,1] else 1
(m[0,2] - m[2,0])/s,
(m[1,0] - m[0,1])/s,
])
elif m[0,0] > m[1,1] and m[0,0] > m[2,2]: return cls( np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2))
t = m[0,0] - m[1,1] - m[2,2] + 1.0
s = 2.0*math.sqrt(t)
return cls(
[ (m[2,1] - m[1,2])/s,
s*0.25,
(m[0,1] + m[1,0])/s,
(m[2,0] + m[0,2])/s,
])
elif m[1,1] > m[2,2]:
t = -m[0,0] + m[1,1] - m[2,2] + 1.0
s = 2.0*math.sqrt(t)
return cls(
[ (m[0,2] - m[2,0])/s,
(m[0,1] + m[1,0])/s,
s*0.25,
(m[1,2] + m[2,1])/s,
])
else:
t = -m[0,0] - m[1,1] + m[2,2] + 1.0
s = 2.0*math.sqrt(t)
return cls(
[ (m[1,0] - m[0,1])/s,
(m[2,0] + m[0,2])/s,
(m[1,2] + m[2,1])/s,
s*0.25,
])
@classmethod @classmethod
@ -663,7 +605,7 @@ class Symmetry:
quaternion, quaternion,
who = []): who = []):
"""List of symmetrically equivalent quaternions based on own symmetry.""" """List of symmetrically equivalent quaternions based on own symmetry."""
return [quaternion*q for q in self.symmetryQuats(who)] return [q*quaternion for q in self.symmetryQuats(who)]
def inFZ(self,R): def inFZ(self,R):
@ -829,7 +771,7 @@ class Orientation:
else: else:
self.quaternion = Quaternion.fromRandom(randomSeed=random) self.quaternion = Quaternion.fromRandom(randomSeed=random)
elif isinstance(Eulers, np.ndarray) and Eulers.shape == (3,): # based on given Euler angles elif isinstance(Eulers, np.ndarray) and Eulers.shape == (3,): # based on given Euler angles
self.quaternion = Quaternion.fromEulers(Eulers,type='bunge',degrees=degrees) self.quaternion = Quaternion.fromEulers(Eulers,degrees=degrees)
elif isinstance(matrix, np.ndarray) : # based on given rotation matrix elif isinstance(matrix, np.ndarray) : # based on given rotation matrix
self.quaternion = Quaternion.fromMatrix(matrix) self.quaternion = Quaternion.fromMatrix(matrix)
elif isinstance(angleAxis, np.ndarray) and angleAxis.shape == (4,): # based on given angle and rotation axis elif isinstance(angleAxis, np.ndarray) and angleAxis.shape == (4,): # based on given angle and rotation axis
@ -855,16 +797,15 @@ class Orientation:
return 'Symmetry: %s\n' % (self.symmetry) + \ return 'Symmetry: %s\n' % (self.symmetry) + \
'Quaternion: %s\n' % (self.quaternion) + \ 'Quaternion: %s\n' % (self.quaternion) + \
'Matrix:\n%s\n' % ( '\n'.join(['\t'.join(map(str,self.asMatrix()[i,:])) for i in range(3)]) ) + \ '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('bunge',degrees=True))) ) 'Bunge Eulers / deg: %s' % ('\t'.join(map(str,self.asEulers(degrees=True))) )
def asQuaternion(self): def asQuaternion(self):
return self.quaternion.asList() return self.quaternion.asList()
def asEulers(self, def asEulers(self,
type = 'bunge',
degrees = False, degrees = False,
standardRange = False): ):
return self.quaternion.asEulers(type, degrees, standardRange) return self.quaternion.asEulers(degrees)
eulers = property(asEulers) eulers = property(asEulers)
def asRodrigues(self): def asRodrigues(self):
@ -912,13 +853,13 @@ class Orientation:
""" """
if self.symmetry != other.symmetry: raise TypeError('disorientation between different symmetry classes not supported yet.') if self.symmetry != other.symmetry: raise TypeError('disorientation between different symmetry classes not supported yet.')
misQ = self.quaternion.conjugated()*other.quaternion misQ = other.quaternion*self.quaternion.conjugated()
mySymQs = self.symmetry.symmetryQuats() if SST else self.symmetry.symmetryQuats()[:1] # take all or only first sym operation mySymQs = self.symmetry.symmetryQuats() if SST else self.symmetry.symmetryQuats()[:1] # take all or only first sym operation
otherSymQs = other.symmetry.symmetryQuats() otherSymQs = other.symmetry.symmetryQuats()
for i,sA in enumerate(mySymQs): for i,sA in enumerate(mySymQs):
for j,sB in enumerate(otherSymQs): for j,sB in enumerate(otherSymQs):
theQ = sA.conjugated()*misQ*sB theQ = sB*misQ*sA.conjugated()
for k in range(2): for k in range(2):
theQ.conjugate() theQ.conjugate()
breaker = self.symmetry.inFZ(theQ) \ breaker = self.symmetry.inFZ(theQ) \
@ -939,10 +880,10 @@ class Orientation:
"""Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)""" """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)"""
if SST: # pole requested to be within SST if SST: # pole requested to be within SST
for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)): # test all symmetric equivalent quaternions for i,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)): # test all symmetric equivalent quaternions
pole = q.conjugated()*axis # align crystal direction to axis pole = q*axis # align crystal direction to axis
if self.symmetry.inSST(pole,proper): break # found SST version if self.symmetry.inSST(pole,proper): break # found SST version
else: else:
pole = self.quaternion.conjugated()*axis # align crystal direction to axis pole = self.quaternion*axis # align crystal direction to axis
return (pole,i if SST else 0) return (pole,i if SST else 0)
@ -951,7 +892,7 @@ class Orientation:
color = np.zeros(3,'d') color = np.zeros(3,'d')
for q in self.symmetry.equivalentQuaternions(self.quaternion): for q in self.symmetry.equivalentQuaternions(self.quaternion):
pole = q.conjugated()*axis # align crystal direction to axis pole = q*axis # align crystal direction to axis
inSST,color = self.symmetry.inSST(pole,color=True) inSST,color = self.symmetry.inSST(pole,color=True)
if inSST: break if inSST: break

2
processing/post/addMises.py
View File

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

2
processing/post/addOrientations.py
View File

@ -169,7 +169,7 @@ for name in filenames:
for output in options.output: for output in options.output:
if output == 'quaternion': table.data_append(o.asQuaternion()) if output == 'quaternion': table.data_append(o.asQuaternion())
elif output == 'rodrigues': table.data_append(o.asRodrigues()) elif output == 'rodrigues': table.data_append(o.asRodrigues())
elif output == 'eulers': table.data_append(o.asEulers('Bunge', degrees=options.degrees)) elif output == 'eulers': table.data_append(o.asEulers(degrees=options.degrees))
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------

26
processing/post/addSchmidfactors.py
View File

@ -190,27 +190,27 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.')
(options.quaternion,4,'quaternion'), (options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested ][np.where(input)[0][0]] # select input label that was requested
c_direction = np.zeros((len(slipSystems[options.lattice]),3),'f') slip_direction = np.zeros((len(slipSystems[options.lattice]),3),'f')
c_normal = np.zeros_like(c_direction) slip_normal = np.zeros_like(slip_direction)
if options.lattice in latticeChoices[:2]: if options.lattice in latticeChoices[:2]:
c_direction = slipSystems[options.lattice][:,:3] slip_direction = slipSystems[options.lattice][:,:3]
c_normal = slipSystems[options.lattice][:,3:] slip_normal = slipSystems[options.lattice][:,3:]
elif options.lattice == latticeChoices[2]: elif options.lattice == latticeChoices[2]:
# convert 4 Miller index notation of hex to orthogonal 3 Miller index notation # convert 4 Miller index notation of hex to orthogonal 3 Miller index notation
for i in range(len(c_direction)): for i in range(len(slip_direction)):
c_direction[i] = np.array([slipSystems['hex'][i,0]*1.5, slip_direction[i] = np.array([slipSystems['hex'][i,0]*1.5,
(slipSystems['hex'][i,0] + 2.*slipSystems['hex'][i,1])*0.5*np.sqrt(3), (slipSystems['hex'][i,0] + 2.*slipSystems['hex'][i,1])*0.5*np.sqrt(3),
slipSystems['hex'][i,3]*options.CoverA, slipSystems['hex'][i,3]*options.CoverA,
]) ])
c_normal[i] = np.array([slipSystems['hex'][i,4], slip_normal[i] = np.array([slipSystems['hex'][i,4],
(slipSystems['hex'][i,4] + 2.*slipSystems['hex'][i,5])/np.sqrt(3), (slipSystems['hex'][i,4] + 2.*slipSystems['hex'][i,5])/np.sqrt(3),
slipSystems['hex'][i,7]/options.CoverA, slipSystems['hex'][i,7]/options.CoverA,
]) ])
c_direction /= np.tile(np.linalg.norm(c_direction,axis=1),(3,1)).T slip_direction /= np.tile(np.linalg.norm(slip_direction,axis=1),(3,1)).T
c_normal /= np.tile(np.linalg.norm(c_normal ,axis=1),(3,1)).T slip_normal /= np.tile(np.linalg.norm(slip_normal ,axis=1),(3,1)).T
# --- loop over input files ------------------------------------------------------------------------ # --- loop over input files ------------------------------------------------------------------------
@ -244,7 +244,7 @@ for name in filenames:
.format( id = i+1, .format( id = i+1,
normal = theNormal, normal = theNormal,
direction = theDirection, direction = theDirection,
) for i,(theNormal,theDirection) in enumerate(zip(c_normal,c_direction))]) ) for i,(theNormal,theDirection) in enumerate(zip(slip_normal,slip_direction))])
table.head_write() table.head_write()
# ------------------------------------------ process data ------------------------------------------ # ------------------------------------------ process data ------------------------------------------
@ -262,9 +262,9 @@ for name in filenames:
elif inputtype == 'quaternion': 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]))),)
rotForce = o.quaternion.conjugated() * force
rotNormal = o.quaternion.conjugated() * normal table.data_append( np.abs( np.sum(slip_direction * (o.quaternion * force) ,axis=1) \
table.data_append(np.abs(np.sum(c_direction*rotForce,axis=1) * np.sum(c_normal*rotNormal,axis=1))) * np.sum(slip_normal * (o.quaternion * normal),axis=1)))
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------

2
processing/pre/seeds_fromRandom.py
View File

@ -127,7 +127,7 @@ options.fraction = np.array(options.fraction)
options.grid = np.array(options.grid) options.grid = np.array(options.grid)
gridSize = options.grid.prod() 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 np.random.seed(options.randomSeed) # init random generators
random.seed(options.randomSeed) random.seed(options.randomSeed)

220
src/IO.f90
View File

@ -33,7 +33,6 @@ module IO
IO_write_jobIntFile, & IO_write_jobIntFile, &
IO_read_realFile, & IO_read_realFile, &
IO_read_intFile, & IO_read_intFile, &
IO_hybridIA, &
IO_isBlank, & IO_isBlank, &
IO_getTag, & IO_getTag, &
IO_stringPos, & IO_stringPos, &
@ -583,223 +582,6 @@ logical function IO_abaqus_hasNoPart(fileUnit)
620 end function IO_abaqus_hasNoPart 620 end function IO_abaqus_hasNoPart
#endif #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 !> @brief identifies strings without content
@ -1758,7 +1540,6 @@ integer(pInt) function IO_verifyIntValue (string,validChars,myName)
validChars, & !< valid characters in string validChars, & !< valid characters in string
myName !< name of caller function (for debugging) myName !< name of caller function (for debugging)
integer(pInt) :: readStatus, invalidWhere integer(pInt) :: readStatus, invalidWhere
!character(len=len(trim(string))) :: trimmed does not work with ifort 14.0.1
IO_verifyIntValue = 0_pInt IO_verifyIntValue = 0_pInt
@ -1788,7 +1569,6 @@ real(pReal) function IO_verifyFloatValue (string,validChars,myName)
myName !< name of caller function (for debugging) myName !< name of caller function (for debugging)
integer(pInt) :: readStatus, invalidWhere integer(pInt) :: readStatus, invalidWhere
!character(len=len(trim(string))) :: trimmed does not work with ifort 14.0.1
IO_verifyFloatValue = 0.0_pReal IO_verifyFloatValue = 0.0_pReal

129
src/material.f90
View File

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