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

using class 

still a very complex script
This commit is contained in:
Martin Diehl 2019-05-26 22:49:05 +02:00
parent 0da39b0c69
commit c8dfba89e5
2 changed files with 76 additions and 102 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

172
processing/pre/geom_fromVoronoiTessellation.py
View File

@ -1,11 +1,14 @@
#!/usr/bin/env python3 #!/usr/bin/env python3
# -*- coding: UTF-8 no BOM -*-
import os,sys,math import os
import numpy as np import sys
import multiprocessing import multiprocessing
from io import StringIO
from optparse import OptionParser,OptionGroup from optparse import OptionParser,OptionGroup
from scipy import spatial
import numpy as np
from scipy import spatial
import damask import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
@ -29,74 +32,74 @@ def meshgrid2(*arrs):
ans.insert(0,arr2) ans.insert(0,arr2)
return tuple(ans) return tuple(ans)
def findClosestSeed(fargs):
def laguerreTessellation(undeformed, coords, weights, grains, periodic = True, cpus = 2):
def findClosestSeed(fargs):
point, seeds, myWeights = fargs point, seeds, myWeights = fargs
tmp = np.repeat(point.reshape(3,1), len(seeds), axis=1).T tmp = np.repeat(point.reshape(3,1), len(seeds), axis=1).T
dist = np.sum((tmp - seeds)**2,axis=1) -myWeights dist = np.sum((tmp - seeds)**2,axis=1) -myWeights
return np.argmin(dist) # seed point closest to point return np.argmin(dist) # seed point closest to point
copies = \
np.array([
[ -1,-1,-1 ],
[ 0,-1,-1 ],
[ 1,-1,-1 ],
[ -1, 0,-1 ],
[ 0, 0,-1 ],
[ 1, 0,-1 ],
[ -1, 1,-1 ],
[ 0, 1,-1 ],
[ 1, 1,-1 ],
[ -1,-1, 0 ],
[ 0,-1, 0 ],
[ 1,-1, 0 ],
[ -1, 0, 0 ],
[ 0, 0, 0 ],
[ 1, 0, 0 ],
[ -1, 1, 0 ],
[ 0, 1, 0 ],
[ 1, 1, 0 ],
[ -1,-1, 1 ],
[ 0,-1, 1 ],
[ 1,-1, 1 ],
[ -1, 0, 1 ],
[ 0, 0, 1 ],
[ 1, 0, 1 ],
[ -1, 1, 1 ],
[ 0, 1, 1 ],
[ 1, 1, 1 ],
]).astype(float)*info['size'] if periodic else \
np.array([
[ 0, 0, 0 ],
]).astype(float)
def laguerreTessellation(undeformed, coords, weights, grains, periodic = True, cpus = 2): repeatweights = np.tile(weights,len(copies)).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3)
for i,vec in enumerate(copies): # periodic copies of seed points ...
try: seeds = np.append(seeds, coords+vec, axis=0) # ... (1+a,2+a,3+a,...,1+z,2+z,3+z)
except NameError: seeds = coords+vec
copies = \ if (repeatweights == 0.0).all(): # standard Voronoi (no weights, KD tree)
np.array([ myKDTree = spatial.cKDTree(seeds)
[ -1,-1,-1 ], devNull,closestSeeds = myKDTree.query(undeformed)
[ 0,-1,-1 ], else:
[ 1,-1,-1 ], damask.util.croak('...using {} cpu{}'.format(options.cpus, 's' if options.cpus > 1 else ''))
[ -1, 0,-1 ], arguments = [[arg,seeds,repeatweights] for arg in list(undeformed)]
[ 0, 0,-1 ],
[ 1, 0,-1 ],
[ -1, 1,-1 ],
[ 0, 1,-1 ],
[ 1, 1,-1 ],
[ -1,-1, 0 ],
[ 0,-1, 0 ],
[ 1,-1, 0 ],
[ -1, 0, 0 ],
[ 0, 0, 0 ],
[ 1, 0, 0 ],
[ -1, 1, 0 ],
[ 0, 1, 0 ],
[ 1, 1, 0 ],
[ -1,-1, 1 ],
[ 0,-1, 1 ],
[ 1,-1, 1 ],
[ -1, 0, 1 ],
[ 0, 0, 1 ],
[ 1, 0, 1 ],
[ -1, 1, 1 ],
[ 0, 1, 1 ],
[ 1, 1, 1 ],
]).astype(float)*info['size'] if periodic else \
np.array([
[ 0, 0, 0 ],
]).astype(float)
repeatweights = np.tile(weights,len(copies)).flatten(order='F') # Laguerre weights (1,2,3,1,2,3,...,1,2,3) if cpus > 1: # use multithreading
for i,vec in enumerate(copies): # periodic copies of seed points ... pool = multiprocessing.Pool(processes = cpus) # initialize workers
try: seeds = np.append(seeds, coords+vec, axis=0) # ... (1+a,2+a,3+a,...,1+z,2+z,3+z) result = pool.map_async(findClosestSeed, arguments) # evaluate function in parallel
except NameError: seeds = coords+vec pool.close()
pool.join()
if (repeatweights == 0.0).all(): # standard Voronoi (no weights, KD tree) closestSeeds = np.array(result.get()).flatten()
myKDTree = spatial.cKDTree(seeds)
devNull,closestSeeds = myKDTree.query(undeformed)
else: else:
damask.util.croak('...using {} cpu{}'.format(options.cpus, 's' if options.cpus > 1 else '')) closestSeeds = np.zeros(len(arguments),dtype='i')
arguments = [[arg,seeds,repeatweights] for arg in list(undeformed)] for i,arg in enumerate(arguments):
closestSeeds[i] = findClosestSeed(arg)
if cpus > 1: # use multithreading
pool = multiprocessing.Pool(processes = cpus) # initialize workers
result = pool.map_async(findClosestSeed, arguments) # evaluate function in parallel
pool.close()
pool.join()
closestSeeds = np.array(result.get()).flatten()
else:
closestSeeds = np.zeros(len(arguments),dtype='i')
for i,arg in enumerate(arguments):
closestSeeds[i] = findClosestSeed(arg)
# closestSeed is modulo number of original seed points (i.e. excluding periodic copies) # closestSeed is modulo number of original seed points (i.e. excluding periodic copies)
return grains[closestSeeds%coords.shape[0]] return grains[closestSeeds%coords.shape[0]]
# -------------------------------------------------------------------- # --------------------------------------------------------------------
# MAIN # MAIN
@ -220,10 +223,7 @@ parser.set_defaults(pos = 'pos',
if filenames == []: filenames = [None] if filenames == []: filenames = [None]
for name in filenames: for name in filenames:
try: table = damask.ASCIItable(name = name, table = damask.ASCIItable(name = name, readonly = True)
outname = os.path.splitext(name)[0]+'.geom' if name else name,
buffered = False)
except: continue
damask.util.report(scriptName,name) damask.util.report(scriptName,name)
# --- read header ---------------------------------------------------------------------------- # --- read header ----------------------------------------------------------------------------
@ -281,7 +281,7 @@ for name in filenames:
else table.data[:,table.label_indexrange(options.pos)] - info['origin'] else table.data[:,table.label_indexrange(options.pos)] - info['origin']
eulers = table.data[:,table.label_indexrange(options.eulers)] if hasEulers \ eulers = table.data[:,table.label_indexrange(options.eulers)] if hasEulers \
else np.zeros(3*len(coords)) else np.zeros(3*len(coords))
grains = table.data[:,table.label_indexrange(options.microstructure)].astype('i') if hasGrains \ grains = table.data[:,table.label_indexrange(options.microstructure)].astype(int) if hasGrains \
else 1+np.arange(len(coords)) else 1+np.arange(len(coords))
weights = table.data[:,table.label_indexrange(options.weight)] if hasWeights \ weights = table.data[:,table.label_indexrange(options.weight)] if hasWeights \
else np.zeros(len(coords)) else np.zeros(len(coords))
@ -299,20 +299,8 @@ for name in filenames:
grid = np.vstack(meshgrid2(x, y, z)).reshape(3,-1).T grid = np.vstack(meshgrid2(x, y, z)).reshape(3,-1).T
indices = laguerreTessellation(grid, coords, weights, grains, options.periodic, options.cpus) indices = laguerreTessellation(grid, coords, weights, grains, options.periodic, options.cpus)
# --- write header ------------------------------------------------------------------------
usedGrainIDs = np.intersect1d(grainIDs,indices)
info['microstructures'] = len(usedGrainIDs)
if info['homogenization'] == 0: info['homogenization'] = options.homogenization
damask.util.report_geom(info,['grid','size','origin','homogenization',])
damask.util.croak(['microstructures: {}{}'.format(info['microstructures'],
(' out of {}'.format(NgrainIDs) if NgrainIDs != info['microstructures'] else '')),
])
config_header = [] config_header = []
formatwidth = 1+int(math.log10(NgrainIDs)) formatwidth = 1+int(np.log10(NgrainIDs))
if options.config: if options.config:
config_header += ['<microstructure>'] config_header += ['<microstructure>']
@ -331,24 +319,10 @@ for name in filenames:
] + theAxes ] + theAxes
config_header += ['<!skip>'] config_header += ['<!skip>']
table.labels_clear() geom = damask.Geom(indices.reshape(info['grid'],order='F'),info['size'],options.homogenization,comments=config_header)
table.info_clear() damask.util.croak(geom)
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']),
"microstructures\t{}".format(info['microstructures']),
config_header,
])
table.head_write()
# --- write microstructure information ------------------------------------------------------------ if name is None:
sys.stdout.write(str(geom.show()))
table.data = indices.reshape(info['grid'][1]*info['grid'][2],info['grid'][0]) else:
table.data_writeArray('%%%ii'%(formatwidth),delimiter=' ') geom.to_file(os.path.splitext(name)[0]+'.geom')
#--- output finalization --------------------------------------------------------------------------
table.close()

2
python/damask/geom.py
View File

@ -12,7 +12,7 @@ class Geom():
if len(microstructure.shape) != 3: if len(microstructure.shape) != 3:
raise ValueError('Invalid microstructure shape {}'.format(*microstructure.shape)) raise ValueError('Invalid microstructure shape {}'.format(*microstructure.shape))
elif microstructure.dtype not in ['int','float']: elif microstructure.dtype not in [int,float]:
raise TypeError('Invalid data type {} for microstructure'.format(microstructure.dtype)) raise TypeError('Invalid data type {} for microstructure'.format(microstructure.dtype))
else: else:
self.microstructure = microstructure self.microstructure = microstructure