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

This commit is contained in:
Martin Diehl 2019-02-26 07:54:03 +01:00
parent b1bb68d523 f23eedb45e
commit 56e2c1264b
41 changed files with 5026 additions and 5098 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
.gitlab-ci.yml
View File

@ -207,7 +207,9 @@ Post_ParaviewRelated:
Post_OrientationConversion: Post_OrientationConversion:
stage: postprocessing stage: postprocessing
script: OrientationConversion/test.py script:
- OrientationConversion/test.py
- OrientationConversion/test2.py
except: except:
- master - master
- release - release
@ -508,7 +510,7 @@ Processing:
stage: createDocumentation stage: createDocumentation
script: script:
- cd $DAMASKROOT/processing/pre - cd $DAMASKROOT/processing/pre
- rm 3DRVEfrom2Dang.py abq_addUserOutput.py marc_addUserOutput.py - rm abq_addUserOutput.py marc_addUserOutput.py
- $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py - $DAMASKROOT/PRIVATE/documenting/scriptHelpToWiki.py --debug *.py
- cd $DAMASKROOT/processing/post - cd $DAMASKROOT/processing/post
- rm marc_to_vtk.py vtk2ang.py - rm marc_to_vtk.py vtk2ang.py

2
PRIVATE

@ -1 +1 @@
Subproject commit 75fbf8c1d9eb9b08fa15b55b7caaa4c4f7c167e0 Subproject commit def4081e837539dba7c4760abbb340553be66d3c

2
VERSION
View File

@ -1 +1 @@
v2.0.2-1837-g3bec76e7 v2.0.2-1937-ge401c212

1
examples/ConfigFiles/Crystallite_All.config
View File

@ -3,7 +3,6 @@
(output) texture (output) texture
(output) volume (output) volume
(output) orientation # quaternion (output) orientation # quaternion
(output) eulerangles # orientation as Bunge triple in degree
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates (output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) in crystal reference coordinates
(output) f # deformation gradient tensor (output) f # deformation gradient tensor
(output) fe # elastic deformation gradient tensor (output) fe # elastic deformation gradient tensor

1
examples/SpectralMethod/Polycrystal/material.config
View File

@ -13,7 +13,6 @@ mech none
(output) texture (output) texture
(output) volume (output) volume
(output) orientation # quaternion (output) orientation # quaternion
(output) eulerangles # orientation as Bunge triple
(output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4) (output) grainrotation # deviation from initial orientation as axis (1-3) and angle in degree (4)
(output) f # deformation gradient tensor; synonyms: "defgrad" (output) f # deformation gradient tensor; synonyms: "defgrad"
(output) fe # elastic deformation gradient tensor (output) fe # elastic deformation gradient tensor

84
processing/misc/DREAM3D_toTable.py
View File

@ -1,84 +0,0 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,h5py
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog [dream3dfile[s]]', description = """
Convert DREAM3D file to ASCIItable. Works for 3D datasets, but, hey, its not called DREAM2D ;)
""", version = scriptID)
parser.add_option('-d','--data',
dest = 'data',
action = 'extend', metavar = '<string LIST>',
help = 'data to extract from DREAM3D file')
parser.add_option('-c','--container',
dest = 'container', metavar = 'string',
help = 'root container(group) in which data is stored [%default]')
parser.set_defaults(container="ImageDataContainer",
)
(options, filenames) = parser.parse_args()
if options.data is None:
parser.error('No data selected')
rootDir ='DataContainers/'+options.container
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: parser.error('no input file specified.')
for name in filenames:
try:
table = damask.ASCIItable(outname = os.path.splitext(name)[0]+'.txt',
buffered = False
)
except: continue
damask.util.report(scriptName,name)
inFile = h5py.File(name, 'r')
try:
grid = inFile[rootDir+'/_SIMPL_GEOMETRY/DIMENSIONS'][...]
except:
damask.util.croak('Group {} not found'.format(options.container))
table.close(dismiss = True)
continue
# --- read comments --------------------------------------------------------------------------------
coords = (np.mgrid[0:grid[2], 0:grid[1], 0: grid[0]]).reshape(3, -1).T
table.data = (np.fliplr(coords)*inFile[rootDir+'/_SIMPL_GEOMETRY/SPACING'][...] \
+ inFile[rootDir+'/_SIMPL_GEOMETRY/ORIGIN'][...] \
+ inFile[rootDir+'/_SIMPL_GEOMETRY/SPACING'][...]*0.5)
labels = ['1_pos','2_pos','3_pos']
for data in options.data:
try:
l = np.prod(inFile[rootDir+'/CellData/'+data].shape[3:])
labels+=['{}_{}'.format(i+1,data.replace(' ','')) for i in range(l)] if l >1 else [data.replace(' ','')]
except KeyError:
damask.util.croak('Data {} not found'.format(data))
pass
table.data = np.hstack((table.data,
inFile[rootDir+'/CellData/'+data][...].reshape(grid.prod(),l)))
# ------------------------------------------ assemble header ---------------------------------------
table.labels_clear()
table.labels_append(labels,reset = True)
table.head_write()
# ------------------------------------------ finalize output ---------------------------------------
table.data_writeArray()
table.close()

48
processing/misc/OIMgrainFile_toTable.py
View File

@ -1,48 +0,0 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog [file[s]]', description = """
Adds header to OIM grain file type 1 to make it accesible as ASCII table
""", version = scriptID)
(options, filenames) = parser.parse_args()
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name,
buffered = False,
labeled = False)
except: continue
damask.util.report(scriptName,name)
table.head_read()
data = []
while table.data_read():
data.append(table.data[0:9])
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append(['1_euler','2_euler','3_euler','1_pos','2_pos','IQ','CI','Fit','GrainID'])
table.head_write()
for i in data:
table.data = i
table.data_write()
# --- output finalization --------------------------------------------------------------------------
table.close() # close ASCII table

176
processing/post/addGrainID.py
View File

@ -1,176 +0,0 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,copy
import numpy as np
import damask
from optparse import OptionParser
from scipy import spatial
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
Add grain index based on similiarity of crystal lattice orientation.
""", version = scriptID)
parser.add_option('-r',
'--radius',
dest = 'radius',
type = 'float', metavar = 'float',
help = 'search radius')
parser.add_option('-d',
'--disorientation',
dest = 'disorientation',
type = 'float', metavar = 'float',
help = 'disorientation threshold in degrees [%default]')
parser.add_option('-s',
'--symmetry',
dest = 'symmetry', type = 'choice', choices = damask.Symmetry.lattices[1:],
metavar = 'string',
help = 'crystal symmetry [%default] {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
parser.add_option('-o',
'--orientation',
dest = 'quaternion',
metavar = 'string',
help = 'label of crystal orientation given as unit quaternion [%default]')
parser.add_option('-p',
'--pos', '--position',
dest = 'pos',
metavar = 'string',
help = 'label of coordinates [%default]')
parser.add_option('--quiet',
dest='verbose',
action = 'store_false',
help = 'hide status bar (useful when piping to file)')
parser.set_defaults(disorientation = 5,
verbose = True,
quaternion = 'orientation',
symmetry = damask.Symmetry.lattices[-1],
pos = 'pos',
)
(options, filenames) = parser.parse_args()
if options.radius is None:
parser.error('no radius specified.')
cos_disorientation = np.cos(np.radians(options.disorientation/2.)) # cos of half the disorientation angle
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try: table = damask.ASCIItable(name = name,
buffered = False)
except: continue
damask.util.report(scriptName,name)
# ------------------------------------------ read header -------------------------------------------
table.head_read()
# ------------------------------------------ sanity checks -----------------------------------------
errors = []
remarks = []
if not 3 >= table.label_dimension(options.pos) >= 1:
errors.append('coordinates "{}" need to have one, two, or three dimensions.'.format(options.pos))
if not np.all(table.label_dimension(options.quaternion) == 4):
errors.append('input "{}" does not have dimension 4.'.format(options.quaternion))
else: column = table.label_index(options.quaternion)
if remarks != []: damask.util.croak(remarks)
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
# ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_append('grainID_{}@{:g}'.format(options.quaternion,options.disorientation)) # report orientation source and disorientation
table.head_write()
# ------------------------------------------ build KD tree -----------------------------------------
table.data_readArray(options.pos) # read position vectors
grainID = -np.ones(len(table.data),dtype=int)
Npoints = table.data.shape[0]
kdtree = spatial.KDTree(copy.deepcopy(table.data))
# ------------------------------------------ assign grain IDs --------------------------------------
orientations = [] # quaternions found for grain
memberCounts = [] # number of voxels in grain
p = 0 # point counter
g = 0 # grain counter
matchedID = -1
lastDistance = np.dot(kdtree.data[-1]-kdtree.data[0],kdtree.data[-1]-kdtree.data[0]) # (arbitrarily) use diagonal of cloud
table.data_rewind()
while table.data_read(): # read next data line of ASCII table
if options.verbose and Npoints > 100 and p%(Npoints//100) == 0: # report in 1% steps if possible and avoid modulo by zero
damask.util.progressBar(iteration=p,total=Npoints)
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced()
matched = False
alreadyChecked = {}
candidates = []
bestDisorientation = damask.Quaternion([0,0,0,1]) # initialize to 180 deg rotation as worst case
for i in kdtree.query_ball_point(kdtree.data[p],options.radius): # check all neighboring points
gID = grainID[i]
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.q > cos_disorientation: # within threshold ...
candidates.append(gID) # remember as potential candidate
if disorientation.quaternion.q >= bestDisorientation.q: # ... and better than current best?
matched = True
matchedID = gID # remember that grain
bestDisorientation = disorientation.quaternion
if matched: # did match existing grain
memberCounts[matchedID] += 1
if len(candidates) > 1: # ambiguity in grain identification?
largestGrain = sorted(candidates,key=lambda x:memberCounts[x])[-1] # find largest among potential candidate grains
matchedID = largestGrain
for c in [c for c in candidates if c != largestGrain]: # loop over smaller candidates
memberCounts[largestGrain] += memberCounts[c] # reassign member count of smaller to largest
memberCounts[c] = 0
grainID = np.where(np.in1d(grainID,candidates), largestGrain, grainID) # relabel grid points of smaller candidates as largest one
else: # no match -> new grain found
orientations += [o] # initialize with current orientation
memberCounts += [1] # start new membership counter
matchedID = g
g += 1 # increment grain counter
grainID[p] = matchedID # remember grain index assigned to point
p += 1 # increment point
grainIDs = np.where(np.array(memberCounts) > 0)[0] # identify "live" grain identifiers
packingMap = dict(zip(list(grainIDs),range(len(grainIDs)))) # map to condense into consecutive IDs
table.data_rewind()
outputAlive = True
p = 0
damask.util.progressBar(iteration=1,total=1)
while outputAlive and table.data_read(): # read next data line of ASCII table
table.data_append(1+packingMap[grainID[p]]) # add (condensed) grain ID
outputAlive = table.data_write() # output processed line
p += 1
# ------------------------------------------ output finalization -----------------------------------
table.close() # close ASCII tables

8
processing/post/addIPFcolor.py
View File

@ -41,6 +41,10 @@ parser.set_defaults(pole = (0.0,0.0,1.0),
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
# damask.Orientation requires Bravais lattice, but we are only interested in symmetry
symmetry2lattice={'cubic':'bcc','hexagonal':'hex','tetragonal':'bct'}
lattice = symmetry2lattice[options.symmetry]
pole = np.array(options.pole) pole = np.array(options.pole)
pole /= np.linalg.norm(pole) pole /= np.linalg.norm(pole)
@ -78,8 +82,8 @@ for name in filenames:
outputAlive = True outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), o = damask.Orientation(np.array(list(map(float,table.data[column:column+4]))),
symmetry = options.symmetry).reduced() lattice = lattice).reduced()
table.data_append(o.IPFcolor(pole)) table.data_append(o.IPFcolor(pole))
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line

5
processing/post/addMises.py
View File

@ -38,9 +38,12 @@ parser.add_option('-s','--stress',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'heading(s) of columns containing stress tensors') help = 'heading(s) of columns containing stress tensors')
parser.set_defaults(strain = [],
stress = [],
)
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
if options.stress is None and options.strain is None: if options.stress is [] and options.strain is []:
parser.error('no data column specified...') parser.error('no data column specified...')
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------

88
processing/post/addOrientations.py
View File

@ -9,31 +9,6 @@ import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# convention conformity checks
# --------------------------------------------------------------------
def check_Eulers(eulers):
if np.any(eulers < 0.0) or np.any(eulers > 2.0*np.pi) or eulers[1] > np.pi: # Euler angles within valid range?
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eulers))
return eulers
def check_quaternion(q):
if q[0] < 0.0: # positive first quaternion component?
raise ValueError('quaternion has negative first component.\n{}'.format(q[0]))
if not np.isclose(np.linalg.norm(q), 1.0): # unit quaternion?
raise ValueError('quaternion is not of unit length.\n{} {} {} {}'.format(*q))
return q
def check_matrix(M):
if not np.isclose(np.linalg.det(M),1.0): # proper rotation?
raise ValueError('matrix is not a proper rotation.\n{}'.format(M))
if not np.isclose(np.dot(M[0],M[1]), 0.0) \
or not np.isclose(np.dot(M[1],M[2]), 0.0) \
or not np.isclose(np.dot(M[2],M[0]), 0.0): # all orthogonal?
raise ValueError('matrix is not orthogonal.\n{}'.format(M))
return M
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# MAIN # MAIN
# -------------------------------------------------------------------- # --------------------------------------------------------------------
@ -46,19 +21,19 @@ Additional (globally fixed) rotations of the lab frame and/or crystal frame can
""", version = scriptID) """, version = scriptID)
outputChoices = { representations = {
'quaternion': ['quat',4], 'quaternion': ['quat',4], #ToDo: Use here Rowenhorst names (qu/ro/om/ax?)
'rodrigues': ['rodr',3], 'rodrigues': ['rodr',4],
'eulers': ['eulr',3], 'eulers': ['eulr',3],
'matrix': ['mtrx',9], 'matrix': ['mtrx',9],
'angleaxis': ['aaxs',4], 'angleaxis': ['aaxs',4],
} }
parser.add_option('-o', parser.add_option('-o',
'--output', '--output',
dest = 'output', dest = 'output',
action = 'extend', metavar = '<string LIST>', action = 'extend', metavar = '<string LIST>',
help = 'output orientation formats {{{}}}'.format(', '.join(outputChoices))) help = 'output orientation formats {{{}}}'.format(', '.join(representations)))
parser.add_option('-d', parser.add_option('-d',
'--degrees', '--degrees',
dest = 'degrees', dest = 'degrees',
@ -104,15 +79,15 @@ parser.add_option('-z',
help = 'label of lab z vector (expressed in crystal coords)') help = 'label of lab z vector (expressed in crystal coords)')
parser.set_defaults(output = [], parser.set_defaults(output = [],
labrotation = (0.,1.,1.,1.), # no rotation about 1,1,1 labrotation = (0.,1.,0.,0.), # no rotation about 1,0,0
crystalrotation = (0.,1.,1.,1.), # no rotation about 1,1,1 crystalrotation = (0.,1.,0.,0.), # no rotation about 1,0,0
) )
(options, filenames) = parser.parse_args() (options, filenames) = parser.parse_args()
options.output = list(map(lambda x: x.lower(), options.output)) options.output = list(map(lambda x: x.lower(), options.output))
if options.output == [] or (not set(options.output).issubset(set(outputChoices))): if options.output == [] or (not set(options.output).issubset(set(representations))):
parser.error('output must be chosen from {}.'.format(', '.join(outputChoices))) parser.error('output must be chosen from {}.'.format(', '.join(representations)))
input = [options.eulers is not None, input = [options.eulers is not None,
options.rodrigues is not None, options.rodrigues is not None,
@ -125,16 +100,18 @@ input = [options.eulers is not None,
if np.sum(input) != 1: parser.error('needs exactly one input format.') if np.sum(input) != 1: parser.error('needs exactly one input format.')
(label,dim,inputtype) = [(options.eulers,3,'eulers'), (label,dim,inputtype) = [(options.eulers,representations['eulers'][1],'eulers'),
(options.rodrigues,3,'rodrigues'), (options.rodrigues,representations['rodrigues'][1],'rodrigues'),
([options.x,options.y,options.z],[3,3,3],'frame'), ([options.x,options.y,options.z],[3,3,3],'frame'),
(options.matrix,9,'matrix'), (options.matrix,representations['matrix'][1],'matrix'),
(options.quaternion,4,'quaternion'), (options.quaternion,representations['quaternion'][1],'quaternion'),
][np.where(input)[0][0]] # select input label that was requested ][np.where(input)[0][0]] # select input label that was requested
toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians crystalrotation = np.array(options.crystalrotation[1:4] + (options.crystalrotation[0],)) # Compatibility hack
r = damask.Quaternion.fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation labrotation = np.array(options.labrotation[1:4] + (options.labrotation[0],)) # Compatibility hack
R = damask.Quaternion.fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation r = damask.Rotation.fromAxisAngle(crystalrotation,options.degrees) # crystal frame rotation
R = damask.Rotation.fromAxisAngle(labrotation,options.degrees) # lab frame rotation
# --- loop over input files ------------------------------------------------------------------------ # --- loop over input files ------------------------------------------------------------------------
@ -168,9 +145,9 @@ for name in filenames:
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:])) table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
for output in options.output: for output in options.output:
if output in outputChoices: if output in representations:
table.labels_append(['{}_{}({})'.format(i+1,outputChoices[output][0],label) \ table.labels_append(['{}_{}({})'.format(i+1,representations[output][0],label) \
for i in range(outputChoices[output][1])]) for i in range(representations[output][1])])
table.head_write() table.head_write()
# ------------------------------------------ process data ------------------------------------------ # ------------------------------------------ process data ------------------------------------------
@ -178,30 +155,35 @@ for name in filenames:
outputAlive = True outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers': if inputtype == 'eulers':
l = representations['eulers'][1]
o = damask.Rotation.fromEulers(list(map(float,table.data[column:column+l])),options.degrees)
o = damask.Orientation(Eulers = check_Eulers(np.array(list(map(float,table.data[column:column+3])))*toRadians))
elif inputtype == 'rodrigues': elif inputtype == 'rodrigues':
o = damask.Orientation(Rodrigues = np.array(list(map(float,table.data[column:column+3])))) l = representations['rodrigues'][1]
elif inputtype == 'matrix': o = damask.Rotation.fromRodrigues(list(map(float,table.data[column:column+l])))
o = damask.Orientation(matrix = check_matrix(np.array(list(map(float,table.data[column:column+9]))).reshape(3,3))) elif inputtype == 'matrix':
l = representations['matrix'][1]
o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+l])))
elif inputtype == 'frame': elif inputtype == 'frame':
M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \ table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3).T table.data[column[2]:column[2]+3]))).reshape(3,3).T
o = damask.Orientation(matrix = check_matrix(M/np.linalg.norm(M,axis=0))) o = damask.Rotation.fromMatrix(M/np.linalg.norm(M,axis=0))
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = check_quaternion(np.array(list(map(float,table.data[column:column+4]))))) elif inputtype == 'quaternion':
l = representations['quaternion'][1]
o = damask.Rotation.fromQuaternion(list(map(float,table.data[column:column+l])))
o.quaternion = r*o.quaternion*R # apply additional lab and crystal frame rotations o= r*o*R # apply additional lab and crystal frame rotations
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(degrees=options.degrees)) elif output == 'eulers': table.data_append(o.asEulers(degrees=options.degrees))
elif output == 'matrix': table.data_append(o.asMatrix()) elif output == 'matrix': table.data_append(o.asMatrix())
elif output == 'angleaxis': table.data_append(o.asAngleAxis(degrees=options.degrees,flat=True)) elif output == 'angleaxis': table.data_append(o.asAxisAngle(degrees=options.degrees))
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------

4
processing/post/addPole.py
View File

@ -75,9 +75,9 @@ for name in filenames:
# ------------------------------------------ process data ------------------------------------------ # ------------------------------------------ process data ------------------------------------------
outputAlive = True outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4])))) o = damask.Rotation(np.array(list(map(float,table.data[column:column+4]))))
rotatedPole = o.quaternion*pole # rotate pole according to crystal orientation rotatedPole = o*pole # rotate pole according to crystal orientation
(x,y) = rotatedPole[0:2]/(1.+abs(pole[2])) # stereographic projection (x,y) = rotatedPole[0:2]/(1.+abs(pole[2])) # stereographic projection
table.data_append([np.sqrt(x*x+y*y),np.arctan2(y,x)] if options.polar else [x,y]) # cartesian coordinates table.data_append([np.sqrt(x*x+y*y),np.arctan2(y,x)] if options.polar else [x,y]) # cartesian coordinates

6
processing/post/addSchmidfactors.py
View File

@ -212,10 +212,10 @@ for name in filenames:
outputAlive = True outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4])))) o = damask.Rotation(list(map(float,table.data[column:column+4])))
table.data_append( np.abs( np.sum(slip_direction * (o.quaternion * force) ,axis=1) \ table.data_append( np.abs( np.sum(slip_direction * (o * force) ,axis=1) \
* np.sum(slip_normal * (o.quaternion * normal),axis=1))) * np.sum(slip_normal * (o * normal),axis=1)))
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------

14
processing/post/rotateData.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys,math import os,sys
import numpy as np import numpy as np
from optparse import OptionParser from optparse import OptionParser
import damask import damask
@ -40,9 +40,8 @@ parser.set_defaults(rotation = (0.,1.,1.,1.),
if options.data is None: if options.data is None:
parser.error('no data column specified.') parser.error('no data column specified.')
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians rotation = np.array(options.rotation[1:4]+(options.rotation[0],)) # Compatibility hack
q = damask.Quaternion().fromAngleAxis(toRadians*options.rotation[0],options.rotation[1:]) r = damask.Rotation.fromAxisAngle(rotation,options.degrees,normalise=True)
R = q.asMatrix()
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------
@ -90,12 +89,11 @@ for name in filenames:
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
for v in active['vector']: for v in active['vector']:
column = table.label_index(v) column = table.label_index(v)
table.data[column:column+3] = q * np.array(list(map(float,table.data[column:column+3]))) table.data[column:column+3] = r * np.array(list(map(float,table.data[column:column+3])))
for t in active['tensor']: for t in active['tensor']:
column = table.label_index(t) column = table.label_index(t)
table.data[column:column+9] = \ table.data[column:column+9] = (r * np.array(list(map(float,table.data[column:column+9]))).reshape((3,3))).reshape(9)
np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+9]))).reshape((3,3)),
R.transpose())).reshape((9))
outputAlive = table.data_write() # output processed line outputAlive = table.data_write() # output processed line
# ------------------------------------------ output finalization ----------------------------------- # ------------------------------------------ output finalization -----------------------------------

11
processing/post/vtk_addGridData.py
View File

@ -53,19 +53,22 @@ parser.set_defaults(data = [],
if not options.vtk: parser.error('No VTK file specified.') if not options.vtk: parser.error('No VTK file specified.')
if not os.path.exists(options.vtk): parser.error('VTK file does not exist.') if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
if os.path.splitext(options.vtk)[1] == '.vtr': vtk_file,vtk_ext = os.path.splitext(options.vtk)
if vtk_ext == '.vtr':
reader = vtk.vtkXMLRectilinearGridReader() reader = vtk.vtkXMLRectilinearGridReader()
reader.SetFileName(options.vtk) reader.SetFileName(options.vtk)
reader.Update() reader.Update()
rGrid = reader.GetOutput() rGrid = reader.GetOutput()
writer = vtk.vtkXMLRectilinearGridWriter() writer = vtk.vtkXMLRectilinearGridWriter()
elif os.path.splitext(options.vtk)[1] == '.vtk': elif vtk_ext == '.vtk':
reader = vtk.vtkGenericDataObjectReader() reader = vtk.vtkGenericDataObjectReader()
reader.SetFileName(options.vtk) reader.SetFileName(options.vtk)
reader.Update() reader.Update()
rGrid = reader.GetRectilinearGridOutput() rGrid = reader.GetRectilinearGridOutput()
writer = vtk.vtkXMLRectilinearGridWriter() writer = vtk.vtkXMLRectilinearGridWriter()
elif os.path.splitext(options.vtk)[1] == '.vtu': vtk_ext = '.vtr'
elif vtk_ext == '.vtu':
reader = vtk.vtkXMLUnstructuredGridReader() reader = vtk.vtkXMLUnstructuredGridReader()
reader.SetFileName(options.vtk) reader.SetFileName(options.vtk)
reader.Update() reader.Update()
@ -74,7 +77,7 @@ elif os.path.splitext(options.vtk)[1] == '.vtu':
else: else:
parser.error('Unsupported VTK file type extension.') parser.error('Unsupported VTK file type extension.')
writer.SetFileName(options.vtk) writer.SetFileName(vtk_file+vtk_ext)
Npoints = rGrid.GetNumberOfPoints() Npoints = rGrid.GetNumberOfPoints()
Ncells = rGrid.GetNumberOfCells() Ncells = rGrid.GetNumberOfCells()

9
processing/post/vtk_addPointCloudData.py
View File

@ -49,16 +49,19 @@ parser.set_defaults(data = [],
if not options.vtk: parser.error('no VTK file specified.') if not options.vtk: parser.error('no VTK file specified.')
if not os.path.exists(options.vtk): parser.error('VTK file does not exist.') if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
if os.path.splitext(options.vtk)[1] == '.vtp': vtk_file,vtk_ext = os.path.splitext(options.vtk)
if vtk_ext == '.vtp':
reader = vtk.vtkXMLPolyDataReader() reader = vtk.vtkXMLPolyDataReader()
reader.SetFileName(options.vtk) reader.SetFileName(options.vtk)
reader.Update() reader.Update()
Polydata = reader.GetOutput() Polydata = reader.GetOutput()
elif os.path.splitext(options.vtk)[1] == '.vtk': elif vtk_ext == '.vtk':
reader = vtk.vtkGenericDataObjectReader() reader = vtk.vtkGenericDataObjectReader()
reader.SetFileName(options.vtk) reader.SetFileName(options.vtk)
reader.Update() reader.Update()
Polydata = reader.GetPolyDataOutput() Polydata = reader.GetPolyDataOutput()
vtk_ext = '.vtp'
else: else:
parser.error('unsupported VTK file type extension.') parser.error('unsupported VTK file type extension.')
@ -149,7 +152,7 @@ for name in filenames:
writer = vtk.vtkXMLPolyDataWriter() writer = vtk.vtkXMLPolyDataWriter()
writer.SetDataModeToBinary() writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib() writer.SetCompressorTypeToZLib()
writer.SetFileName(options.vtk) writer.SetFileName(vtk_file+vtk_ext)
writer.SetInputData(Polydata) writer.SetInputData(Polydata)
writer.Write() writer.Write()

8
processing/post/vtk_addRectilinearGridData.py
View File

@ -53,16 +53,18 @@ parser.set_defaults(data = [],
if not options.vtk: parser.error('no VTK file specified.') if not options.vtk: parser.error('no VTK file specified.')
if not os.path.exists(options.vtk): parser.error('VTK file does not exist.') if not os.path.exists(options.vtk): parser.error('VTK file does not exist.')
if os.path.splitext(options.vtk)[1] == '.vtr': vtk_file,vtk_ext = os.path.splitext(options.vtk)
if vtk_ext == '.vtr':
reader = vtk.vtkXMLRectilinearGridReader() reader = vtk.vtkXMLRectilinearGridReader()
reader.SetFileName(options.vtk) reader.SetFileName(options.vtk)
reader.Update() reader.Update()
rGrid = reader.GetOutput() rGrid = reader.GetOutput()
elif os.path.splitext(options.vtk)[1] == '.vtk': elif vtk_ext == '.vtk':
reader = vtk.vtkGenericDataObjectReader() reader = vtk.vtkGenericDataObjectReader()
reader.SetFileName(options.vtk) reader.SetFileName(options.vtk)
reader.Update() reader.Update()
rGrid = reader.GetRectilinearGridOutput() rGrid = reader.GetRectilinearGridOutput()
vtk_ext = '.vtr'
else: else:
parser.error('unsupported VTK file type extension.') parser.error('unsupported VTK file type extension.')
@ -159,7 +161,7 @@ for name in filenames:
writer = vtk.vtkXMLRectilinearGridWriter() writer = vtk.vtkXMLRectilinearGridWriter()
writer.SetDataModeToBinary() writer.SetDataModeToBinary()
writer.SetCompressorTypeToZLib() writer.SetCompressorTypeToZLib()
writer.SetFileName(options.vtk) writer.SetFileName(vtk_file+vtk_ext)
writer.SetInputData(rGrid) writer.SetInputData(rGrid)
writer.Write() writer.Write()

119
processing/pre/3DRVEfrom2Dang.py
View File

@ -1,119 +0,0 @@
#!/usr/bin/env python2.7
# -*- coding: UTF-8 no BOM -*-
import os,sys,math
from optparse import OptionParser
import damask
import pipes
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]',
description ='generate 3D RVE from .ang files of EBSD slices .',
version = scriptID)
parser.add_option('--offset',
dest='offset',
type='float',
help='offset of EBSD slices [%default]',
metavar='float')
parser.add_option('--outname',
dest='outName',
type='string',
help='output file name [%default]', metavar='string')
parser.add_option('--vtr',
action="store_true",
dest='vtr')
parser.add_option('--geom',
action="store_true",
dest='geom')
parser.set_defaults(offset = 1.0,
outName = 'RVE3D')
(options,filenames) = parser.parse_args()
numFiles = len(filenames)
formatwidth = 1+int(math.log10(numFiles))
# copy original files to tmp files to not alter originals
for i in range(numFiles):
sliceID = 'slice' + str(i).zfill(formatwidth) + '.tmp'
strCommand = 'cp ' + pipes.quote(filenames[i]) + ' ' + sliceID
os.system(strCommand)
# modify tmp files
print('Add z-coordinates')
for i in range(numFiles):
sliceID = 'slice' + str(i).zfill(formatwidth) + '.tmp'
strCommand = 'OIMgrainFile_toTable ' + sliceID
os.system(strCommand)
strCommand = 'addCalculation --label 3Dpos --formula "np.array(#pos#.tolist()+[' + str(i*options.offset) + '])" ' + sliceID
os.system(strCommand)
# join temp files into one
print('\n Colocate files')
fileOut = open(options.outName + '.ang','w')
# take header information from 1st file
sliceID = 'slice' + str(0).zfill(formatwidth) + '.tmp'
fileRead = open(sliceID)
data = fileRead.readlines()
fileRead.close()
headerLines = int(data[0].split()[0])
fileOut.write(str(headerLines+1) + '\t header\n')
for line in data[1:headerLines]:
fileOut.write(line)
fileOut.write(scriptID + '\t' + ' '.join(sys.argv[1:]) + '\n')
for line in data[headerLines:]:
fileOut.write(line)
# append other files content without header
for i in range(numFiles-1):
sliceID = 'slice' + str(i+1).zfill(formatwidth) + '.tmp'
fileRead = open(sliceID)
data = fileRead.readlines()
fileRead.close()
headerLines = int(data[0].split()[0])
for line in data[headerLines+1:]:
fileOut.write(line)
fileOut.close()
# tidy up and add phase column
print('\n Remove temp data and add phase info')
strCommand = 'filterTable --black pos ' + options.outName + '.ang'
os.system(strCommand)
strCommand = 'reLabel --label 3Dpos --substitute pos ' + options.outName + '.ang'
os.system(strCommand)
strCommand = 'addCalculation -l phase -f 1 ' + options.outName + '.ang'
os.system(strCommand)
# create geom file when asked for
if options.geom:
print('\n Build geometry file')
strCommand = 'geom_fromTable --phase phase --eulers euler --coordinates pos ' + pipes.quote(options.outName) + '.ang'
os.system(strCommand)
# create paraview file when asked for
if options.vtr:
print('\n Build Paraview file')
strCommand = 'addIPFcolor --eulers euler --pole 0.0 0.0 1.0 ' + options.outName + '.ang'
os.system(strCommand)
strCommand = 'vtk_rectilinearGrid ' + pipes.quote(options.outName) + '.ang'
os.system(strCommand)
os.rename(pipes.quote(options.outName) + '_pos(cell)'+'.vtr', pipes.quote(options.outName) + '.vtr')
strCommand = 'vtk_addRectilinearGridData --vtk '+ pipes.quote(options.outName) + '.vtr --color IPF_001_cubic '\
+ pipes.quote(options.outName) + '.ang'
os.system(strCommand)
# delete tmp files
for i in range(numFiles):
sliceID = 'slice' + str(i).zfill(formatwidth) + '.tmp'
os.remove(sliceID)

15
processing/pre/geom_addPrimitive.py
View File

@ -43,7 +43,7 @@ parser.add_option('-f', '--fill', dest='fill', type='int', metavar = 'int'
help='grain index to fill primitive. "0" selects maximum microstructure index + 1 [%default]') help='grain index to fill primitive. "0" selects maximum microstructure index + 1 [%default]')
parser.add_option('-q', '--quaternion', dest='quaternion', type='float', nargs = 4, metavar=' '.join(['float']*4), parser.add_option('-q', '--quaternion', dest='quaternion', type='float', nargs = 4, metavar=' '.join(['float']*4),
help = 'rotation of primitive as quaternion') help = 'rotation of primitive as quaternion')
parser.add_option('-a', '--angleaxis', dest='angleaxis', nargs = 4, metavar=' '.join(['float']*4), parser.add_option('-a', '--angleaxis', dest='angleaxis', nargs = 4, metavar=' '.join(['float']*4), type=float,
help = 'angle,x,y,z clockwise rotation of primitive about axis by angle') help = 'angle,x,y,z clockwise rotation of primitive about axis by angle')
parser.add_option( '--degrees', dest='degrees', action='store_true', parser.add_option( '--degrees', dest='degrees', action='store_true',
help = 'angle is given in degrees [%default]') help = 'angle is given in degrees [%default]')
@ -63,14 +63,12 @@ parser.set_defaults(center = (.0,.0,.0),
if options.dimension is None: if options.dimension is None:
parser.error('no dimension specified.') parser.error('no dimension specified.')
if options.angleaxis is not None: if options.angleaxis is not None:
options.angleaxis = list(map(float,options.angleaxis)) ax = np.array(options.angleaxis[1:4] + (options.angleaxis[0],)) # Compatibility hack
rotation = damask.Quaternion.fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0], rotation = damask.Rotation.fromAxisAngle(ax,options.degrees,normalise=True)
options.angleaxis[1:4])
elif options.quaternion is not None: elif options.quaternion is not None:
options.quaternion = list(map(float,options.quaternion)) rotation = damask.Rotation.fromQuaternion(options.quaternion)
rotation = damask.Quaternion(quat=options.quaternion)
else: else:
rotation = damask.Quaternion() rotation = damask.Rotation()
options.center = np.array(options.center) options.center = np.array(options.center)
options.dimension = np.array(options.dimension) options.dimension = np.array(options.dimension)
@ -159,8 +157,7 @@ for name in filenames:
X -= options.center[0] - 0.5 X -= options.center[0] - 0.5
Y -= options.center[1] - 0.5 Y -= options.center[1] - 0.5
Z -= options.center[2] - 0.5 Z -= options.center[2] - 0.5
# and then by applying the quaternion # and then by applying the rotation
# this should be rotation.conjugate() * (X,Y,Z), but it is this way for backwards compatibility with the older version of this script
(X, Y, Z) = rotation * (X, Y, Z) (X, Y, Z) = rotation * (X, Y, Z)
# and finally by scaling (we don't worry about options.dimension being negative, np.abs occurs on the microstructure = np.where... line) # and finally by scaling (we don't worry about options.dimension being negative, np.abs occurs on the microstructure = np.where... line)
X /= options.dimension[0] * 0.5 X /= options.dimension[0] * 0.5

4
processing/pre/geom_check.sh
View File

@ -18,8 +18,8 @@ do
< $geom \ < $geom \
| \ | \
vtk_addRectilinearGridData \ vtk_addRectilinearGridData \
--vtk ${geom%.*}.vtk \
--data microstructure \ --data microstructure \
--inplace \
--vtk ${geom%.*}.vtk
rm ${geom%.*}.vtk rm ${geom%.*}.vtk
done done

189
processing/pre/geom_fromDREAM3D.py Executable file
View File

@ -0,0 +1,189 @@
#!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,h5py
import numpy as np
from optparse import OptionParser
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
#--------------------------------------------------------------------------------------------------
# MAIN
#--------------------------------------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog [dream3dfile[s]]', description = """
Convert DREAM3D file to geometry file. This can be done from cell data (direct pointwise takeover) or
from grain data (individual grains are segmented). Requires orientation data as quaternion.
""", version = scriptID)
parser.add_option('-b','--basegroup',
dest = 'basegroup', metavar = 'string',
help = 'name of the group in "DataContainers" that contains all the data')
parser.add_option('-p','--pointwise',
dest = 'pointwise', metavar = 'string',
help = 'name of the group in "DataContainers/<basegroup>" that contains pointwise data [%default]')
parser.add_option('-a','--average',
dest = 'average', metavar = 'string',
help = 'name of the group in "DataContainers</basegroup>" that contains grain average data. '\
+ 'Leave empty for pointwise data')
parser.add_option('--phase',
dest = 'phase',
type = 'string', metavar = 'string',
help = 'name of the dataset containing pointwise/average phase IDs [%default]')
parser.add_option('--microstructure',
dest = 'microstructure',
type = 'string', metavar = 'string',
help = 'name of the dataset connecting pointwise and average data [%default]')
parser.add_option('-q', '--quaternion',
dest = 'quaternion',
type = 'string', metavar='string',
help = 'name of the dataset containing pointwise/average orientation as quaternion [%default]')
parser.set_defaults(pointwise = 'CellData',
quaternion = 'Quats',
phase = 'Phases',
microstructure = 'FeatureIds',
crystallite = 1,
)
(options, filenames) = parser.parse_args()
if options.basegroup is None:
parser.error('No base group selected')
rootDir ='DataContainers'
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: parser.error('no input file specified.')
for name in filenames:
try:
table = damask.ASCIItable(outname = os.path.splitext(name)[0]+'.geom',
buffered = False, labeled=False,
)
except: continue
damask.util.report(scriptName,name)
errors = []
info = {}
ori = []
inFile = h5py.File(name, 'r')
group_geom = os.path.join(rootDir,options.basegroup,'_SIMPL_GEOMETRY')
try:
info['size'] = inFile[os.path.join(group_geom,'DIMENSIONS')][...] \
* inFile[os.path.join(group_geom,'SPACING')][...]
info['grid'] = inFile[os.path.join(group_geom,'DIMENSIONS')][...]
info['origin'] = inFile[os.path.join(group_geom,'ORIGIN')][...]
except:
errors.append('Geometry data ({}) not found'.format(group_geom))
group_pointwise = os.path.join(rootDir,options.basegroup,options.pointwise)
if options.average is None:
label = 'point'
N_microstructure = np.product(info['grid'])
dataset = os.path.join(group_pointwise,options.quaternion)
try:
quats = np.reshape(inFile[dataset][...],(N_microstructure,3))
except:
errors.append('Pointwise orientation data ({}) not found'.format(dataset))
texture = [damask.Rotation.fromQuaternion(q,True,P=+1) for q in quats]
dataset = os.path.join(group_pointwise,options.phase)
try:
phase = np.reshape(inFile[dataset][...],(N_microstructure))
except:
errors.append('Pointwise phase data ({}) not found'.format(dataset))
else:
label = 'grain'
dataset = os.path.join(group_pointwise,options.microstructure)
try:
microstructure = np.reshape(inFile[dataset][...],(np.product(info['grid'])))
N_microstructure = np.max(microstructure)
except:
errors.append('Link between pointwise and grain average data ({}) not found'.format(dataset))
group_average = os.path.join(rootDir,options.basegroup,options.average)
dataset = os.path.join(group_average,options.quaternion)
try:
texture = [damask.Rotation.fromQuaternion(q,True,P=+1) for q in inFile[dataset][...][1:]] # skip first entry (unindexed)
except:
errors.append('Average orientation data ({}) not found'.format(dataset))
dataset = os.path.join(group_average,options.phase)
try:
phase = [i[0] for i in inFile[dataset][...]][1:] # skip first entry (unindexed)
except:
errors.append('Average phase data ({}) not found'.format(dataset))
if errors != []:
damask.util.croak(errors)
table.close(dismiss = True)
continue
mat = damask.Material()
mat.verbose = False
# dummy <homogenization>
h = damask.config.material.Homogenization()
mat.add_section('Homogenization','none',h)
info['homogenization'] = 1
# <crystallite> placeholder (same for all microstructures at the moment)
c = damask.config.material.Crystallite()
mat.add_section('Crystallite','tbd',c)
# <phase> placeholders
for i in range(np.max(phase)):
p = damask.config.material.Phase()
mat.add_section('phase','phase{}-tbd'.format(i+1),p)
# <texture>
for i,o in enumerate(texture):
t = damask.config.material.Texture()
t.add_component('gauss',{'eulers':o.asEulers(degrees=True)})
mat.add_section(part='texture', section='{}{}'.format(label,i+1),initialData=t)
# <microstructure>
for i in range(N_microstructure):
m = damask.config.material.Microstructure()
mat.add_section('microstructure','{}{}'.format(label,i+1),m)
mat.add_microstructure('{}{}'.format(label,i+1),
{'phase': 'phase{}-tbd'.format(phase[i]),
'texture':'{}{}'.format(label,i+1),
'crystallite':'tbd',
'fraction':1
})
table.info_append([
scriptID + ' ' + ' '.join(sys.argv[1:]),
"grid\ta {}\tb {}\tc {}".format(*info['grid']),
"size\tx {}\ty {}\tz {}".format(*info['size']),
"origin\tx {}\ty {}\tz {}".format(*info['origin']),
"homogenization\t{}".format(info['homogenization']),
str(mat).split('\n')
])
table.head_write()
if options.average is None:
table.data = [1, 'to', format(N_microstructure)]
table.data_write()
else:
table.data = microstructure.reshape(info['grid'][1]*info['grid'][2],info['grid'][0])
table.data_writeArray()
table.close()

155
processing/pre/geom_fromTable.py
View File

@ -1,8 +1,8 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*- # -*- coding: UTF-8 no BOM -*-
import os,sys,math,time import os,sys,math
import scipy.spatial, numpy as np import numpy as np
from optparse import OptionParser from optparse import OptionParser
import damask import damask
@ -32,34 +32,6 @@ parser.add_option('--microstructure',
dest = 'microstructure', dest = 'microstructure',
type = 'string', metavar = 'string', type = 'string', metavar = 'string',
help = 'microstructure label') help = 'microstructure label')
parser.add_option('-t', '--tolerance',
dest = 'tolerance',
type = 'float', metavar = 'float',
help = 'angular tolerance for orientation squashing [%default]')
parser.add_option('-e', '--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'Euler angles label')
parser.add_option('-d', '--degrees',
dest = 'degrees',
action = 'store_true',
help = 'all angles are in degrees')
parser.add_option('-m', '--matrix',
dest = 'matrix',
type = 'string', metavar = 'string',
help = 'orientation matrix label')
parser.add_option('-a',
dest='a',
type = 'string', metavar = 'string',
help = 'crystal frame a vector label')
parser.add_option('-b',
dest='b',
type = 'string', metavar = 'string',
help = 'crystal frame b vector label')
parser.add_option('-c',
dest = 'c',
type = 'string', metavar='string',
help = 'crystal frame c vector label')
parser.add_option('-q', '--quaternion', parser.add_option('-q', '--quaternion',
dest = 'quaternion', dest = 'quaternion',
type = 'string', metavar='string', type = 'string', metavar='string',
@ -67,11 +39,8 @@ parser.add_option('-q', '--quaternion',
parser.add_option('--axes', parser.add_option('--axes',
dest = 'axes', dest = 'axes',
type = 'string', nargs = 3, metavar = ' '.join(['string']*3), type = 'string', nargs = 3, metavar = ' '.join(['string']*3),
help = 'orientation coordinate frame in terms of position coordinate frame [same]') help = 'orientation coordinate frame in terms of position coordinate frame [+x +y +z]')
parser.add_option('-s', '--symmetry',
dest = 'symmetry',
action = 'extend', metavar = '<string LIST>',
help = 'crystal symmetry of each phase %default {{{}}} '.format(', '.join(damask.Symmetry.lattices[1:])))
parser.add_option('--homogenization', parser.add_option('--homogenization',
dest = 'homogenization', dest = 'homogenization',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
@ -80,27 +49,16 @@ parser.add_option('--crystallite',
dest = 'crystallite', dest = 'crystallite',
type = 'int', metavar = 'int', type = 'int', metavar = 'int',
help = 'crystallite index to be used [%default]') help = 'crystallite index to be used [%default]')
parser.add_option('--verbose',
dest = 'verbose', action = 'store_true',
help = 'output extra info')
parser.set_defaults(symmetry = [damask.Symmetry.lattices[-1]],
tolerance = 0.0, parser.set_defaults(homogenization = 1,
degrees = False,
homogenization = 1,
crystallite = 1, crystallite = 1,
verbose = False,
pos = 'pos', pos = 'pos',
) )
(options,filenames) = parser.parse_args() (options,filenames) = parser.parse_args()
input = [options.eulers is not None, input = [ options.quaternion is not None,
options.a is not None and \
options.b is not None and \
options.c is not None,
options.matrix is not None,
options.quaternion is not None,
options.microstructure is not None, options.microstructure is not None,
] ]
@ -109,14 +67,9 @@ if np.sum(input) != 1:
if options.axes is not None and not set(options.axes).issubset(set(['x','+x','-x','y','+y','-y','z','+z','-z'])): if options.axes is not None and not set(options.axes).issubset(set(['x','+x','-x','y','+y','-y','z','+z','-z'])):
parser.error('invalid axes {} {} {}.'.format(*options.axes)) parser.error('invalid axes {} {} {}.'.format(*options.axes))
(label,dim,inputtype) = [(options.eulers,3,'eulers'), (label,dim,inputtype) = [(options.quaternion,4,'quaternion'),
([options.a,options.b,options.c],[3,3,3],'frame'),
(options.matrix,9,'matrix'),
(options.quaternion,4,'quaternion'),
(options.microstructure,1,'microstructure'), (options.microstructure,1,'microstructure'),
][np.where(input)[0][0]] # select input label that was requested ][np.where(input)[0][0]] # select input label that was requested
toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale all angles to radians
threshold = np.cos(options.tolerance/2.*toRadians) # cosine of (half of) tolerance angle
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------
@ -146,7 +99,7 @@ for name in filenames:
if options.phase and table.label_dimension(options.phase) != 1: if options.phase and table.label_dimension(options.phase) != 1:
errors.append('phase column "{}" is not scalar.'.format(options.phase)) errors.append('phase column "{}" is not scalar.'.format(options.phase))
if errors != []: if errors != []:
damask.util.croak(errors) damask.util.croak(errors)
table.close(dismiss = True) table.close(dismiss = True)
continue continue
@ -157,10 +110,8 @@ for name in filenames:
if coordDim == 2: if coordDim == 2:
table.data = np.insert(table.data,2,np.zeros(len(table.data)),axis=1) # add zero z coordinate for two-dimensional input table.data = np.insert(table.data,2,np.zeros(len(table.data)),axis=1) # add zero z coordinate for two-dimensional input
if options.verbose: damask.util.croak('extending to 3D...')
if options.phase is None: if options.phase is None:
table.data = np.column_stack((table.data,np.ones(len(table.data)))) # add single phase if no phase column given table.data = np.column_stack((table.data,np.ones(len(table.data)))) # add single phase if no phase column given
if options.verbose: damask.util.croak('adding dummy phase info...')
# --------------- figure out size and grid --------------------------------------------------------- # --------------- figure out size and grid ---------------------------------------------------------
@ -196,17 +147,10 @@ for name in filenames:
grain = table.data[:,colOri] grain = table.data[:,colOri]
nGrains = len(np.unique(grain)) nGrains = len(np.unique(grain))
else: elif inputtype == 'quaternion':
if options.verbose: bg = damask.util.backgroundMessage(); bg.start() # start background messaging
colPhase = -1 # column of phase data comes last colPhase = -1 # column of phase data comes last
if options.verbose: bg.set_message('sorting positions...')
index = np.lexsort((table.data[:,0],table.data[:,1],table.data[:,2])) # index of position when sorting x fast, z slow index = np.lexsort((table.data[:,0],table.data[:,1],table.data[:,2])) # index of position when sorting x fast, z slow
if options.verbose: bg.set_message('building KD tree...')
KDTree = scipy.spatial.KDTree((table.data[index,:3]-mincorner) / delta) # build KDTree with dX = dY = dZ = 1 and origin 0,0,0
statistics = {'global': 0, 'local': 0}
grain = -np.ones(N,dtype = 'int32') # initialize empty microstructure grain = -np.ones(N,dtype = 'int32') # initialize empty microstructure
orientations = [] # orientations orientations = [] # orientations
multiplicity = [] # orientation multiplicity (number of group members) multiplicity = [] # orientation multiplicity (number of group members)
@ -215,87 +159,26 @@ for name in filenames:
existingGrains = np.arange(nGrains) existingGrains = np.arange(nGrains)
myPos = 0 # position (in list) of current grid point myPos = 0 # position (in list) of current grid point
tick = time.clock()
if options.verbose: bg.set_message('assigning grain IDs...')
for z in range(grid[2]): for z in range(grid[2]):
for y in range(grid[1]): for y in range(grid[1]):
for x in range(grid[0]): for x in range(grid[0]):
if (myPos+1)%(N/500.) < 1:
time_delta = (time.clock()-tick) * (N - myPos) / myPos
if options.verbose: bg.set_message('(%02i:%02i:%02i) processing point %i of %i (grain count %i)...'
%(time_delta//3600,time_delta%3600//60,time_delta%60,myPos,N,nGrains))
myData = table.data[index[myPos]] # read data for current grid point myData = table.data[index[myPos]] # read data for current grid point
myPhase = int(myData[colPhase]) myPhase = int(myData[colPhase])
mySym = options.symmetry[min(myPhase,len(options.symmetry))-1] # take last specified option for all with higher index
o = damask.Rotation(myData[colOri:colOri+4])
if inputtype == 'eulers':
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),
symmetry = mySym)
elif inputtype == 'frame':
o = damask.Orientation(matrix = np.hstack((myData[colOri[0]:colOri[0]+3],
myData[colOri[1]:colOri[1]+3],
myData[colOri[2]:colOri[2]+3],
)).reshape(3,3),
symmetry = mySym)
elif inputtype == 'quaternion':
o = damask.Orientation(quaternion = myData[colOri:colOri+4],
symmetry = mySym)
cos_disorientations = -np.ones(1,dtype=float) # largest possible disorientation grain[myPos] = nGrains # assign new grain to me ...
closest_grain = -1 # invalid neighbor nGrains += 1 # ... and update counter
orientations.append(o) # store new orientation for future comparison
if options.tolerance > 0.0: # only try to compress orientations if asked to multiplicity.append(1) # having single occurrence so far
neighbors = np.array(KDTree.query_ball_point([x,y,z], 3)) # point indices within radius phases.append(myPhase) # store phase info for future reporting
# filter neighbors: skip myself, anyone further ahead (cannot yet have a grain ID), and other phases existingGrains = np.arange(nGrains) # update list of existing grains
neighbors = neighbors[(neighbors < myPos) & \
(table.data[index[neighbors],colPhase] == myPhase)]
grains = np.unique(grain[neighbors]) # unique grain IDs among valid neighbors
if len(grains) > 0: # check immediate neighborhood first
cos_disorientations = np.array([o.disorientation(orientations[grainID],
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'
if cos_disorientations[closest_grain] < threshold: # orientation not close enough?
grains = existingGrains[np.atleast_1d( (np.array(phases) == myPhase ) & \
(np.in1d(existingGrains,grains,invert=True)))] # other already identified grains (of my phase)
if len(grains) > 0:
cos_disorientations = np.array([o.disorientation(orientations[grainID],
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'
if cos_disorientations[closest_grain] >= threshold: # orientation now close enough?
grainID = grains[closest_grain]
grain[myPos] = grainID # assign myself to that grain ...
orientations[grainID] = damask.Orientation.average([orientations[grainID],o],
[multiplicity[grainID],1]) # update average orientation of best matching grain
multiplicity[grainID] += 1
statistics[match] += 1
else:
grain[myPos] = nGrains # assign new grain to me ...
nGrains += 1 # ... and update counter
orientations.append(o) # store new orientation for future comparison
multiplicity.append(1) # having single occurrence so far
phases.append(myPhase) # store phase info for future reporting
existingGrains = np.arange(nGrains) # update list of existing grains
myPos += 1 myPos += 1
if options.verbose:
bg.stop()
bg.join()
damask.util.croak("{} seconds total.\n{} local and {} global matches.".\
format(time.clock()-tick,statistics['local'],statistics['global']))
grain += 1 # offset from starting index 0 to 1 grain += 1 # offset from starting index 0 to 1

15
processing/pre/geom_rotate.py
View File

@ -52,13 +52,14 @@ parser.set_defaults(degrees = False,
if sum(x is not None for x in [options.rotation,options.eulers,options.matrix,options.quaternion]) != 1: if sum(x is not None for x in [options.rotation,options.eulers,options.matrix,options.quaternion]) != 1:
parser.error('not exactly one rotation specified...') parser.error('not exactly one rotation specified...')
eulers = np.array(damask.orientation.Orientation( if options.quaternion is not None:
quaternion = np.array(options.quaternion) if options.quaternion else None, eulers = damask.Rotation.fromQuaternion(np.array(options.quaternion)).asEulers(degrees=True)
angleAxis = np.array(options.rotation) if options.rotation else None, if options.rotation is not None:
matrix = np.array(options.matrix) if options.matrix else None, eulers = damask.Rotation.fromAxisAngle(np.array(options.rotation,degrees=True)).asEulers(degrees=True)
Eulers = np.array(options.eulers) if options.eulers else None, if options.matrix is not None:
degrees = options.degrees, eulers = damask.Rotation.fromMatrix(np.array(options.Matrix)).asEulers(degrees=True)
).asEulers(degrees=True)) if options.eulers is not None:
eulers = damask.Rotation.fromEulers(np.array(options.eulers),degrees=True).asEulers(degrees=True)
# --- loop over input files ------------------------------------------------------------------------- # --- loop over input files -------------------------------------------------------------------------

14
processing/pre/seeds_fromRandom.py
View File

@ -90,11 +90,7 @@ group.add_option( '-s',
'--selective', '--selective',
action = 'store_true', action = 'store_true',
dest = 'selective', dest = 'selective',
help = 'selective picking of seed points from random seed points [%default]') help = 'selective picking of seed points from random seed points')
group.add_option( '--force',
action = 'store_true',
dest = 'force',
help = 'try selective picking despite large seed point number [%default]')
group.add_option( '--distance', group.add_option( '--distance',
dest = 'distance', dest = 'distance',
type = 'float', metavar = 'float', type = 'float', metavar = 'float',
@ -115,7 +111,6 @@ parser.set_defaults(randomSeed = None,
sigma = 0.05, sigma = 0.05,
microstructure = 1, microstructure = 1,
selective = False, selective = False,
force = False,
distance = 0.2, distance = 0.2,
numCandidates = 10, numCandidates = 10,
format = None, format = None,
@ -148,10 +143,11 @@ for name in filenames:
errors = [] errors = []
if gridSize == 0: if gridSize == 0:
errors.append('zero grid dimension for {}.'.format(', '.join([['a','b','c'][x] for x in np.where(options.grid == 0)[0]]))) errors.append('zero grid dimension for {}.'.format(', '.join([['a','b','c'][x] for x in np.where(options.grid == 0)[0]])))
if options.N > gridSize/10.: errors.append('seed count exceeds 0.1 of grid points.') if options.N > gridSize/10.:
remarks.append('seed count exceeds 0.1 of grid points.')
if options.selective and 4./3.*math.pi*(options.distance/2.)**3*options.N > 0.5: if options.selective and 4./3.*math.pi*(options.distance/2.)**3*options.N > 0.5:
(remarks if options.force else errors).append('maximum recommended seed point count for given distance is {}.{}'. remarks.append('maximum recommended seed point count for given distance is {}.{}'.
format(int(3./8./math.pi/(options.distance/2.)**3),'..'*options.force)) format(int(3./8./math.pi/(options.distance/2.)**3)))
if remarks != []: damask.util.croak(remarks) if remarks != []: damask.util.croak(remarks)
if errors != []: if errors != []:

125
python/damask/Lambert.py Normal file
View File

@ -0,0 +1,125 @@
# -*- coding: UTF-8 no BOM -*-
####################################################################################################
# Code below available according to the followin conditions on https://github.com/MarDiehl/3Drotations
####################################################################################################
# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification, are
# permitted provided that the following conditions are met:
#
# - Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
# - Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or
# other materials provided with the distribution.
# - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
# of its contributors may be used to endorse or promote products derived from
# this software without specific prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
# USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
####################################################################################################
import numpy as np
sc = np.pi**(1./6.)/6.**(1./6.)
beta = np.pi**(5./6.)/6.**(1./6.)/2.
R1 = (3.*np.pi/4.)**(1./3.)
def CubeToBall(cube):
if np.abs(np.max(cube))>np.pi**(2./3.) * 0.5:
raise ValueError
# transform to the sphere grid via the curved square, and intercept the zero point
if np.allclose(cube,0.0,rtol=0.0,atol=1.0e-300):
ball = np.zeros(3)
else:
# get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cube)
XYZ = cube[p] * sc
# intercept all the points along the z-axis
if np.allclose(XYZ[0:2],0.0,rtol=0.0,atol=1.0e-300):
ball = np.array([0.0, 0.0, np.sqrt(6.0/np.pi) * XYZ[2]])
else:
order = [1,0] if np.abs(XYZ[1]) <= np.abs(XYZ[0]) else [0,1]
q = np.pi/12.0 * XYZ[order[0]]/XYZ[order[1]]
c = np.cos(q)
s = np.sin(q)
q = R1*2.0**0.25/beta * XYZ[order[1]] / np.sqrt(np.sqrt(2.0)-c)
T = np.array([ (np.sqrt(2.0)*c - 1.0), np.sqrt(2.0) * s]) * q
# transform to sphere grid (inverse Lambert)
# note that there is no need to worry about dividing by zero, since XYZ[2] can not become zero
c = np.sum(T**2)
s = c * np.pi/24.0 /XYZ[2]**2
c = c * np.sqrt(np.pi/24.0)/XYZ[2]
q = np.sqrt( 1.0 - s )
ball = np.array([ T[order[1]] * q, T[order[0]] * q, np.sqrt(6.0/np.pi) * XYZ[2] - c ])
# reverse the coordinates back to the regular order according to the original pyramid number
ball = ball[p]
return ball
def BallToCube(ball):
rs = np.linalg.norm(ball)
if rs > R1:
raise ValueError
if np.allclose(ball,0.0,rtol=0.0,atol=1.0e-300):
cube = np.zeros(3)
else:
p = GetPyramidOrder(ball)
xyz3 = ball[p]
# inverse M_3
xyz2 = xyz3[0:2] * np.sqrt( 2.0*rs/(rs+np.abs(xyz3[2])) )
# inverse M_2
qxy = np.sum(xyz2**2)
if np.isclose(qxy,0.0,rtol=0.0,atol=1.0e-300):
Tinv = np.zeros(2)
else:
q2 = qxy + np.max(np.abs(xyz2))**2
sq2 = np.sqrt(q2)
q = (beta/np.sqrt(2.0)/R1) * np.sqrt(q2*qxy/(q2-np.max(np.abs(xyz2))*sq2))
tt = np.clip((np.min(np.abs(xyz2))**2+np.max(np.abs(xyz2))*sq2)/np.sqrt(2.0)/qxy,-1.0,1.0)
Tinv = np.array([1.0,np.arccos(tt)/np.pi*12.0]) if np.abs(xyz2[1]) <= np.abs(xyz2[0]) else \
np.array([np.arccos(tt)/np.pi*12.0,1.0])
Tinv = q * np.where(xyz2<0.0,-Tinv,Tinv)
# inverse M_1
cube = np.array([ Tinv[0], Tinv[1], (-1.0 if xyz3[2] < 0.0 else 1.0) * rs / np.sqrt(6.0/np.pi) ]) /sc
# reverse the coordinates back to the regular order according to the original pyramid number
cube = cube[p]
return cube
def GetPyramidOrder(xyz):
if (abs(xyz[0])<= xyz[2]) and (abs(xyz[1])<= xyz[2]) or \
(abs(xyz[0])<=-xyz[2]) and (abs(xyz[1])<=-xyz[2]):
return [0,1,2]
elif (abs(xyz[2])<= xyz[0]) and (abs(xyz[1])<= xyz[0]) or \
(abs(xyz[2])<=-xyz[0]) and (abs(xyz[1])<=-xyz[0]):
return [1,2,0]
elif (abs(xyz[0])<= xyz[1]) and (abs(xyz[2])<= xyz[1]) or \
(abs(xyz[0])<=-xyz[1]) and (abs(xyz[2])<=-xyz[1]):
return [2,0,1]

2
python/damask/__init__.py
View File

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

23
python/damask/config/material.py
View File

@ -77,18 +77,6 @@ class Texture(Section):
) )
) )
if multiKey == 'fiber':
self.add_multiKey(multiKey,'alpha1 %g\talpha2 %g\tbeta1 %g\tbeta2 %g\tscatter %g\tfraction %g'%(
properties['eulers'][0],
properties['eulers'][1],
properties['eulers'][2],
properties['eulers'][3],
scatter,
fraction,
)
)
class Material(): class Material():
"""Reads, manipulates and writes material.config files""" """Reads, manipulates and writes material.config files"""
@ -97,10 +85,10 @@ class Material():
"""Generates ordered list of parts""" """Generates ordered list of parts"""
self.parts = [ self.parts = [
'homogenization', 'homogenization',
'microstructure',
'crystallite', 'crystallite',
'phase', 'phase',
'texture', 'texture',
'microstructure',
] ]
self.data = {\ self.data = {\
'homogenization': {'__order__': []}, 'homogenization': {'__order__': []},
@ -117,15 +105,12 @@ class Material():
for part in self.parts: for part in self.parts:
if self.verbose: print('processing <{}>'.format(part)) if self.verbose: print('processing <{}>'.format(part))
me += ['', me += ['',
'#-----------------------------#', '#'*100,
'<{}>'.format(part), '<{}>'.format(part),
'#-----------------------------#', '#'*100,
] ]
for section in self.data[part]['__order__']: for section in self.data[part]['__order__']:
me += ['', me += ['[{}] {}'.format(section,'#'+'-'*max(0,96-len(section)))]
'[{}] {}'.format(section,'#'*max(0,27-len(section))),
'',
]
for key in self.data[part][section]['__order__']: for key in self.data[part][section]['__order__']:
if key.startswith('(') and key.endswith(')'): # multiple (key) if key.startswith('(') and key.endswith(')'): # multiple (key)
me += ['{}\t{}'.format(key,' '.join(values)) for values in self.data[part][section][key]] me += ['{}\t{}'.format(key,' '.join(values)) for values in self.data[part][section][key]]

1917
python/damask/orientation.py
View File

File diff suppressed because it is too large Load Diff

16
src/CMakeLists.txt
View File

@ -62,6 +62,18 @@ add_library(MATH OBJECT "math.f90")
add_dependencies(MATH NUMERICS) add_dependencies(MATH NUMERICS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATH>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATH>)
add_library(QUATERNIONS OBJECT "quaternions.f90")
add_dependencies(QUATERNIONS MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:QUATERNIONS>)
add_library(LAMBERT OBJECT "Lambert.f90")
add_dependencies(LAMBERT MATH)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:LAMBERT>)
add_library(ROTATIONS OBJECT "rotations.f90")
add_dependencies(ROTATIONS LAMBERT QUATERNIONS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:ROTATIONS>)
add_library(MESH_BASE OBJECT "mesh_base.f90") add_library(MESH_BASE OBJECT "mesh_base.f90")
add_dependencies(MESH_BASE ELEMENT) add_dependencies(MESH_BASE ELEMENT)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH_BASE>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH_BASE>)
@ -75,13 +87,13 @@ elseif (PROJECT_NAME STREQUAL "DAMASK_FEM")
add_library(FEZoo OBJECT "FEM_zoo.f90") add_library(FEZoo OBJECT "FEM_zoo.f90")
add_dependencies(FEZoo IO) add_dependencies(FEZoo IO)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:FEZoo>)
add_library(MESH OBJECT "mesh_FEM.f90") add_library(MESH OBJECT "mesh_FEM.f90")
add_dependencies(MESH FEZoo MESH_BASE MATH FEsolving) add_dependencies(MESH FEZoo MESH_BASE MATH FEsolving)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:MESH>)
endif() endif()
add_library(MATERIAL OBJECT "material.f90") add_library(MATERIAL OBJECT "material.f90")
add_dependencies(MATERIAL MESH DAMASK_CONFIG) add_dependencies(MATERIAL MESH DAMASK_CONFIG ROTATIONS)
list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATERIAL>) list(APPEND OBJECTFILES $<TARGET_OBJECTS:MATERIAL>)
add_library(LATTICE OBJECT "lattice.f90") add_library(LATTICE OBJECT "lattice.f90")

217
src/Lambert.f90 Normal file
View File

@ -0,0 +1,217 @@
! ###################################################################
! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University
! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
! - Redistributions of source code must retain the above copyright notice, this list
! of conditions and the following disclaimer.
! - Redistributions in binary form must reproduce the above copyright notice, this
! list of conditions and the following disclaimer in the documentation and/or
! other materials provided with the distribution.
! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
! of its contributors may be used to endorse or promote products derived from
! this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Mapping homochoric <-> cubochoric
!
!> @details
!> D. Rosca, A. Morawiec, and M. De Graef. “A new method of constructing a grid
!> in the space of 3D rotations and its applications to texture analysis”.
!> Modeling and Simulations in Materials Science and Engineering 22, 075013 (2014).
!--------------------------------------------------------------------------
module Lambert
use math
use prec, only: &
pReal
implicit none
private
real(pReal), parameter, private :: &
SPI = sqrt(PI), &
PREF = sqrt(6.0_pReal/PI), &
A = PI**(5.0_pReal/6.0_pReal)/6.0_pReal**(1.0_pReal/6.0_pReal), &
AP = PI**(2.0_pReal/3.0_pReal), &
SC = A/AP, &
BETA = A/2.0_pReal, &
R1 = (3.0_pReal*PI/4.0_pReal)**(1.0_pReal/3.0_pReal), &
R2 = sqrt(2.0_pReal), &
PI12 = PI/12.0_pReal, &
PREK = R1 * 2.0_pReal**(1.0_pReal/4.0_pReal)/BETA
public :: &
LambertCubeToBall, &
LambertBallToCube
private :: &
GetPyramidOrder
contains
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D cubic grid to 3D ball
!--------------------------------------------------------------------------
function LambertCubeToBall(cube) result(ball)
use, intrinsic :: IEEE_ARITHMETIC
use prec, only: &
pInt, &
dEq0
implicit none
real(pReal), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ
real(pReal) :: T(2), c, s, q
real(pReal), parameter :: eps = 1.0e-8_pReal
integer(pInt), dimension(3) :: p
integer(pInt), dimension(2) :: order
if (maxval(abs(cube)) > AP/2.0+eps) then
ball = IEEE_value(cube,IEEE_positive_inf)
return
end if
! transform to the sphere grid via the curved square, and intercept the zero point
center: if (all(dEq0(cube))) then
ball = 0.0_pReal
else center
! get pyramide and scale by grid parameter ratio
p = GetPyramidOrder(cube)
XYZ = cube(p) * sc
! intercept all the points along the z-axis
special: if (all(dEq0(XYZ(1:2)))) then
LamXYZ = [ 0.0_pReal, 0.0_pReal, pref * XYZ(3) ]
else special
order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ
q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger
c = cos(q)
s = sin(q)
q = prek * XYZ(order(2))/ sqrt(R2-c)
T = [ (R2*c - 1.0), R2 * s] * q
! transform to sphere grid (inverse Lambert)
! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero]
c = sum(T**2)
s = Pi * c/(24.0*XYZ(3)**2)
c = sPi * c / sqrt(24.0_pReal) / XYZ(3)
q = sqrt( 1.0 - s )
LamXYZ = [ T(order(2)) * q, T(order(1)) * q, pref * XYZ(3) - c ]
endif special
! reverse the coordinates back to the regular order according to the original pyramid number
ball = LamXYZ(p)
endif center
end function LambertCubeToBall
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D ball to 3D cubic grid
!--------------------------------------------------------------------------
pure function LambertBallToCube(xyz) result(cube)
use, intrinsic :: IEEE_ARITHMETIC, only:&
IEEE_positive_inf, &
IEEE_value
use prec, only: &
pInt, &
dEq0
implicit none
real(pReal), intent(in), dimension(3) :: xyz
real(pReal), dimension(3) :: cube, xyz1, xyz3
real(pReal), dimension(2) :: Tinv, xyz2
real(pReal) :: rs, qxy, q2, sq2, q, tt
integer(pInt), dimension(3) :: p
rs = norm2(xyz)
if (rs > R1) then
cube = IEEE_value(cube,IEEE_positive_inf)
return
endif
center: if (all(dEq0(xyz))) then
cube = 0.0_pReal
else center
p = GetPyramidOrder(xyz)
xyz3 = xyz(p)
! inverse M_3
xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) )
! inverse M_2
qxy = sum(xyz2**2)
special: if (dEq0(qxy)) then
Tinv = 0.0
else special
q2 = qxy + maxval(abs(xyz2))**2
sq2 = sqrt(q2)
q = (beta/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2))
tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy
Tinv = q * sign(1.0,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], &
[ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], &
abs(xyz2(2)) <= abs(xyz2(1)))
endif special
! inverse M_1
xyz1 = [ Tinv(1), Tinv(2), sign(1.0,xyz3(3)) * rs / pref ] /sc
! reverst the coordinates back to the regular order according to the original pyramid number
cube = xyz1(p)
endif center
end function LambertBallToCube
!--------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief determine to which pyramid a point in a cubic grid belongs
!--------------------------------------------------------------------------
pure function GetPyramidOrder(xyz)
use prec, only: &
pInt
implicit none
real(pReal),intent(in),dimension(3) :: xyz
integer(pInt), dimension(3) :: GetPyramidOrder
if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. &
((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then
GetPyramidOrder = [1,2,3]
else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. &
((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then
GetPyramidOrder = [2,3,1]
else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. &
((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then
GetPyramidOrder = [3,1,2]
else
GetPyramidOrder = -1 ! should be impossible, but might simplify debugging
end if
end function GetPyramidOrder
end module Lambert

3
src/commercialFEM_fileList.f90
View File

@ -11,6 +11,9 @@
#include "HDF5_utilities.f90" #include "HDF5_utilities.f90"
#endif #endif
#include "math.f90" #include "math.f90"
#include "quaternions.f90"
#include "Lambert.f90"
#include "rotations.f90"
#include "FEsolving.f90" #include "FEsolving.f90"
#include "element.f90" #include "element.f90"
#include "mesh_base.f90" #include "mesh_base.f90"

132
src/constitutive.f90
View File

@ -48,11 +48,8 @@ subroutine constitutive_init()
use IO, only: & use IO, only: &
IO_error, & IO_error, &
IO_open_file, & IO_open_file, &
IO_checkAndRewind, &
IO_open_jobFile_stat, & IO_open_jobFile_stat, &
IO_write_jobFile IO_write_jobFile
use config, only: &
config_phase
use config, only: & use config, only: &
material_Nphase, & material_Nphase, &
material_localFileExt, & material_localFileExt, &
@ -132,40 +129,29 @@ subroutine constitutive_init()
nonlocalConstitutionPresent = .false. nonlocalConstitutionPresent = .false.
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! open material.config ! initialized plasticity
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file
!--------------------------------------------------------------------------------------------------
! parse plasticities from config file
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init
if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init
if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) call plastic_nonlocal_init
call plastic_nonlocal_init(FILEUNIT)
call plastic_nonlocal_stateInit()
endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parse source mechanisms from config file ! initialize source mechanisms
call IO_checkAndRewind(FILEUNIT)
if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init
if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init
if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init
if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init
if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init
if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! parse kinematic mechanisms from config file ! initialize kinematic mechanisms
call IO_checkAndRewind(FILEUNIT)
if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init
if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init
if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init
close(FILEUNIT)
call config_deallocate('material.config/phase') call config_deallocate('material.config/phase')
@ -337,7 +323,7 @@ end function constitutive_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different constitutive models !> @brief calls microstructure function of the different constitutive models
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) subroutine constitutive_microstructure(Fe, Fp, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use material, only: & use material, only: &
@ -352,7 +338,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
PLASTICITY_disloucla_ID, & PLASTICITY_disloucla_ID, &
PLASTICITY_nonlocal_ID PLASTICITY_nonlocal_ID
use plastic_nonlocal, only: & use plastic_nonlocal, only: &
plastic_nonlocal_microstructure plastic_nonlocal_dependentState
use plastic_dislotwin, only: & use plastic_dislotwin, only: &
plastic_dislotwin_dependentState plastic_dislotwin_dependentState
use plastic_disloUCLA, only: & use plastic_disloUCLA, only: &
@ -370,8 +356,6 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
instance, of instance, of
real(pReal), intent(in), dimension(:,:,:,:) :: &
orientations !< crystal orientations as quaternions
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
@ -386,7 +370,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_disloUCLA_dependentState(instance,of) call plastic_disloUCLA_dependentState(instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_microstructure (Fe,Fp,ip,el) call plastic_nonlocal_dependentState (Fe,Fp,ip,el)
end select plasticityType end select plasticityType
end subroutine constitutive_microstructure end subroutine constitutive_microstructure
@ -394,15 +378,15 @@ end subroutine constitutive_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: Discuss wheter it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, el) subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33
math_6toSym33, &
math_sym33to6, &
math_99to3333
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -418,6 +402,8 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_DISLOUCLA_ID, & PLASTICITY_DISLOUCLA_ID, &
PLASTICITY_NONLOCAL_ID PLASTICITY_NONLOCAL_ID
use mesh, only: &
mesh_ipVolume
use plastic_isotropic, only: & use plastic_isotropic, only: &
plastic_isotropic_LpAndItsTangent plastic_isotropic_LpAndItsTangent
use plastic_phenopowerlaw, only: & use plastic_phenopowerlaw, only: &
@ -436,9 +422,8 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola-Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
@ -447,11 +432,8 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
dLp_dFi !< derivative of Lp with respect to Fi dLp_dFi !< derivative of Lp with respect to Fi
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(9,9) :: &
dLp_dMp99 !< derivative of Lp with respect to Mstar (matrix notation)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mp, & !< Mandel stress work conjugate with Lp Mp !< Mandel stress work conjugate with Lp
S !< 2nd Piola-Kirchhoff stress
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme !< thermal member position tme !< thermal member position
@ -461,7 +443,6 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
S = math_6toSym33(S6)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -486,9 +467,8 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_sym33to6(Mp), & call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, &
temperature(ho)%p(tme),ip,el) temperature(ho)%p(tme),mesh_ipVolume(ip,el),ip,el)
dLp_dMp = math_99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
@ -523,15 +503,15 @@ end subroutine constitutive_LpAndItsTangents
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi? ! ToDo: MD: S is Mi?
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el) subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
S, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_I3, & math_I3, &
math_inv33, & math_inv33, &
math_det33, & math_det33, &
math_mul33x33, & math_mul33x33
math_6toSym33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -558,8 +538,8 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(3,3) :: &
S6 !< 2nd Piola-Kirchhoff stress (vector notation) S !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
@ -588,7 +568,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
case (PLASTICITY_isotropic_ID) plasticityType case (PLASTICITY_isotropic_ID) plasticityType
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6),instance,of) call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of)
case default plasticityType case default plasticityType
my_Li = 0.0_pReal my_Li = 0.0_pReal
my_dLi_dS = 0.0_pReal my_dLi_dS = 0.0_pReal
@ -600,9 +580,9 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, e
KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el))
kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6), ipc, ip, el) call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, math_6toSym33(S6), ipc, ip, el) call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ipc, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el) call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
case default kinematicsType case default kinematicsType
@ -634,9 +614,7 @@ pure function constitutive_initialFi(ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_I3, & math_I3
math_inv33, &
math_mul33x33
use material, only: & use material, only: &
material_phase, & material_phase, &
material_homog, & material_homog, &
@ -710,7 +688,8 @@ end subroutine constitutive_SandItsTangents
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermeidate deformation gradients using Hookes law !> the elastic and intermeidate deformation gradients using Hookes law
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only : & use math, only : &
@ -774,7 +753,7 @@ end subroutine constitutive_hooke_SandItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @brief contains the constitutive equation for calculating the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfracArray,ipc, ip, el) subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip, el)
use prec, only: & use prec, only: &
pReal, & pReal, &
pLongInt pLongInt
@ -784,8 +763,6 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_mul33x33, & math_mul33x33, &
math_6toSym33, &
math_sym33to6, &
math_mul33x33 math_mul33x33
use mesh, only: & use mesh, only: &
theMesh theMesh
@ -839,27 +816,25 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
el !< element el !< element
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), intent(in), dimension(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
subfracArray !< subfraction of timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
FeArray, & !< elastic deformation gradient FeArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(3,3) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation) S !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mp Mp
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
s, & !< counter in source loop i, & !< counter in source loop
instance, of instance, of
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_6toSym33(S6)) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -889,16 +864,16 @@ subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfrac
call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of) call plastic_disloucla_dotState (Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (math_sym33to6(Mp),FeArray,FpArray,temperature(ho)%p(tme), & call plastic_nonlocal_dotState (Mp,FeArray,FpArray,temperature(ho)%p(tme), &
subdt,subfracArray,ip,el) subdt,ip,el)
end select plasticityType end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) SourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
case (SOURCE_damage_anisoBrittle_ID) sourceType case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState (math_6toSym33(S6), ipc, ip, el) !< correct stress? call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState ( ipc, ip, el) call source_damage_isoDuctile_dotState ( ipc, ip, el)
@ -928,7 +903,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: & use math, only: &
math_sym33to6, &
math_mul33x33 math_mul33x33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
@ -972,7 +946,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
call plastic_kinehardening_deltaState(Mp,instance,of) call plastic_kinehardening_deltaState(Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(math_sym33to6(Mp),ip,el) call plastic_nonlocal_deltaState(Mp,ip,el)
end select plasticityType end select plasticityType
@ -994,14 +968,11 @@ end subroutine constitutive_collectDeltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results !> @brief returns array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el) function constitutive_postResults(S, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
math_6toSym33, &
math_mul33x33 math_mul33x33
use mesh, only: &
theMesh
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticityInstance, & phase_plasticityInstance, &
@ -1014,7 +985,6 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
material_homogenizationAt, & material_homogenizationAt, &
temperature, & temperature, &
thermalMapping, & thermalMapping, &
homogenization_maxNgrains, &
PLASTICITY_NONE_ID, & PLASTICITY_NONE_ID, &
PLASTICITY_ISOTROPIC_ID, & PLASTICITY_ISOTROPIC_ID, &
PLASTICITY_PHENOPOWERLAW_ID, & PLASTICITY_PHENOPOWERLAW_ID, &
@ -1057,10 +1027,8 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
constitutive_postResults constitutive_postResults
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: & real(pReal), intent(in), dimension(3,3) :: &
FeArray !< elastic deformation gradient S !< 2nd Piola Kirchhoff stress
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Mp !< Mandel stress Mp !< Mandel stress
integer(pInt) :: & integer(pInt) :: &
@ -1068,11 +1036,11 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
s, of, instance !< counter in source loop i, of, instance !< counter in source loop
constitutive_postResults = 0.0_pReal constitutive_postResults = 0.0_pReal
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_6toSym33(S6)) Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
@ -1113,14 +1081,14 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_nonlocal_postResults (S6,FeArray,ip,el) plastic_nonlocal_postResults (Mp,ip,el)
end select plasticityType end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el)) SourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
startPos = endPos + 1_pInt startPos = endPos + 1_pInt
endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(s)%sizePostResults endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(i)%sizePostResults
of = phasememberAt(ipc,ip,el) of = phasememberAt(ipc,ip,el)
sourceType: select case (phase_source(s,material_phase(ipc,ip,el))) sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
case (SOURCE_damage_isoBrittle_ID) sourceType case (SOURCE_damage_isoBrittle_ID) sourceType
constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(material_phase(ipc,ip,el),of) constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(material_phase(ipc,ip,el),of)
case (SOURCE_damage_isoDuctile_ID) sourceType case (SOURCE_damage_isoDuctile_ID) sourceType

1191
src/constitutive.f90.orig
View File

File diff suppressed because it is too large Load Diff

84
src/crystallite.f90
View File

@ -9,14 +9,16 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module crystallite module crystallite
use prec, only: &
pReal, &
pInt
use rotations, only: &
rotation
use FEsolving, only: & use FEsolving, only: &
FEsolving_execElem, & FEsolving_execElem, &
FEsolving_execIP FEsolving_execIP
use material, only: & use material, only: &
homogenization_Ngrains homogenization_Ngrains
use prec, only: &
pReal, &
pInt
implicit none implicit none
@ -40,10 +42,9 @@ module crystallite
crystallite_Tstar_v, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) ToDo: Should be called S, 3x3 crystallite_Tstar_v, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step) ToDo: Should be called S, 3x3
crystallite_Tstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc ToDo: Should be called S, 3x3 crystallite_Tstar0_v, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc ToDo: Should be called S, 3x3
crystallite_partionedTstar0_v !< 2nd Piola-Kirchhoff stress vector at start of homog inc ToDo: Should be called S, 3x3 crystallite_partionedTstar0_v !< 2nd Piola-Kirchhoff stress vector at start of homog inc ToDo: Should be called S, 3x3
real(pReal), dimension(:,:,:,:), allocatable, private :: & type(rotation), dimension(:,:,:), allocatable, private :: &
crystallite_orientation, & !< orientation as quaternion crystallite_orientation, & !< orientation
crystallite_orientation0, & !< initial orientation as quaternion crystallite_orientation0 !< initial orientation
crystallite_rotation !< grain rotation away from initial orientation as axis-angle (in degrees) in crystal reference frame
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
crystallite_Fe, & !< current "elastic" def grad (end of converged time step) crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
crystallite_P !< 1st Piola-Kirchhoff stress per grain crystallite_P !< 1st Piola-Kirchhoff stress per grain
@ -89,7 +90,6 @@ module crystallite
volume_ID, & volume_ID, &
orientation_ID, & orientation_ID, &
grainrotation_ID, & grainrotation_ID, &
eulerangles_ID, &
defgrad_ID, & defgrad_ID, &
fe_ID, & fe_ID, &
fp_ID, & fp_ID, &
@ -233,9 +233,8 @@ subroutine crystallite_init
allocate(crystallite_subdt(cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subdt(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subFrac(cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subFrac(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_subStep(cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_subStep(cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_orientation(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_orientation(cMax,iMax,eMax))
allocate(crystallite_orientation0(4,cMax,iMax,eMax), source=0.0_pReal) allocate(crystallite_orientation0(cMax,iMax,eMax))
allocate(crystallite_rotation(4,cMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.) allocate(crystallite_localPlasticity(cMax,iMax,eMax), source=.true.)
allocate(crystallite_requested(cMax,iMax,eMax), source=.false.) allocate(crystallite_requested(cMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(cMax,iMax,eMax), source=.false.) allocate(crystallite_todo(cMax,iMax,eMax), source=.false.)
@ -284,9 +283,7 @@ subroutine crystallite_init
crystallite_outputID(o,c) = orientation_ID crystallite_outputID(o,c) = orientation_ID
case ('grainrotation') outputName case ('grainrotation') outputName
crystallite_outputID(o,c) = grainrotation_ID crystallite_outputID(o,c) = grainrotation_ID
case ('eulerangles') outputName case ('defgrad','f') outputName ! ToDo: no alias (f only)
crystallite_outputID(o,c) = eulerangles_ID
case ('defgrad','f') outputName
crystallite_outputID(o,c) = defgrad_ID crystallite_outputID(o,c) = defgrad_ID
case ('fe') outputName case ('fe') outputName
crystallite_outputID(o,c) = fe_ID crystallite_outputID(o,c) = fe_ID
@ -300,13 +297,13 @@ subroutine crystallite_init
crystallite_outputID(o,c) = li_ID crystallite_outputID(o,c) = li_ID
case ('p','firstpiola','1stpiola') outputName case ('p','firstpiola','1stpiola') outputName
crystallite_outputID(o,c) = p_ID crystallite_outputID(o,c) = p_ID
case ('s','tstar','secondpiola','2ndpiola') outputName case ('s','tstar','secondpiola','2ndpiola') outputName ! ToDo: no alias (s only)
crystallite_outputID(o,c) = s_ID crystallite_outputID(o,c) = s_ID
case ('elasmatrix') outputName case ('elasmatrix') outputName
crystallite_outputID(o,c) = elasmatrix_ID crystallite_outputID(o,c) = elasmatrix_ID
case ('neighboringip') outputName case ('neighboringip') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh
crystallite_outputID(o,c) = neighboringip_ID crystallite_outputID(o,c) = neighboringip_ID
case ('neighboringelement') outputName case ('neighboringelement') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh
crystallite_outputID(o,c) = neighboringelement_ID crystallite_outputID(o,c) = neighboringelement_ID
case default outputName case default outputName
call IO_error(105_pInt,ext_msg=trim(str(o))//' (Crystallite)') call IO_error(105_pInt,ext_msg=trim(str(o))//' (Crystallite)')
@ -322,8 +319,6 @@ subroutine crystallite_init
mySize = 1_pInt mySize = 1_pInt
case(orientation_ID,grainrotation_ID) case(orientation_ID,grainrotation_ID)
mySize = 4_pInt mySize = 4_pInt
case(eulerangles_ID)
mySize = 3_pInt
case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID) case(defgrad_ID,fe_ID,fp_ID,fi_ID,lp_ID,li_ID,p_ID,s_ID)
mySize = 9_pInt mySize = 9_pInt
case(elasmatrix_ID) case(elasmatrix_ID)
@ -394,8 +389,7 @@ subroutine crystallite_init
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
call constitutive_microstructure(crystallite_orientation, & call constitutive_microstructure(crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fp(1:3,1:3,c,i,e), & crystallite_Fp(1:3,1:3,c,i,e), &
c,i,e) ! update dependent state variables to be consistent with basic states c,i,e) ! update dependent state variables to be consistent with basic states
enddo enddo
@ -426,7 +420,7 @@ end subroutine crystallite_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P) !> @brief calculate stress (P)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_stress(a) function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
use prec, only: & use prec, only: &
tol_math_check, & tol_math_check, &
dNeq0 dNeq0
@ -462,14 +456,11 @@ function crystallite_stress(a)
sourceState, & sourceState, &
phase_Nsources, & phase_Nsources, &
phaseAt, phasememberAt phaseAt, phasememberAt
use constitutive, only: &
constitutive_SandItsTangents, &
constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents
implicit none implicit none
logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress
real(pReal), intent(in), optional :: a !ToDo: for some reason this prevents an internal compiler error in GNU. Very strange real(pReal), intent(in), optional :: &
dummyArgumentToPreventInternalCompilerErrorWithGCC
real(pReal) :: & real(pReal) :: &
formerSubStep formerSubStep
integer(pInt) :: & integer(pInt) :: &
@ -764,7 +755,7 @@ subroutine crystallite_stressTangent()
crystallite_Fe(1:3,1:3,c,i,e), & crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, &
crystallite_Tstar_v(1:6,c,i,e), & math_6toSym33(crystallite_Tstar_v(1:6,c,i,e)), &
crystallite_Fi(1:3,1:3,c,i,e), & crystallite_Fi(1:3,1:3,c,i,e), &
c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration
@ -793,7 +784,7 @@ subroutine crystallite_stressTangent()
endif endif
call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, & call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
crystallite_Tstar_v(1:6,c,i,e), & math_6toSym33(crystallite_Tstar_v(1:6,c,i,e)), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
@ -894,9 +885,7 @@ subroutine crystallite_orientations
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e)) do c = 1_pInt,homogenization_Ngrains(mesh_element(3,e))
crystallite_orientation(1:4,c,i,e) = math_RtoQ(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e)))) call crystallite_orientation(c,i,e)%fromRotationMatrix(transpose(math_rotationalPart33(crystallite_Fe(1:3,1:3,c,i,e))))
crystallite_rotation(1:4,c,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,c,i,e), &! active rotation from initial
crystallite_orientation(1:4,c,i,e)) ! to current orientation (with no symmetry)
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -969,6 +958,8 @@ function crystallite_postResults(ipc, ip, el)
use constitutive, only: & use constitutive, only: &
constitutive_homogenizedC, & constitutive_homogenizedC, &
constitutive_postResults constitutive_postResults
use rotations, only: &
rotation
implicit none implicit none
integer(pInt), intent(in):: & integer(pInt), intent(in):: &
@ -988,6 +979,7 @@ function crystallite_postResults(ipc, ip, el)
crystID, & crystID, &
mySize, & mySize, &
n n
type(rotation) :: rot
crystID = microstructure_crystallite(mesh_element(4,el)) crystID = microstructure_crystallite(mesh_element(4,el))
@ -1011,15 +1003,12 @@ function crystallite_postResults(ipc, ip, el)
/ real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute) / real(homogenization_Ngrains(mesh_element(3,el)),pReal) ! grain volume (not fraction but absolute)
case (orientation_ID) case (orientation_ID)
mySize = 4_pInt mySize = 4_pInt
crystallite_postResults(c+1:c+mySize) = crystallite_orientation(1:4,ipc,ip,el) ! grain orientation as quaternion crystallite_postResults(c+1:c+mySize) = crystallite_orientation(ipc,ip,el)%asQuaternion()
case (eulerangles_ID)
mySize = 3_pInt
crystallite_postResults(c+1:c+mySize) = inDeg &
* math_qToEuler(crystallite_orientation(1:4,ipc,ip,el)) ! grain orientation as Euler angles in degree
case (grainrotation_ID) case (grainrotation_ID)
rot = crystallite_orientation0(ipc,ip,el)%misorientation(crystallite_orientation(ipc,ip,el))
mySize = 4_pInt mySize = 4_pInt
crystallite_postResults(c+1:c+mySize) = & crystallite_postResults(c+1:c+mySize) = rot%asAxisAnglePair()
math_qToEulerAxisAngle(crystallite_rotation(1:4,ipc,ip,el)) ! grain rotation away from initial orientation as axis-angle in sample reference coordinates
crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree crystallite_postResults(c+4) = inDeg * crystallite_postResults(c+4) ! angle in degree
! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33 ! remark: tensor output is of the form 11,12,13, 21,22,23, 31,32,33
@ -1078,8 +1067,8 @@ function crystallite_postResults(ipc, ip, el)
c = c + 1_pInt c = c + 1_pInt
if (size(crystallite_postResults)-c > 0_pInt) & if (size(crystallite_postResults)-c > 0_pInt) &
crystallite_postResults(c+1:size(crystallite_postResults)) = & crystallite_postResults(c+1:size(crystallite_postResults)) = &
constitutive_postResults(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fi(1:3,1:3,ipc,ip,el), & constitutive_postResults(math_6toSym33(crystallite_Tstar_v(1:6,ipc,ip,el)), crystallite_Fi(1:3,1:3,ipc,ip,el), &
crystallite_Fe, ipc, ip, el) ipc, ip, el)
end function crystallite_postResults end function crystallite_postResults
@ -1289,7 +1278,7 @@ logical function integrateStress(&
!* calculate plastic velocity gradient and its tangent from constitutive law !* calculate plastic velocity gradient and its tangent from constitutive law
call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, & call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
math_sym33to6(S), Fi_new, ipc, ip, el) S, Fi_new, ipc, ip, el)
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
@ -1390,7 +1379,7 @@ logical function integrateStress(&
!* calculate intermediate velocity gradient and its tangent from constitutive law !* calculate intermediate velocity gradient and its tangent from constitutive law
call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, & call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
math_sym33to6(S), Fi_new, ipc, ip, el) S, Fi_new, ipc, ip, el)
#ifdef DEBUG #ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
@ -2206,8 +2195,7 @@ subroutine update_dependentState()
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,homogenization_Ngrains(mesh_element(3,e)) do g = 1,homogenization_Ngrains(mesh_element(3,e))
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) & if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_dependentState(crystallite_orientation, & call constitutive_dependentState(crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) g, i, e)
enddo; enddo; enddo enddo; enddo; enddo
@ -2271,6 +2259,8 @@ end subroutine update_state
subroutine update_dotState(timeFraction) subroutine update_dotState(timeFraction)
use, intrinsic :: & use, intrinsic :: &
IEEE_arithmetic IEEE_arithmetic
use math, only: &
math_6toSym33 !ToDo: Temporarly needed until T_star_v is called S and stored as matrix
use material, only: & use material, only: &
plasticState, & plasticState, &
sourceState, & sourceState, &
@ -2303,11 +2293,11 @@ subroutine update_dotState(timeFraction)
do g = 1,homogenization_Ngrains(mesh_element(3,e)) do g = 1,homogenization_Ngrains(mesh_element(3,e))
!$OMP FLUSH(nonlocalStop) !$OMP FLUSH(nonlocalStop)
if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then if ((crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) .and. .not. nonlocalStop) then
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), & call constitutive_collectDotState(math_6toSym33(crystallite_Tstar_v(1:6,g,i,e)), &
crystallite_Fe, & crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), & crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, & crystallite_Fp, &
crystallite_subdt(g,i,e)*timeFraction, crystallite_subFrac, g,i,e) crystallite_subdt(g,i,e)*timeFraction, g,i,e)
p = phaseAt(g,i,e); c = phasememberAt(g,i,e) p = phaseAt(g,i,e); c = phasememberAt(g,i,e)
NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c))) NaN = any(IEEE_is_NaN(plasticState(p)%dotState(:,c)))
do s = 1_pInt, phase_Nsources(p) do s = 1_pInt, phase_Nsources(p)

169
src/kinematics_thermal_expansion.f90.orig
View File

@ -1,169 +0,0 @@
!--------------------------------------------------------------------------------------------------
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine incorporating kinematics resulting from thermal expansion
!> @details to be done
!--------------------------------------------------------------------------------------------------
module kinematics_thermal_expansion
use prec, only: &
pReal, &
pInt
implicit none
private
<<<<<<< HEAD
!type, private :: tParameters
! real(pReal), allocatable, dimension(:) :: &
!end type tParameters
=======
integer(pInt), dimension(:), allocatable, public, protected :: &
kinematics_thermal_expansion_sizePostResults, & !< cumulative size of post results
kinematics_thermal_expansion_offset, & !< which kinematics is my current damage mechanism?
kinematics_thermal_expansion_instance !< instance of damage kinematics mechanism
integer(pInt), dimension(:,:), allocatable, target, public :: &
kinematics_thermal_expansion_sizePostResult !< size of each post result output
character(len=64), dimension(:,:), allocatable, target, public :: &
kinematics_thermal_expansion_output !< name of each post result output
integer(pInt), dimension(:), allocatable, target, public :: &
kinematics_thermal_expansion_Noutput !< number of outputs per instance of this damage
enum, bind(c) ! ToDo kinematics need state machinery to deal with sizePostResult
enumerator :: undefined_ID, & ! possible remedy is to decouple having state vars from having output
thermalexpansionrate_ID ! which means to separate user-defined types tState + tOutput...
end enum
>>>>>>> development
public :: &
kinematics_thermal_expansion_init, &
kinematics_thermal_expansion_initialStrain, &
kinematics_thermal_expansion_LiAndItsTangent
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_init()
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
use, intrinsic :: iso_fortran_env, only: &
compiler_version, &
compiler_options
#endif
use debug, only: &
debug_level,&
debug_constitutive,&
debug_levelBasic
use IO, only: &
IO_timeStamp
use material, only: &
phase_kinematics, &
KINEMATICS_thermal_expansion_label, &
KINEMATICS_thermal_expansion_ID
use config, only: &
config_phase
implicit none
integer(pInt) :: &
Ninstance, &
p
write(6,'(/,a)') ' <<<+- kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init -+>>>'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
Ninstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt)
if (Ninstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',Ninstance
do p = 1_pInt, size(phase_kinematics)
if (all(phase_kinematics(:,p) /= KINEMATICS_thermal_expansion_ID)) cycle
enddo
end subroutine kinematics_thermal_expansion_init
!--------------------------------------------------------------------------------------------------
!> @brief report initial thermal strain based on current temperature deviation from reference
!--------------------------------------------------------------------------------------------------
pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset)
use material, only: &
temperature
use lattice, only: &
lattice_thermalExpansion33, &
lattice_referenceTemperature
implicit none
integer(pInt), intent(in) :: &
phase, &
homog, offset
real(pReal), dimension(3,3) :: &
kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though)
kinematics_thermal_expansion_initialStrain = &
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**1 / 1. * &
lattice_thermalExpansion33(1:3,1:3,1,phase) + & ! constant coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**2 / 2. * &
lattice_thermalExpansion33(1:3,1:3,2,phase) + & ! linear coefficient
(temperature(homog)%p(offset) - lattice_referenceTemperature(phase))**3 / 3. * &
lattice_thermalExpansion33(1:3,1:3,3,phase) ! quadratic coefficient
end function kinematics_thermal_expansion_initialStrain
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)
use material, only: &
material_phase, &
material_homog, &
temperature, &
temperatureRate, &
thermalMapping
use lattice, only: &
lattice_thermalExpansion33, &
lattice_referenceTemperature
implicit none
integer(pInt), intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
integer(pInt) :: &
phase, &
homog, offset
real(pReal) :: &
T, TRef, TDot
phase = material_phase(ipc,ip,el)
homog = material_homog(ip,el)
offset = thermalMapping(homog)%p(ip,el)
T = temperature(homog)%p(offset)
TDot = temperatureRate(homog)%p(offset)
TRef = lattice_referenceTemperature(phase)
Li = TDot * ( &
lattice_thermalExpansion33(1:3,1:3,1,phase)*(T - TRef)**0 & ! constant coefficient
+ lattice_thermalExpansion33(1:3,1:3,2,phase)*(T - TRef)**1 & ! linear coefficient
+ lattice_thermalExpansion33(1:3,1:3,3,phase)*(T - TRef)**2 & ! quadratic coefficient
) / &
(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. &
)
dLi_dTstar = 0.0_pReal
end subroutine kinematics_thermal_expansion_LiAndItsTangent
end module kinematics_thermal_expansion

716
src/lattice.f90
View File

File diff suppressed because it is too large Load Diff

28
src/math.f90
View File

@ -70,6 +70,10 @@ module math
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Provide deprecated names for compatibility ! Provide deprecated names for compatibility
interface math_cross
module procedure math_crossproduct
end interface math_cross
! ToDo MD: Our naming scheme was a little bit odd: We use essentially the re-ordering according to Nye ! ToDo MD: Our naming scheme was a little bit odd: We use essentially the re-ordering according to Nye
! (convenient because Abaqus and Marc want to have 12 on position 4) ! (convenient because Abaqus and Marc want to have 12 on position 4)
! but weight the shear components according to Mandel (convenient for matrix multiplications) ! but weight the shear components according to Mandel (convenient for matrix multiplications)
@ -98,23 +102,13 @@ module math
module procedure math_99to3333 module procedure math_99to3333
end interface math_Plain99to3333 end interface math_Plain99to3333
interface math_Mandel3333to66
module procedure math_sym3333to66
end interface math_Mandel3333to66
interface math_Mandel66to3333
module procedure math_66toSym3333
end interface math_Mandel66to3333
public :: & public :: &
math_Plain33to9, & math_Plain33to9, &
math_Plain9to33, & math_Plain9to33, &
math_Mandel33to6, & math_Mandel33to6, &
math_Mandel6to33, & math_Mandel6to33, &
math_Plain3333to99, & math_Plain3333to99, &
math_Plain99to3333, & math_Plain99to3333
math_Mandel3333to66, &
math_Mandel66to3333
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
public :: & public :: &
@ -129,6 +123,7 @@ module math
math_identity4th, & math_identity4th, &
math_civita, & math_civita, &
math_delta, & math_delta, &
math_cross, &
math_crossproduct, & math_crossproduct, &
math_tensorproduct33, & math_tensorproduct33, &
math_mul3x3, & math_mul3x3, &
@ -1887,7 +1882,6 @@ function math_sampleGaussOri(center,FWHM)
math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center))) math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center)))
endif endif
end function math_sampleGaussOri end function math_sampleGaussOri
@ -1960,11 +1954,11 @@ real(pReal) function math_sampleGaussVar(meanvalue, stddev, width)
tol_math_check tol_math_check
implicit none implicit none
real(pReal), intent(in) :: meanvalue, & ! meanvalue of gauss distribution real(pReal), intent(in) :: meanvalue, & ! meanvalue of gauss distribution
stddev ! standard deviation of gauss distribution stddev ! standard deviation of gauss distribution
real(pReal), intent(in), optional :: width ! width of considered values as multiples of standard deviation real(pReal), intent(in), optional :: width ! width of considered values as multiples of standard deviation
real(pReal), dimension(2) :: rnd ! random numbers real(pReal), dimension(2) :: rnd ! random numbers
real(pReal) :: scatter, & ! normalized scatter around meanvalue real(pReal) :: scatter, & ! normalized scatter around meanvalue
myWidth myWidth
if (abs(stddev) < tol_math_check) then if (abs(stddev) < tol_math_check) then

2872
src/plastic_nonlocal.f90
View File

File diff suppressed because it is too large Load Diff

4
src/prec.f90
View File

@ -29,8 +29,8 @@ module prec
real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation) real(pReal), parameter, public :: tol_math_check = 1.0e-8_pReal !< tolerance for internal math self-checks (rotation)
integer(pInt), allocatable, dimension(:) :: realloc_lhs_test integer(pInt), allocatable, dimension(:) :: realloc_lhs_test
type, public :: group_float !< variable length datatype used for storage of state type, public :: group_float !< variable length datatype used for storage of state
real(pReal), dimension(:), pointer :: p real(pReal), dimension(:), pointer :: p
end type group_float end type group_float

443
src/quaternions.f90 Normal file
View File

@ -0,0 +1,443 @@
! ###################################################################
! Copyright (c) 2013-2015, Marc De Graef/Carnegie Mellon University
! Modified 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without modification, are
! permitted provided that the following conditions are met:
!
! - Redistributions of source code must retain the above copyright notice, this list
! of conditions and the following disclaimer.
! - Redistributions in binary form must reproduce the above copyright notice, this
! list of conditions and the following disclaimer in the documentation and/or
! other materials provided with the distribution.
! - Neither the names of Marc De Graef, Carnegie Mellon University nor the names
! of its contributors may be used to endorse or promote products derived from
! this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
! ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
! LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE
! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
! ###################################################################
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief general quaternion math, not limited to unit quaternions
!> @details w is the real part, (x, y, z) are the imaginary parts.
!---------------------------------------------------------------------------------------------------
module quaternions
use prec, only: &
pReal
implicit none
public
real(pReal), parameter, public :: P = -1.0_pReal !< parameter for orientation conversion.
type, public :: quaternion
real(pReal) :: w = 0.0_pReal
real(pReal) :: x = 0.0_pReal
real(pReal) :: y = 0.0_pReal
real(pReal) :: z = 0.0_pReal
contains
procedure, private :: add__
procedure, private :: pos__
generic, public :: operator(+) => add__,pos__
procedure, private :: sub__
procedure, private :: neg__
generic, public :: operator(-) => sub__,neg__
procedure, private :: mul_quat__
procedure, private :: mul_scal__
generic, public :: operator(*) => mul_quat__, mul_scal__
procedure, private :: div_quat__
procedure, private :: div_scal__
generic, public :: operator(/) => div_quat__, div_scal__
procedure, private :: eq__
generic, public :: operator(==) => eq__
procedure, private :: neq__
generic, public :: operator(/=) => neq__
procedure, private :: pow_quat__
procedure, private :: pow_scal__
generic, public :: operator(**) => pow_quat__, pow_scal__
procedure, private :: abs__
procedure, private :: dot_product__
procedure, private :: conjg__
procedure, private :: exp__
procedure, private :: log__
procedure, public :: homomorphed => quat_homomorphed
end type
interface assignment (=)
module procedure assign_quat__
module procedure assign_vec__
end interface assignment (=)
interface quaternion
module procedure init__
end interface quaternion
interface abs
procedure abs__
end interface abs
interface dot_product
procedure dot_product__
end interface dot_product
interface conjg
module procedure conjg__
end interface conjg
interface exp
module procedure exp__
end interface exp
interface log
module procedure log__
end interface log
contains
!---------------------------------------------------------------------------------------------------
!> constructor for a quaternion from a 4-vector
!---------------------------------------------------------------------------------------------------
type(quaternion) pure function init__(array)
implicit none
real(pReal), intent(in), dimension(4) :: array
init__%w=array(1)
init__%x=array(2)
init__%y=array(3)
init__%z=array(4)
end function init__
!---------------------------------------------------------------------------------------------------
!> assing a quaternion
!---------------------------------------------------------------------------------------------------
elemental subroutine assign_quat__(self,other)
implicit none
type(quaternion), intent(out) :: self
type(quaternion), intent(in) :: other
self%w = other%w
self%x = other%x
self%y = other%y
self%z = other%z
end subroutine assign_quat__
!---------------------------------------------------------------------------------------------------
!> assing a 4-vector
!---------------------------------------------------------------------------------------------------
pure subroutine assign_vec__(self,other)
implicit none
type(quaternion), intent(out) :: self
real(pReal), intent(in), dimension(4) :: other
self%w = other(1)
self%x = other(2)
self%y = other(3)
self%z = other(4)
end subroutine assign_vec__
!---------------------------------------------------------------------------------------------------
!> addition of two quaternions
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function add__(self,other)
implicit none
class(quaternion), intent(in) :: self,other
add__%w = self%w + other%w
add__%x = self%x + other%x
add__%y = self%y + other%y
add__%z = self%z + other%z
end function add__
!---------------------------------------------------------------------------------------------------
!> unary positive operator
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pos__(self)
implicit none
class(quaternion), intent(in) :: self
pos__%w = self%w
pos__%x = self%x
pos__%y = self%y
pos__%z = self%z
end function pos__
!---------------------------------------------------------------------------------------------------
!> subtraction of two quaternions
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function sub__(self,other)
implicit none
class(quaternion), intent(in) :: self,other
sub__%w = self%w - other%w
sub__%x = self%x - other%x
sub__%y = self%y - other%y
sub__%z = self%z - other%z
end function sub__
!---------------------------------------------------------------------------------------------------
!> unary positive operator
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function neg__(self)
implicit none
class(quaternion), intent(in) :: self
neg__%w = -self%w
neg__%x = -self%x
neg__%y = -self%y
neg__%z = -self%z
end function neg__
!---------------------------------------------------------------------------------------------------
!> multiplication of two quaternions
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function mul_quat__(self,other)
implicit none
class(quaternion), intent(in) :: self, other
mul_quat__%w = self%w*other%w - self%x*other%x - self%y*other%y - self%z*other%z
mul_quat__%x = self%w*other%x + self%x*other%w + P * (self%y*other%z - self%z*other%y)
mul_quat__%y = self%w*other%y + self%y*other%w + P * (self%z*other%x - self%x*other%z)
mul_quat__%z = self%w*other%z + self%z*other%w + P * (self%x*other%y - self%y*other%x)
end function mul_quat__
!---------------------------------------------------------------------------------------------------
!> multiplication of quaternions with scalar
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function mul_scal__(self,scal)
implicit none
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal
mul_scal__%w = self%w*scal
mul_scal__%x = self%x*scal
mul_scal__%y = self%y*scal
mul_scal__%z = self%z*scal
end function mul_scal__
!---------------------------------------------------------------------------------------------------
!> division of two quaternions
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function div_quat__(self,other)
implicit none
class(quaternion), intent(in) :: self, other
div_quat__ = self * (conjg(other)/(abs(other)**2.0_pReal))
end function div_quat__
!---------------------------------------------------------------------------------------------------
!> divisiont of quaternions by scalar
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function div_scal__(self,scal)
implicit none
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: scal
div_scal__ = [self%w,self%x,self%y,self%z]/scal
end function div_scal__
!---------------------------------------------------------------------------------------------------
!> equality of two quaternions
!---------------------------------------------------------------------------------------------------
logical elemental function eq__(self,other)
use prec, only: &
dEq
implicit none
class(quaternion), intent(in) :: self,other
eq__ = all(dEq([ self%w, self%x, self%y, self%z], &
[other%w,other%x,other%y,other%z]))
end function eq__
!---------------------------------------------------------------------------------------------------
!> inequality of two quaternions
!---------------------------------------------------------------------------------------------------
logical elemental function neq__(self,other)
implicit none
class(quaternion), intent(in) :: self,other
neq__ = .not. self%eq__(other)
end function neq__
!---------------------------------------------------------------------------------------------------
!> quaternion to the power of a scalar
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pow_scal__(self,expon)
implicit none
class(quaternion), intent(in) :: self
real(pReal), intent(in) :: expon
pow_scal__ = exp(log(self)*expon)
end function pow_scal__
!---------------------------------------------------------------------------------------------------
!> quaternion to the power of a quaternion
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function pow_quat__(self,expon)
implicit none
class(quaternion), intent(in) :: self
type(quaternion), intent(in) :: expon
pow_quat__ = exp(log(self)*expon)
end function pow_quat__
!---------------------------------------------------------------------------------------------------
!> exponential of a quaternion
!> ToDo: Lacks any check for invalid operations
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function exp__(self)
implicit none
class(quaternion), intent(in) :: self
real(pReal) :: absImag
absImag = norm2([self%x, self%y, self%z])
exp__ = exp(self%w) * [ cos(absImag), &
self%x/absImag * sin(absImag), &
self%y/absImag * sin(absImag), &
self%z/absImag * sin(absImag)]
end function exp__
!---------------------------------------------------------------------------------------------------
!> logarithm of a quaternion
!> ToDo: Lacks any check for invalid operations
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function log__(self)
implicit none
class(quaternion), intent(in) :: self
real(pReal) :: absImag
absImag = norm2([self%x, self%y, self%z])
log__ = [log(abs(self)), &
self%x/absImag * acos(self%w/abs(self)), &
self%y/absImag * acos(self%w/abs(self)), &
self%z/absImag * acos(self%w/abs(self))]
end function log__
!---------------------------------------------------------------------------------------------------
!> norm of a quaternion
!---------------------------------------------------------------------------------------------------
real(pReal) elemental function abs__(a)
implicit none
class(quaternion), intent(in) :: a
abs__ = norm2([a%w,a%x,a%y,a%z])
end function abs__
!---------------------------------------------------------------------------------------------------
!> dot product of two quaternions
!---------------------------------------------------------------------------------------------------
real(pReal) elemental function dot_product__(a,b)
implicit none
class(quaternion), intent(in) :: a,b
dot_product__ = a%w*b%w + a%x*b%x + a%y*b%y + a%z*b%z
end function dot_product__
!---------------------------------------------------------------------------------------------------
!> conjugate complex of a quaternion
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function conjg__(a)
implicit none
class(quaternion), intent(in) :: a
conjg__ = quaternion([a%w, -a%x, -a%y, -a%z])
end function conjg__
!---------------------------------------------------------------------------------------------------
!> homomorphed quaternion of a quaternion
!---------------------------------------------------------------------------------------------------
type(quaternion) elemental function quat_homomorphed(a)
implicit none
class(quaternion), intent(in) :: a
quat_homomorphed = quaternion(-[a%w,a%x,a%y,a%z])
end function quat_homomorphed
end module quaternions

1198
src/rotations.f90 Normal file
View File

File diff suppressed because it is too large Load Diff