298 lines
12 KiB
Python
Executable File
298 lines
12 KiB
Python
Executable File
#!/usr/bin/env python
|
|
|
|
##
|
|
# This script will read in all the seeds and partition the space
|
|
# using scipy.spatial.Delaunay triangulation.
|
|
# The unknown location will be then interpolated through Barycentric
|
|
# interpolation method, which relies on the triangulation.
|
|
# A rim will be automatically added to the patch, which will help
|
|
# improve the compatibility with the spectral solver as well as
|
|
# maintain meaningful microstructure(reduce artifacts).
|
|
|
|
|
|
import os
|
|
import numpy as np
|
|
import argparse
|
|
from scipy.spatial import Delaunay
|
|
import damask
|
|
|
|
scriptName = os.path.splitext(os.path.basename(__file__))[0]
|
|
scriptID = ' '.join([scriptName,damask.version])
|
|
|
|
|
|
OFFSET = 0.1 #resize the seeded volume to give space for rim/pan
|
|
PHANTOM_ID = -1 #grain ID for phantom seeds
|
|
|
|
def d_print(info, data, separator=False):
|
|
"""quickly print debug information"""
|
|
if(separator): print "*"*80
|
|
print info
|
|
print data
|
|
|
|
|
|
def meshgrid2(*arrs):
|
|
"""code inspired by http://stackoverflow.com/questions/1827489/numpy-meshgrid-in-3d"""
|
|
arrs = tuple(reversed(arrs))
|
|
arrs = tuple(arrs)
|
|
lens = np.array(map(len, arrs))
|
|
dim = len(arrs)
|
|
ans = []
|
|
for i, arr in enumerate(arrs):
|
|
slc = np.ones(dim,'i')
|
|
slc[i] = lens[i]
|
|
arr2 = np.asarray(arr).reshape(slc)
|
|
for j, sz in enumerate(lens):
|
|
if j != i:
|
|
arr2 = arr2.repeat(sz, axis=j)
|
|
|
|
ans.insert(0,arr2)
|
|
return tuple(ans)
|
|
|
|
|
|
#prepare command line interface
|
|
parser = argparse.ArgumentParser(prog="geoFromBarycentic",
|
|
description='''Generate geom file through \
|
|
Barycentric interpolating seeds file.''',
|
|
epilog="requires numpy, and scipy.")
|
|
parser.add_argument("seeds",
|
|
help="seeds file in DAMASK format:\
|
|
http://damask.mpie.de/Documentation/AsciiTableFormat",
|
|
default="test.seeds")
|
|
parser.add_argument("-v", "--version",
|
|
action="version",
|
|
version="%(prog)s 0.1")
|
|
parser.add_argument("-g", "--grid",
|
|
nargs=3,
|
|
help="grid size(mesh resolution, recommend using 2^x)",
|
|
default=[32,32,32],
|
|
type=int)
|
|
parser.add_argument("-s", "--size",
|
|
help="physical size of the target volume.",
|
|
nargs=3,
|
|
default=[1.0,1.0,1.0],
|
|
type=float)
|
|
parser.add_argument("-o", "--origin",
|
|
help="lower left corner of the patch.",
|
|
nargs=3,
|
|
default=[0.0,0.0,0.0],
|
|
type=float)
|
|
parser.add_argument('-m', '--homogenization',
|
|
help='homogenization index to be used',
|
|
default=1,
|
|
type=int)
|
|
parser.add_argument('-c', '--crystallite',
|
|
help='crystallite index to be used',
|
|
default=1,
|
|
type=int)
|
|
parser.add_argument('-p', '--phase',
|
|
help='phase index to be used',
|
|
default=1,
|
|
type=int)
|
|
parser.add_argument('-F', '--Favg',
|
|
help='reshape the periodicity, not useful for RIM method',
|
|
nargs=9,
|
|
default=[1.0,0.0,0.0,
|
|
0.0,1.0,0.0,
|
|
0.0,0.0,1.0],
|
|
type=float)
|
|
parser.add_argument("-G", "--geomFile",
|
|
help='the name of the output geom file',
|
|
default='seeds.geom',
|
|
type=str)
|
|
parser.add_argument("-C", "--configFile",
|
|
help='output dummy material.config file',
|
|
action='store_true',
|
|
default=False)
|
|
parser.add_argument("-d", "--debug",
|
|
help="start debugging script",
|
|
action='store_true',
|
|
default=False)
|
|
parser.add_argument("-S", "--seedsFile",
|
|
help="write out resized seeds file",
|
|
action='store_true',
|
|
default=False)
|
|
parser.add_argument("-r", '--addRim',
|
|
help="add rim and provide control of face lifting point",
|
|
action='store_true',
|
|
default=False)
|
|
args = parser.parse_args() # get all the arguments right after
|
|
|
|
#quick help to user
|
|
print "*"*80
|
|
parser.print_help()
|
|
print """Sample usage:
|
|
./geoFromBarycentic.py 20grains.seeds -g 128 128 128 -S -r; geom_check seeds.geom; seeds_check new_seed.seeds.
|
|
"""
|
|
print "*"*80
|
|
if (args.debug):
|
|
d_print("args are:", parser.parse_args(),separator=True)
|
|
|
|
#/\/\/\/\/#
|
|
# m a i n #
|
|
#\/\/\/\/\#
|
|
print "only work for 3D case now, 2D support coming soon..."
|
|
print "reading seeds file: {}".format(args.seeds)
|
|
|
|
with open(args.seeds, 'r') as f:
|
|
rawtext = f.readlines()
|
|
n_header = int(rawtext.pop(0).split()[0])
|
|
#record all the seeds position
|
|
if (args.addRim):
|
|
grid_shift = np.array(args.size) * np.array([OFFSET,OFFSET,OFFSET*2])
|
|
s_coords = np.array([[np.array(float(item))*(1 - OFFSET*2)
|
|
for item in line.split()[:3]] + grid_shift
|
|
for line in rawtext[n_header:]])
|
|
else:
|
|
#no need for shifting with periodicity
|
|
s_coords = np.array([[np.array(float(item))
|
|
for item in line.split()[:3]]
|
|
for line in rawtext[n_header:]])
|
|
|
|
#record ID of the seeds: int/EulerAngles
|
|
if 'microstructure' in rawtext[n_header-1]:
|
|
s_id = [int(line.split()[-1]) for line in rawtext[n_header:]]
|
|
else:
|
|
print "WARNING:"
|
|
print "THIS SCRIPT DOES NOT UNDERSTAND HOW TO GROUP CRYSTALLITES."
|
|
print "ALL CRYSTAL ORIENTATIONS ARE CONSIDERED TO BE UNIQUE."
|
|
print "FOR MORE ACCURATE CONTROL OF SEEDS GROUPING, USE MICROSTRUCTURE ID."
|
|
s_id = range(len(s_coords))
|
|
#s_eulers here is just a quick book keeping
|
|
s_eulers = np.array([[float(item) for item in line.split()[3:]]
|
|
for line in rawtext[n_header:]])
|
|
|
|
if(args.debug):
|
|
print d_print("resize point cloud to make space for rim/pan:",
|
|
s_coords)
|
|
|
|
if(args.addRim):
|
|
#add binding box to create rim/pan for the volume where the ID of the seeds is
|
|
#unknown
|
|
print "Shrining the seeds to {}x in each direction".format(1 - OFFSET*2)
|
|
x,y,z = args.size[0],args.size[1],args.size[2]
|
|
print "Use circumscribed sphere to place phantom seeds."
|
|
r = np.sqrt(x**2+y**2+z**2)/2.0
|
|
BINDBOX = [[0,0,0],[x,0,0],[0,y,0],[x,y,0],
|
|
[0,0,z],[x,0,z],[0,y,z],[x,y,z],
|
|
[x/2.0+r,y/2, z/2], [x/2.0-r, y/2, z/2],
|
|
[x/2, y/2.0+r, z/2], [x/2, y/2.0-r, z/2],
|
|
[x/2, y/2, z/2.0-r]] #8 corners + 5 face centers (no top)
|
|
print "Adding phantom seeds for RIM generation:"
|
|
for point in BINDBOX:
|
|
print point
|
|
s_coords = np.vstack([s_coords,point])
|
|
s_id.append(PHANTOM_ID)
|
|
else:
|
|
#The idea here is that we read in each seed point, than duplicate in 3D (make a few copies),
|
|
#move on to the next seed point, repeat the same procedure. As for the ID list, we can just use the
|
|
#same one. The trick here is use the floor division to find the correct id since we pretty much duplicate
|
|
#the same point several times.
|
|
Favg = np.array(args.Favg).reshape((3,3))
|
|
x,y,z = args.size[0],args.size[1],args.size[2]
|
|
tmp = []
|
|
for seed in s_coords:
|
|
tmp += [np.dot(Favg, np.array(seed) + np.array([dx,dy,dz]))
|
|
for dz in [-z, 0, z]
|
|
for dy in [-y, 0, y]
|
|
for dx in [-x, 0, x]]
|
|
s_coords = tmp
|
|
|
|
if (args.seedsFile):
|
|
with open("new_seed.seeds", "w") as f:
|
|
outstr = "4\theader\n"
|
|
outstr += "grid\ta {}\tb {}\tc {}\n".format(args.grid[0],
|
|
args.grid[1],
|
|
args.grid[2])
|
|
outstr += "microstructures {}\n".format(len(set(s_id)))
|
|
outstr += "x\ty\tz\tmicrostructure"
|
|
if (args.addRim):
|
|
for i in range(len(s_id)):
|
|
outstr += "{}\t{}\t{}\t{}\n".format(s_coords[i][0],
|
|
s_coords[i][1],
|
|
s_coords[i][2],
|
|
s_id[i])
|
|
else:
|
|
for i in range(len(s_coords)):
|
|
outstr += "{}\t{}\t{}\t{}\n".format(s_coords[i][0],
|
|
s_coords[i][1],
|
|
s_coords[i][2],
|
|
s_id[i//3**3])
|
|
f.write(outstr)
|
|
|
|
#triangulate space with given point-clouds
|
|
tri = Delaunay(s_coords)
|
|
|
|
if(args.debug):
|
|
d_print("simplices:", tri.simplices, separator=True)
|
|
d_print("vertices:", s_coords[tri.simplices])
|
|
|
|
#populate grid points (only 3D for now)
|
|
'''
|
|
#populating grid using meshgrid2
|
|
x = (np.arange(args.grid[0])+0.5)*args.size[0]/args.grid[0]
|
|
y = (np.arange(args.grid[1])+0.5)*args.size[1]/args.grid[1]
|
|
z = (np.arange(args.grid[2])+0.5)*args.size[2]/args.grid[2]
|
|
mesh_pts = np.transpose(np.vstack(map(np.ravel, meshgrid2(x, y, z))))
|
|
print mesh_pts
|
|
'''
|
|
#this is actually faster than using meshgrid2
|
|
mesh_pts = [[(i+0.5)*args.size[0]/args.grid[0],
|
|
(j+0.5)*args.size[1]/args.grid[1],
|
|
(k+0.5)*args.size[2]/args.grid[2]]
|
|
for k in range(args.grid[2])
|
|
for j in range(args.grid[1])
|
|
for i in range(args.grid[0])]
|
|
|
|
mesh_ids = [PHANTOM_ID*2]*len(mesh_pts) #initialize grid
|
|
|
|
#search ID for each grid point
|
|
s_id = np.array(s_id) #allow multi-indexing
|
|
mesh_idx = tri.find_simplex(mesh_pts)
|
|
|
|
for i, pt in enumerate(mesh_pts):
|
|
if mesh_idx[i] < 0:
|
|
continue #didn't find any envelop tetrahedron --> something wrong!
|
|
#calculate Barycentric coordinates
|
|
bary_c = tri.transform[mesh_idx[i],:3,:3].dot(pt-tri.transform[mesh_idx[i],3,:])
|
|
bary_c = np.append(bary_c, 1 - bary_c.sum())
|
|
|
|
if (args.addRim):
|
|
tmp_ids = s_id[tri.simplices[mesh_idx[i]]] #rim method
|
|
else:
|
|
tmp_ids = np.array(s_id[tri.simplices[mesh_idx[i]]//(3**3)]) #kill periodicity through floor division
|
|
#print tmp_ids
|
|
#print tri.simplices[mesh_idx[i]]//(3**3)
|
|
|
|
max_weight = -1960
|
|
for this_id in tmp_ids:
|
|
msk = [item==this_id for item in tmp_ids] #find vertex with the same id
|
|
tmp_weight = sum([bary_c[j] for j in range(len(bary_c)) if msk[j]])
|
|
if tmp_weight > max_weight:
|
|
max_weight = tmp_weight
|
|
mesh_ids[i] = this_id
|
|
if (args.debug):
|
|
d_print("bary_c:",bary_c,separator=True)
|
|
d_print("vertex ID:", tmp_ids)
|
|
d_print("final ID:", mesh_ids[i])
|
|
|
|
mesh_ids = np.reshape(mesh_ids, (-1, args.grid[0]))
|
|
|
|
#write to file
|
|
with open(args.geomFile, "w") as f:
|
|
outstr = "5\theader\n"
|
|
outstr += "grid\ta {}\tb {}\tc {}\n".format(args.grid[0],
|
|
args.grid[1],
|
|
args.grid[2])
|
|
outstr += "size\tx {}\ty {}\tz {}\n".format(args.size[0],
|
|
args.size[1],
|
|
args.size[2])
|
|
outstr += "origin\tx {}\ty {}\tz {}\n".format(args.origin[0],
|
|
args.origin[1],
|
|
args.origin[2])
|
|
outstr += "homogenization\t{}\nmicrostructure\t{}\n".format(args.homogenization,
|
|
len(set(s_id)))
|
|
for row in mesh_ids:
|
|
row = [str(item) for item in list(row)]
|
|
outstr += "\t".join(row) + "\n"
|
|
f.write(outstr)
|