DAMASK_EICMD/processing/pre/geom_fromBarycentric.py

#!/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 sys, os, string
import numpy as np
import argparse
from scipy.spatial import Delaunay

scriptID   = string.replace('$Id$','\n','\\n')
scriptName = os.path.splitext(scriptID.split()[1])[0]


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)