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

Merge branch 'cleaning' into 'development' 

Cleaning

See merge request damask/DAMASK!237
This commit is contained in:
Philip Eisenlohr 2020-09-28 00:28:11 +02:00
parent c982658b11 8ec1014134
commit c545ad869d
12 changed files with 45 additions and 1047 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

57
.gitlab-ci.yml
View File

@ -239,6 +239,16 @@ Compile_Intel_Prepare:
- release
###################################################################################################
Pytest_grid:
stage: grid
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cd pytest
- pytest
except:
- master
- release
Thermal:
stage: grid
script: Thermal/test.py
@ -253,27 +263,6 @@ grid_parsingArguments:
- master
- release
StateIntegration_compareVariants:
stage: grid
script: StateIntegration_compareVariants/test.py
except:
- master
- release
nonlocal_densityConservation:
stage: grid
script: nonlocal_densityConservation/test.py
except:
- master
- release
RGC_DetectChanges:
stage: grid
script: RGC_DetectChanges/test.py
except:
- master
- release
Nonlocal_Damage_DetectChanges:
stage: grid
script: Nonlocal_Damage_DetectChanges/test.py
@ -302,15 +291,6 @@ grid_all_loadCaseRotation:
- master
- release
grid_mech_MPI:
stage: grid
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- grid_mech_MPI/test.py
except:
- master
- release
grid_all_restartMPI:
stage: grid
script:
@ -327,13 +307,6 @@ Plasticity_DetectChanges:
- master
- release
Homogenization:
stage: grid
script: Homogenization/test.py
except:
- master
- release
Phenopowerlaw_singleSlip:
stage: grid
script: Phenopowerlaw_singleSlip/test.py
@ -341,16 +314,6 @@ Phenopowerlaw_singleSlip:
- master
- release
Pytest_grid:
stage: grid
script:
- module load $IntelCompiler $MPICH_Intel $PETSc_MPICH_Intel
- cd pytest
- pytest
except:
- master
- release
###################################################################################################
Marc_compileIfort:
stage: compileMarc

8
Makefile
View File

@ -9,18 +9,10 @@ all: grid mesh processing
.PHONY: grid
grid: build/grid
@(cd build/grid;make -j${DAMASK_NUM_THREADS} all install;)
@rm -f ${DAMASK_ROOT}/bin/DAMASK_spectral > /dev/null || true
@ln -s ${DAMASK_ROOT}/bin/DAMASK_grid ${DAMASK_ROOT}/bin/DAMASK_spectral || true
.PHONY: spectral
spectral: grid
.PHONY: mesh
mesh: build/mesh
@(cd build/mesh; make -j${DAMASK_NUM_THREADS} all install;)
@rm -f ${DAMASK_ROOT}/bin/DAMASK_FEM > /dev/null || true
@ln -s ${DAMASK_ROOT}/bin/DAMASK_mesh ${DAMASK_ROOT}/bin/DAMASK_FEM || true
.PHONY: FEM
FEM: mesh
.PHONY: build/grid
build/grid:

2
PRIVATE

@ -1 +1 @@
Subproject commit 68837540cab7435d8e2a06ae4c74e069e9386f35
Subproject commit 25ce39dd0f5bf49dc5c2bec20767d93d2b76d353

0
processing/post/vtk2ang.py → processing/legacy/vtk2ang.py
View File

161
processing/post/filterTable.py
View File

@ -1,161 +0,0 @@
#!/usr/bin/env python3
import os
import sys
from optparse import OptionParser
import re
import fnmatch
import math # noqa
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
def sortingList(labels,whitelistitems):
indices = []
names = []
for label in labels:
if re.match(r'^\d+_',label):
indices.append(int(label.split('_',1)[0]))
names.append(label.split('_',1)[1])
else:
indices.append(0)
names.append(label)
return [indices,names,whitelistitems]
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [ASCIItable(s)]', description = """
Filter rows according to condition and columns by either white or black listing.
Examples:
Every odd row if x coordinate is positive -- " #ip.x# >= 0.0 and #_row_#%2 == 1 ).
All rows where label 'foo' equals 'bar' -- " #s#foo# == 'bar' "
""", version = scriptID)
parser.add_option('-w','--white',
dest = 'whitelist',
action = 'extend', metavar = '<string LIST>',
help = 'whitelist of column labels (a,b,c,...)')
parser.add_option('-b','--black',
dest = 'blacklist',
action = 'extend', metavar='<string LIST>',
help = 'blacklist of column labels (a,b,c,...)')
parser.add_option('-c','--condition',
dest = 'condition', metavar='string',
help = 'condition to filter rows')
parser.set_defaults(condition = None,
)
(options,filenames) = parser.parse_args()
# --- loop over input files -------------------------------------------------------------------------
if filenames == []: filenames = [None]
for name in filenames:
try:
table = damask.ASCIItable(name = name)
except IOError:
continue
damask.util.report(scriptName,name)
# ------------------------------------------ assemble info ---------------------------------------
table.head_read()
# ------------------------------------------ process data ---------------------------------------
specials = { \
'_row_': 0,
}
labels = []
positions = []
for position,label in enumerate(table.labels(raw = True)):
if (options.whitelist is None or any([ position in table.label_indexrange(needle) \
or fnmatch.fnmatch(label,needle) for needle in options.whitelist])) \
and (options.blacklist is None or not any([ position in table.label_indexrange(needle) \
or fnmatch.fnmatch(label,needle) for needle in options.blacklist])): # a label to keep?
labels.append(label) # remember name...
positions.append(position) # ...and position
if len(labels) > 0 and options.whitelist is not None and options.blacklist is None: # check whether reordering is possible
whitelistitem = np.zeros(len(labels),dtype=int)
for i,label in enumerate(labels): # check each selected label
match = [ positions[i] in table.label_indexrange(needle) \
or fnmatch.fnmatch(label,needle) for needle in options.whitelist] # which whitelist items do match it
whitelistitem[i] = match.index(True) if np.sum(match) == 1 else -1 # unique match to a whitelist item --> store which
order = range(len(labels)) if np.any(whitelistitem < 0) \
else np.lexsort(sortingList(labels,whitelistitem)) # reorder if unique, i.e. no "-1" in whitelistitem
else:
order = range(len(labels)) # maintain original order of labels
# --------------------------------------- evaluate condition ---------------------------------------
if options.condition is not None:
condition = options.condition # copy per file, since might be altered inline
breaker = False
for position,(all,marker,column) in enumerate(set(re.findall(r'#(([s]#)?(.+?))#',condition))): # find three groups
idx = table.label_index(column)
dim = table.label_dimension(column)
if idx < 0 and column not in specials:
damask.util.croak('column "{}" not found.'.format(column))
breaker = True
else:
if column in specials:
replacement = 'specials["{}"]'.format(column)
elif dim == 1: # scalar input
replacement = '{}(table.data[{}])'.format({ '':'float',
's#':'str'}[marker],idx) # take float or string value of data column
elif dim > 1: # multidimensional input (vector, tensor, etc.)
replacement = 'np.array(table.data[{}:{}],dtype=float)'.format(idx,idx+dim) # use (flat) array representation
condition = condition.replace('#'+all+'#',replacement)
if breaker: continue # found mistake in condition evaluation --> next file
# ------------------------------------------ assemble header ---------------------------------------
table.info_append(scriptID + '\t' + ' '.join(sys.argv[1:]))
table.labels_clear()
table.labels_append(np.array(labels)[order]) # update with new label set
table.head_write()
# ------------------------------------------ process and output data ------------------------------------------
positions = np.array(positions)[order]
atOnce = options.condition is None
if atOnce: # read full array and filter columns
try:
table.data_readArray(positions+1) # read desired columns (indexed 1,...)
table.data_writeArray() # directly write out
except Exception:
table.data_rewind()
atOnce = False # data contains items that prevent array chunking
if not atOnce: # read data line by line
outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table
specials['_row_'] += 1 # count row
if options.condition is None or eval(condition): # valid row ?
table.data = [table.data[position] for position in positions] # retain filtered columns
outputAlive = table.data_write() # output processed line
# ------------------------------------------ finalize output -----------------------------------------
table.close() # close input ASCII table (works for stdin)

226
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -1,226 +0,0 @@
#!/usr/bin/env python3
import os
import sys
import multiprocessing
from io import StringIO
from functools import partial
from optparse import OptionParser,OptionGroup
import numpy as np
from scipy import spatial
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
def findClosestSeed(seeds, weights, point):
return np.argmin(np.sum((np.broadcast_to(point,(len(seeds),3))-seeds)**2,axis=1) - weights)
def Laguerre_tessellation(grid, size, seeds, weights, origin = np.zeros(3), periodic = True, cpus = 2):
if periodic:
weights_p = np.tile(weights.squeeze(),27) # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
seeds_p = np.vstack((seeds -np.array([size[0],0.,0.]),seeds, seeds +np.array([size[0],0.,0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,size[1],0.]),seeds_p,seeds_p+np.array([0.,size[1],0.])))
seeds_p = np.vstack((seeds_p-np.array([0.,0.,size[2]]),seeds_p,seeds_p+np.array([0.,0.,size[2]])))
coords = damask.grid_filters.cell_coord0(grid*3,size*3,-origin-size).reshape(-1,3)
else:
weights_p = weights.squeeze()
seeds_p = seeds
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3)
if cpus > 1:
pool = multiprocessing.Pool(processes = cpus)
result = pool.map_async(partial(findClosestSeed,seeds_p,weights_p), [coord for coord in coords])
pool.close()
pool.join()
closest_seed = np.array(result.get()).reshape(-1,3)
else:
closest_seed= np.array([findClosestSeed(seeds_p,weights_p,coord) for coord in coords])
if periodic:
closest_seed = closest_seed.reshape(grid*3)
return closest_seed[grid[0]:grid[0]*2,grid[1]:grid[1]*2,grid[2]:grid[2]*2]%seeds.shape[0]
else:
return closest_seed
def Voronoi_tessellation(grid, size, seeds, origin = np.zeros(3), periodic = True):
coords = damask.grid_filters.cell_coord0(grid,size,-origin).reshape(-1,3)
KDTree = spatial.cKDTree(seeds,boxsize=size) if periodic else spatial.cKDTree(seeds)
devNull,closest_seed = KDTree.query(coords)
return closest_seed
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog option(s) [seedfile(s)]', description = """
Generate geometry description and material configuration by tessellation of given seeds file.
""", version = scriptID)
group = OptionGroup(parser, "Tessellation","")
group.add_option('-l',
'--laguerre',
dest = 'laguerre',
action = 'store_true',
help = 'use Laguerre (weighted Voronoi) tessellation')
group.add_option('--cpus',
dest = 'cpus',
type = 'int', metavar = 'int',
help = 'number of parallel processes to use for Laguerre tessellation [%default]')
group.add_option('--nonperiodic',
dest = 'periodic',
action = 'store_false',
help = 'nonperiodic tessellation')
parser.add_option_group(group)
group = OptionGroup(parser, "Geometry","")
group.add_option('-g',
'--grid',
dest = 'grid',
type = 'int', nargs = 3, metavar = ' '.join(['int']*3),
help = 'a,b,c grid of hexahedral box')
group.add_option('-s',
'--size',
dest = 'size',
type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'x,y,z size of hexahedral box [1.0 1.0 1.0]')
group.add_option('-o',
'--origin',
dest = 'origin',
type = 'float', nargs = 3, metavar=' '.join(['float']*3),
help = 'origin of grid [0.0 0.0 0.0]')
parser.add_option_group(group)
group = OptionGroup(parser, "Seeds","")
group.add_option('-p',
'--pos', '--seedposition',
dest = 'pos',
type = 'string', metavar = 'string',
help = 'label of coordinates [%default]')
group.add_option('-w',
'--weight',
dest = 'weight',
type = 'string', metavar = 'string',
help = 'label of weights [%default]')
group.add_option('-m',
'--microstructure',
dest = 'microstructure',
type = 'string', metavar = 'string',
help = 'label of microstructures [%default]')
group.add_option('-e',
'--eulers',
dest = 'eulers',
type = 'string', metavar = 'string',
help = 'label of Euler angles [%default]')
group.add_option('--axes',
dest = 'axes',
type = 'string', nargs = 3, metavar = ' '.join(['string']*3),
help = 'orientation coordinate frame in terms of position coordinate frame')
parser.add_option_group(group)
group = OptionGroup(parser, "Configuration","")
group.add_option('--without-config',
dest = 'config',
action = 'store_false',
help = 'omit material configuration header')
group.add_option('--phase',
dest = 'phase',
type = 'int', metavar = 'int',
help = 'phase index to be used [%default]')
parser.add_option_group(group)
parser.set_defaults(pos = 'pos',
weight = 'weight',
microstructure = 'microstructure',
eulers = 'euler',
phase = 1,
cpus = 2,
laguerre = False,
periodic = True,
config = True,
)
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
for name in filenames:
damask.util.report(scriptName,name)
table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
size = np.ones(3)
origin = np.zeros(3)
for line in table.comments:
items = line.lower().strip().split()
key = items[0] if items else ''
if key == 'grid':
grid = np.array([ int(dict(zip(items[1::2],items[2::2]))[i]) for i in ['a','b','c']])
elif key == 'size':
size = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
elif key == 'origin':
origin = np.array([float(dict(zip(items[1::2],items[2::2]))[i]) for i in ['x','y','z']])
if options.grid: grid = np.array(options.grid)
if options.size: size = np.array(options.size)
if options.origin: origin = np.array(options.origin)
seeds = table.get(options.pos)
grains = table.get(options.microstructure) if options.microstructure in table.labels else np.arange(len(seeds))+1
grainIDs = np.unique(grains).astype('i')
if options.eulers in table.labels:
eulers = table.get(options.eulers)
if options.laguerre:
indices = grains[Laguerre_tessellation(grid,size,seeds,table.get(options.weight),origin,
options.periodic,options.cpus)]
else:
indices = grains[Voronoi_tessellation (grid,size,seeds,origin,options.periodic)]
config_header = []
if options.config:
if options.eulers in table.labels:
config_header += ['<texture>']
for ID in grainIDs:
eulerID = np.nonzero(grains == ID)[0][0] # find first occurrence of this grain id
config_header += ['[Grain{}]'.format(ID),
'(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*eulers[eulerID])
]
if options.axes: config_header += ['axes\t{} {} {}'.format(*options.axes)]
config_header += ['<microstructure>']
for ID in grainIDs:
config_header += ['[Grain{}]'.format(ID),
'(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(options.phase,ID)
]
config_header += ['<!skip>']
header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\
+ config_header
geom = damask.Geom(indices.reshape(grid),size,origin,
comments=header)
damask.util.croak(geom)
geom.save_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',compress=False)

349
processing/pre/hybridIA_linODFsampling.py
View File

@ -1,349 +0,0 @@
#!/usr/bin/env python3
import os
import sys
import math
import random
from io import StringIO
from optparse import OptionParser
import numpy as np
import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version])
# --- helper functions ---
def integerFactorization(i):
j = int(math.floor(math.sqrt(float(i))))
while j>1 and int(i)%j != 0:
j -= 1
return j
def binAsBins(bin,intervals):
"""Explode compound bin into 3D bins list."""
bins = [0]*3
bins[0] = (bin//(intervals[1] * intervals[2])) % intervals[0]
bins[1] = (bin//intervals[2]) % intervals[1]
bins[2] = bin % intervals[2]
return bins
def binsAsBin(bins,intervals):
"""Implode 3D bins into compound bin."""
return (bins[0]*intervals[1] + bins[1])*intervals[2] + bins[2]
def EulersAsBins(Eulers,intervals,deltas,center):
"""Return list of Eulers translated into 3D bins list."""
return [int((euler+(0.5-center)*delta)//delta)%interval \
for euler,delta,interval in zip(Eulers,deltas,intervals) \
]
def binAsEulers(bin,intervals,deltas,center):
"""Compound bin number translated into list of Eulers."""
Eulers = [0.0]*3
Eulers[2] = (bin%intervals[2] + center)*deltas[2]
Eulers[1] = (bin//intervals[2]%intervals[1] + center)*deltas[1]
Eulers[0] = (bin//(intervals[2]*intervals[1]) + center)*deltas[0]
return Eulers
def directInvRepetitions(probability,scale):
"""Calculate number of samples drawn by direct inversion."""
nDirectInv = 0
for bin in range(len(probability)): # loop over bins
nDirectInv += int(round(probability[bin]*scale)) # calc repetition
return nDirectInv
# ---------------------- sampling methods -----------------------------------------------------------------------
# ----- efficient algorithm ---------
def directInversion (ODF,nSamples):
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
nOptSamples = max(ODF['nNonZero'],nSamples) # random subsampling if too little samples requested
nInvSamples = 0
repetition = [None]*ODF['nBins']
scaleLower = 0.0
nInvSamplesLower = 0
scaleUpper = float(nOptSamples)
incFactor = 1.0
nIter = 0
nInvSamplesUpper = directInvRepetitions(ODF['dV_V'],scaleUpper)
while (\
(scaleUpper-scaleLower > scaleUpper*1e-15 or nInvSamplesUpper < nOptSamples) and \
nInvSamplesUpper != nOptSamples \
): # closer match required?
if nInvSamplesUpper < nOptSamples:
scaleLower,scaleUpper = scaleUpper,scaleUpper+incFactor*(scaleUpper-scaleLower)/2.0
incFactor *= 2.0
nInvSamplesLower,nInvSamplesUpper = nInvSamplesUpper,directInvRepetitions(ODF['dV_V'],scaleUpper)
else:
scaleUpper = (scaleLower+scaleUpper)/2.0
incFactor = 1.0
nInvSamplesUpper = directInvRepetitions(ODF['dV_V'],scaleUpper)
nIter += 1
damask.util.croak('%i:(%12.11f,%12.11f) %i <= %i <= %i'%(nIter,scaleLower,scaleUpper,
nInvSamplesLower,nOptSamples,nInvSamplesUpper))
nInvSamples = nInvSamplesUpper
scale = scaleUpper
damask.util.croak('created set of %i samples (%12.11f) with scaling %12.11f delivering %i'%(nInvSamples,
float(nInvSamples)/nOptSamples-1.0,
scale,nSamples))
repetition = [None]*ODF['nBins'] # preallocate and clear
for bin in range(ODF['nBins']): # loop over bins
repetition[bin] = int(round(ODF['dV_V'][bin]*scale)) # calc repetition
# build set
set = [None]*nInvSamples
i = 0
for bin in range(ODF['nBins']):
set[i:i+repetition[bin]] = [bin]*repetition[bin] # fill set with bin, i.e. orientation
i += repetition[bin] # advance set counter
orientations = np.zeros((nSamples,3),'f')
reconstructedODF = np.zeros(ODF['nBins'],'f')
unitInc = 1.0/nSamples
for j in range(nSamples):
if (j == nInvSamples-1): ex = j
else: ex = int(round(random.uniform(j+0.5,nInvSamples-0.5)))
bin = set[ex]
Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
orientations[j] = np.degrees(Eulers)
reconstructedODF[bin] += unitInc
set[ex] = set[j] # exchange orientations
return orientations, reconstructedODF
# ----- trial and error algorithms ---------
def MonteCarloEulers (ODF,nSamples):
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
countMC = 0
maxdV_V = max(ODF['dV_V'])
orientations = np.zeros((nSamples,3),'f')
reconstructedODF = np.zeros(ODF['nBins'],'f')
unitInc = 1.0/nSamples
for j in range(nSamples):
MC = maxdV_V*2.0
bin = 0
while MC > ODF['dV_V'][bin]:
countMC += 1
MC = maxdV_V*random.random()
Eulers = [limit*random.random() for limit in ODF['limit']]
bins = EulersAsBins(Eulers,ODF['interval'],ODF['delta'],ODF['center'])
bin = binsAsBin(bins,ODF['interval'])
orientations[j] = np.degrees(Eulers)
reconstructedODF[bin] += unitInc
return orientations, reconstructedODF, countMC
def MonteCarloBins (ODF,nSamples):
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
countMC = 0
maxdV_V = max(ODF['dV_V'])
orientations = np.zeros((nSamples,3),'f')
reconstructedODF = np.zeros(ODF['nBins'],'f')
unitInc = 1.0/nSamples
for j in range(nSamples):
MC = maxdV_V*2.0
bin = 0
while MC > ODF['dV_V'][bin]:
countMC += 1
MC = maxdV_V*random.random()
bin = int(ODF['nBins'] * random.random())
Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
orientations[j] = np.degrees(Eulers)
reconstructedODF[bin] += unitInc
return orientations, reconstructedODF
def TothVanHoutteSTAT (ODF,nSamples):
"""ODF contains 'dV_V' (normalized to 1), 'center', 'intervals', 'limits' (in radians)."""
orientations = np.zeros((nSamples,3),'f')
reconstructedODF = np.zeros(ODF['nBins'],'f')
unitInc = 1.0/nSamples
selectors = [random.random() for i in range(nSamples)]
selectors.sort()
indexSelector = 0
cumdV_V = 0.0
countSamples = 0
for bin in range(ODF['nBins']) :
cumdV_V += ODF['dV_V'][bin]
while indexSelector < nSamples and selectors[indexSelector] < cumdV_V:
Eulers = binAsEulers(bin,ODF['interval'],ODF['delta'],ODF['center'])
orientations[countSamples] = np.degrees(Eulers)
reconstructedODF[bin] += unitInc
countSamples += 1
indexSelector += 1
damask.util.croak('created set of %i when asked to deliver %i'%(countSamples,nSamples))
return orientations, reconstructedODF
# --------------------------------------------------------------------
# MAIN
# --------------------------------------------------------------------
parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', description ="""
Transform linear binned ODF data into given number of orientations.
IA: integral approximation, STAT: Van Houtte, MC: Monte Carlo
""", version = scriptID)
algorithms = ['IA', 'STAT','MC']
parser.add_option('-n', '--nsamples',
dest = 'number',
type = 'int', metavar = 'int',
help = 'number of orientations to be generated [%default]')
parser.add_option('-a','--algorithm',
dest = 'algorithm',
choices = algorithms, metavar = 'string',
help = 'sampling algorithm {%s} [IA]'%(', '.join(algorithms)))
parser.add_option('-p','--phase',
dest = 'phase',
type = 'int', metavar = 'int',
help = 'phase index to be used [%default]')
parser.add_option('-r', '--rnd',
dest = 'randomSeed',
type = 'int', metavar = 'int', \
help = 'seed of random number generator [%default]')
parser.set_defaults(randomSeed = None,
number = 500,
algorithm = 'IA',
phase = 1,
ang = True,
)
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
for name in filenames:
damask.util.report(scriptName,name)
table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name)
randomSeed = int(os.urandom(4).hex(),16) if options.randomSeed is None else options.randomSeed # random seed per file
random.seed(randomSeed)
# --------------- figure out limits (left/right), delta, and interval -----------------------------
ODF = {}
eulers = table.get('euler')
limits = np.array([np.min(eulers,axis=0),np.max(eulers,axis=0)]) # min/max euler angles in degrees
ODF['limit'] = np.radians(limits[1,:]) # right hand limits in radians
ODF['center'] = 0.0 if all(limits[0,:]<1e-8) else 0.5 # vertex or cell centered
ODF['interval'] = np.array(list(map(len,[np.unique(eulers[:,i]) for i in range(3)])),'i') # steps are number of distict values
ODF['nBins'] = ODF['interval'].prod()
ODF['delta'] = np.radians(np.array(limits[1,0:3]-limits[0,0:3])/(ODF['interval']-1)) # step size
if eulers.shape[0] != ODF['nBins']:
damask.util.croak('expecting %i values but got %i'%(ODF['nBins'],eulers.shape[0]))
continue
# ----- build binnedODF array and normalize ------------------------------------------------------
sumdV_V = 0.0
ODF['dV_V'] = [None]*ODF['nBins']
ODF['nNonZero'] = 0
dg = ODF['delta'][0]*2.0*math.sin(ODF['delta'][1]/2.0)*ODF['delta'][2]
intensity = table.get('intensity')
for b in range(ODF['nBins']):
ODF['dV_V'][b] = \
max(0.0,intensity[b,0]) * dg * \
math.sin(((b//ODF['interval'][2])%ODF['interval'][1]+ODF['center'])*ODF['delta'][1])
if ODF['dV_V'][b] > 0.0:
sumdV_V += ODF['dV_V'][b]
ODF['nNonZero'] += 1
for b in range(ODF['nBins']):
ODF['dV_V'][b] /= sumdV_V # normalize dV/V
damask.util.croak(['non-zero fraction: %12.11f (%i/%i)'%(float(ODF['nNonZero'])/ODF['nBins'],
ODF['nNonZero'],
ODF['nBins']),
'Volume integral of ODF: %12.11f\n'%sumdV_V,
'Reference Integral: %12.11f\n'%(ODF['limit'][0]*ODF['limit'][2]*(1-math.cos(ODF['limit'][1]))),
])
Functions = {'IA': 'directInversion', 'STAT': 'TothVanHoutteSTAT', 'MC': 'MonteCarloBins'}
method = Functions[options.algorithm]
Orientations, ReconstructedODF = (globals()[method])(ODF,options.number)
# calculate accuracy of sample
squaredDiff = {'orig':0.0,method:0.0}
squaredRelDiff = {'orig':0.0,method:0.0}
mutualProd = {'orig':0.0,method:0.0}
indivSum = {'orig':0.0,method:0.0}
indivSquaredSum = {'orig':0.0,method:0.0}
for bin in range(ODF['nBins']):
squaredDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[bin])**2
if ODF['dV_V'][bin] > 0.0:
squaredRelDiff[method] += (ODF['dV_V'][bin] - ReconstructedODF[bin])**2/ODF['dV_V'][bin]**2
mutualProd[method] += ODF['dV_V'][bin]*ReconstructedODF[bin]
indivSum[method] += ReconstructedODF[bin]
indivSquaredSum[method] += ReconstructedODF[bin]**2
indivSum['orig'] += ODF['dV_V'][bin]
indivSquaredSum['orig'] += ODF['dV_V'][bin]**2
damask.util.croak(['sqrt(N*)RMSD of ODFs:\t %12.11f'% math.sqrt(options.number*squaredDiff[method]),
'RMSrD of ODFs:\t %12.11f'%math.sqrt(squaredRelDiff[method]),
'rMSD of ODFs:\t %12.11f'%(squaredDiff[method]/indivSquaredSum['orig']),
'nNonZero correlation slope:\t %12.11f'\
%((ODF['nNonZero']*mutualProd[method]-indivSum['orig']*indivSum[method])/\
(ODF['nNonZero']*indivSquaredSum['orig']-indivSum['orig']**2)),
'nNonZero correlation confidence:\t %12.11f'\
%((mutualProd[method]-indivSum['orig']*indivSum[method]/ODF['nNonZero'])/\
(ODF['nNonZero']*math.sqrt((indivSquaredSum['orig']/ODF['nNonZero']-(indivSum['orig']/ODF['nNonZero'])**2)*\
(indivSquaredSum[method]/ODF['nNonZero']-(indivSum[method]/ODF['nNonZero'])**2)))),
])
if method == 'IA' and options.number < ODF['nNonZero']:
strOpt = '(%i)'%ODF['nNonZero']
formatwidth = 1+int(math.log10(options.number))
materialConfig = [
'#' + scriptID + ' ' + ' '.join(sys.argv[1:]),
'# random seed %i'%randomSeed,
'#-------------------#',
'<microstructure>',
'#-------------------#',
]
for i,ID in enumerate(range(options.number)):
materialConfig += ['[Grain%s]'%(str(ID+1).zfill(formatwidth)),
'(constituent) phase %i texture %s fraction 1.0'%(options.phase,str(ID+1).rjust(formatwidth)),
]
materialConfig += [
'#-------------------#',
'<texture>',
'#-------------------#',
]
for ID in range(options.number):
eulers = Orientations[ID]
materialConfig += ['[Grain%s]'%(str(ID+1).zfill(formatwidth)),
'(gauss) phi1 {} Phi {} phi2 {} scatter 0.0 fraction 1.0'.format(*eulers),
]
#--- output finalization --------------------------------------------------------------------------
with (open(os.path.splitext(name)[0]+'_'+method+'_'+str(options.number)+'_material.config','w')) as outfile:
outfile.write('\n'.join(materialConfig)+'\n')

26
processing/pre/seeds_fromDistribution.py
View File

@ -48,7 +48,6 @@ class myThread (threading.Thread):
myBestSeedsVFile = StringIO() # store local copy of best seeds file
perturbedSeedsVFile = StringIO() # perturbed best seeds file
perturbedGeomVFile = StringIO() # tessellated geom file
#--- still not matching desired bin class ----------------------------------------------------------
while bestMatch < options.threshold:
@ -92,15 +91,10 @@ class myThread (threading.Thread):
perturbedSeedsTable.set('pos',coords).save(perturbedSeedsVFile,legacy=True)
#--- do tesselation with perturbed seed file ------------------------------------------------------
perturbedGeomVFile.close()
perturbedGeomVFile = StringIO()
perturbedSeedsVFile.seek(0)
perturbedGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
' -g '+' '.join(list(map(str, options.grid))),streamIn=perturbedSeedsVFile)[0])
perturbedGeomVFile.seek(0)
perturbedGeom = damask.Geom.from_Voronoi_tessellation(options.grid,np.ones(3),coords)
#--- evaluate current seeds file ------------------------------------------------------------------
perturbedGeom = damask.Geom.load_ASCII(perturbedGeomVFile)
myNmaterials = len(np.unique(perturbedGeom.material))
currentData = np.bincount(perturbedGeom.material.ravel())[1:]/points
currentError=[]
@ -227,22 +221,15 @@ for i in range(1,nMaterials+1):
# ----------- create initial seed file or open existing one
bestSeedsVFile = StringIO()
if os.path.isfile(os.path.splitext(options.seedFile)[0]+'.seeds'):
with open(os.path.splitext(options.seedFile)[0]+'.seeds') as initialSeedFile:
for line in initialSeedFile: bestSeedsVFile.write(line)
initial_seeds = damask.Table.load(os.path.splitext(options.seedFile)[0]+'.seeds').get('pos')
else:
bestSeedsVFile.write(damask.util.execute('seeds_fromRandom'+\
' -g '+' '.join(list(map(str, options.grid)))+\
' -r {:d}'.format(options.randomSeed)+\
' -N '+str(nMaterials))[0])
initial_seeds = damask.seeds.from_random(np.ones(3),nMaterials,options.grid,options.randomSeed)
bestSeedsUpdate = time.time()
# ----------- tessellate initial seed file to get and evaluate geom file
bestSeedsVFile.seek(0)
initialGeomVFile = StringIO()
initialGeomVFile.write(damask.util.execute('geom_fromVoronoiTessellation '+
' -g '+' '.join(list(map(str, options.grid))),bestSeedsVFile)[0])
initialGeomVFile.seek(0)
initialGeom = damask.Geom.load_ASCII(initialGeomVFile)
initialGeom = damask.Geom.from_Voronoi_tessellation(options.grid,np.ones(3),initial_seeds)
if len(np.unique(targetGeom.material)) != nMaterials:
damask.util.croak('error. Material count mismatch')
@ -269,7 +256,6 @@ else:
if match >0: damask.util.croak('Stage {:d} cleared'.format(match))
sys.stdout.flush()
initialGeomVFile.close()
# start mulithreaded monte carlo simulation
threads = []

68
processing/pre/seeds_fromGeom.py
View File

@ -1,68 +0,0 @@
#!/usr/bin/env python3
import os
import sys
from io import StringIO
from optparse import OptionParser
import numpy as np
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 options [file[s]]', description = """
Create seed file taking material indices from given geom file.
Indices can be black-listed or white-listed.
""", version = scriptID)
parser.add_option('-w',
'--white',
action = 'extend', metavar = '<int LIST>',
dest = 'whitelist',
help = 'whitelist of grain IDs')
parser.add_option('-b',
'--black',
action = 'extend', metavar = '<int LIST>',
dest = 'blacklist',
help = 'blacklist of grain IDs')
parser.set_defaults(whitelist = [],
blacklist = [],
)
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
options.whitelist = [int(i) for i in options.whitelist]
options.blacklist = [int(i) for i in options.blacklist]
for name in filenames:
damask.util.report(scriptName,name)
geom = damask.Geom.load_ASCII(StringIO(''.join(sys.stdin.read())) if name is None else name)
material = geom.material.reshape((-1,1),order='F')
mask = np.logical_and(np.in1d(material,options.whitelist,invert=False) if options.whitelist else \
np.full(geom.grid.prod(),True,dtype=bool),
np.in1d(material,options.blacklist,invert=True) if options.blacklist else \
np.full(geom.grid.prod(),True,dtype=bool))
seeds = damask.grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
comments = geom.comments \
+ [scriptID + ' ' + ' '.join(sys.argv[1:]),
'grid\ta {}\tb {}\tc {}'.format(*geom.grid),
'size\tx {}\ty {}\tz {}'.format(*geom.size),
'origin\tx {}\ty {}\tz {}'.format(*geom.origin),
]
damask.Table(seeds[mask],{'pos':(3,)},comments)\
.add('material',material[mask].astype(int))\
.save(sys.stdout if name is None else os.path.splitext(name)[0]+'.seeds',legacy=True)

165
processing/pre/seeds_fromRandom.py
View File

@ -1,165 +0,0 @@
#!/usr/bin/env python3
import os
import sys
from optparse import OptionParser,OptionGroup
import numpy as np
from scipy import spatial
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 options', description = """
Distribute given number of points randomly within rectangular cuboid.
Reports positions with random crystal orientations in seeds file format to STDOUT.
""", version = scriptID)
parser.add_option('-N',
dest = 'N',
type = 'int', metavar = 'int',
help = 'number of seed points [%default]')
parser.add_option('-s',
'--size',
dest = 'size',
type = 'float', nargs = 3, metavar = 'float float float',
help='size x,y,z of unit cube to fill %default')
parser.add_option('-g',
'--grid',
dest = 'grid',
type = 'int', nargs = 3, metavar = 'int int int',
help='min a,b,c grid of hexahedral box %default')
parser.add_option('-m',
'--microstructure',
dest = 'microstructure',
type = 'int', metavar = 'int',
help = 'first microstructure index [%default]')
parser.add_option('-r',
'--rnd',
dest = 'randomSeed', type = 'int', metavar = 'int',
help = 'seed of random number generator [%default]')
group = OptionGroup(parser, "Laguerre Tessellation",
"Parameters determining shape of weight distribution of seed points"
)
group.add_option( '-w',
'--weights',
action = 'store_true',
dest = 'weights',
help = 'assign random weights to seed points for Laguerre tessellation [%default]')
group.add_option( '--max',
dest = 'max',
type = 'float', metavar = 'float',
help = 'max of uniform distribution for weights [%default]')
group.add_option( '--mean',
dest = 'mean',
type = 'float', metavar = 'float',
help = 'mean of normal distribution for weights [%default]')
group.add_option( '--sigma',
dest = 'sigma',
type = 'float', metavar = 'float',
help='standard deviation of normal distribution for weights [%default]')
parser.add_option_group(group)
group = OptionGroup(parser, "Selective Seeding",
"More uniform distribution of seed points using Mitchell's Best Candidate Algorithm"
)
group.add_option( '--selective',
action = 'store_true',
dest = 'selective',
help = 'selective picking of seed points from random seed points')
group.add_option( '--distance',
dest = 'distance',
type = 'float', metavar = 'float',
help = 'minimum distance to next neighbor [%default]')
group.add_option( '--numCandidates',
dest = 'numCandidates',
type = 'int', metavar = 'int',
help = 'size of point group to select best distance from [%default]')
parser.add_option_group(group)
parser.set_defaults(randomSeed = None,
grid = (16,16,16),
size = (1.0,1.0,1.0),
N = 20,
weights = False,
max = 0.0,
mean = 0.2,
sigma = 0.05,
microstructure = 1,
selective = False,
distance = 0.2,
numCandidates = 10,
)
(options,filenames) = parser.parse_args()
if filenames == []: filenames = [None]
size = np.array(options.size)
grid = np.array(options.grid)
np.random.seed(int(os.urandom(4).hex(),16) if options.randomSeed is None else options.randomSeed)
for name in filenames:
damask.util.report(scriptName,name)
if options.N > np.prod(grid):
damask.util.croak('More seeds than grid positions.')
sys.exit()
if options.selective and options.distance < min(size/grid):
damask.util.croak('Distance must be larger than grid spacing.')
sys.exit()
if options.selective and options.distance**3*options.N > 0.5*np.prod(size):
damask.util.croak('Number of seeds for given size and distance should be < {}.'\
.format(int(0.5*np.prod(size)/options.distance**3)))
eulers = np.random.rand(options.N,3) # create random Euler triplets
eulers[:,0] *= 360.0 # phi_1 is uniformly distributed
eulers[:,1] = np.degrees(np.arccos(2*eulers[:,1]-1.0)) # cos(Phi) is uniformly distributed
eulers[:,2] *= 360.0 # phi_2 is uniformly distributed
if not options.selective:
coords = damask.grid_filters.cell_coord0(grid,size).reshape(-1,3,order='F')
seeds = coords[np.random.choice(np.prod(grid), options.N, replace=False)] \
+ np.broadcast_to(size/grid,(options.N,3))*(np.random.rand(options.N,3)*.5-.25) # wobble without leaving grid
else:
seeds = np.empty((options.N,3))
seeds[0] = np.random.random(3) * size
i = 1
progress = damask.util._ProgressBar(options.N,'',50)
while i < options.N:
candidates = np.random.rand(options.numCandidates,3)*np.broadcast_to(size,(options.numCandidates,3))
tree = spatial.cKDTree(seeds[:i])
distances, dev_null = tree.query(candidates)
best = distances.argmax()
if distances[best] > options.distance: # require minimum separation
seeds[i] = candidates[best] # maximum separation to existing point cloud
i += 1
progress.update(i)
comments = [scriptID + ' ' + ' '.join(sys.argv[1:]),
'grid\ta {}\tb {}\tc {}'.format(*grid),
'size\tx {}\ty {}\tz {}'.format(*size),
'randomSeed\t{}'.format(options.randomSeed),
]
table = damask.Table(np.hstack((seeds,eulers)),{'pos':(3,),'euler':(3,)},comments)\
.add('microstructure',np.arange(options.microstructure,options.microstructure + options.N,dtype=int))
if options.weights:
weights = np.random.uniform(low = 0, high = options.max, size = options.N) if options.max > 0.0 \
else np.random.normal(loc = options.mean, scale = options.sigma, size = options.N)
table = table.add('weight',weights)
table.save(sys.stdout if name is None else name,legacy=True)

2
python/damask/seeds.py
View File

@ -95,7 +95,7 @@ def from_geom(geom,selection=None,invert=False,average=False,periodic=True):
"""
material = geom.material.reshape((-1,1),order='F')
mask = _np.full(geom.grid.prod(),True,dtype=bool) if selection is None else \
_np.isin(material,selection,invert=invert)
_np.isin(material,selection,invert=invert).flatten()
coords = grid_filters.cell_coord0(geom.grid,geom.size).reshape(-1,3,order='F')
if not average:

28
python/tests/test_seeds.py
View File

@ -3,6 +3,7 @@ import numpy as np
from scipy.spatial import cKDTree
from damask import seeds
from damask import grid_filters
from damask import Geom
class TestSeeds:
@ -25,7 +26,7 @@ class TestSeeds:
cKDTree(coords).query(coords, 2)
assert (0<= coords).all() and (coords<size).all() and np.min(min_dists[:,1])>=distance
def test_from_geom(self):
def test_from_geom_reconstruct(self):
grid = np.random.randint(10,20,3)
N_seeds = np.random.randint(30,300)
size = np.ones(3) + np.random.random(3)
@ -34,3 +35,28 @@ class TestSeeds:
coords,material = seeds.from_geom(geom_1)
geom_2 = Geom.from_Voronoi_tessellation(grid,size,coords,material)
assert (geom_2.material==geom_1.material).all()
@pytest.mark.parametrize('periodic',[True,False])
@pytest.mark.parametrize('average',[True,False])
def test_from_geom_grid(self,periodic,average):
grid = np.random.randint(10,20,3)
size = np.ones(3) + np.random.random(3)
coords = grid_filters.cell_coord0(grid,size).reshape(-1,3)
np.random.shuffle(coords)
geom_1 = Geom.from_Voronoi_tessellation(grid,size,coords)
coords,material = seeds.from_geom(geom_1,average=average,periodic=periodic)
geom_2 = Geom.from_Voronoi_tessellation(grid,size,coords,material)
assert (geom_2.material==geom_1.material).all()
@pytest.mark.parametrize('periodic',[True,False])
@pytest.mark.parametrize('average',[True,False])
@pytest.mark.parametrize('invert',[True,False])
def test_from_geom_selection(self,periodic,average,invert):
grid = np.random.randint(10,20,3)
N_seeds = np.random.randint(30,300)
size = np.ones(3) + np.random.random(3)
coords = seeds.from_random(size,N_seeds,grid)
geom = Geom.from_Voronoi_tessellation(grid,size,coords)
selection=np.random.randint(N_seeds)+1
coords,material = seeds.from_geom(geom,average=average,periodic=periodic,invert=invert,selection=[selection])
assert selection not in material if invert else (selection==material).all()