added options for laguerrre tessellation

processing/pre/geom_fromVoronoiTessellation.py

@@ -30,6 +30,45 @@ def meshgrid2(*arrs):
     ans.insert(0,arr2)
   return tuple(ans)
 
+def laguerreTessellation(undeformed, coords):
+    bestdist = np.ones(len(undeformed)) * np.finfo('d').max
+    bestseed = np.zeros(len(undeformed))
+
+    for i,seed in enumerate(coords):
+      for copy in np.array([[1, 0, 0,    ],
+                            [0, 1, 0,    ],
+                            [0, 0, 1,    ],
+                            [-1, 0, 0,   ],
+                            [0, -1, 0,   ],
+                            [0, 0, -1,   ],
+                            [1, 1, 0,    ],
+                            [1, 0, 1,    ],
+                            [0, 1, 1,    ],
+                            [-1, 1, 0,   ],
+                            [-1, 0, 1,   ],
+                            [0, -1, 1,   ],
+                            [-1, -1, 0,  ],
+                            [-1, 0, -1,  ],
+                            [0, -1, -1,  ],
+                            [1, -1, 0,   ],
+                            [1, 0, -1,   ],
+                            [0, 1, -1,   ],
+                            [1, 1, 1,    ],
+                            [-1, 1, 1,   ],
+                            [1, -1, 1,   ],
+                            [1, 1, -1,   ],
+                            [-1, -1, -1, ],
+                            [1, -1, -1,  ],
+                            [-1, 1, -1,  ],
+                            [-1, -1, 1,  ]]).astype(float):
+
+        diff = undeformed - np.repeat((seed+info['size']*copy).reshape(3,1),len(undeformed),axis=1).T
+        dist = np.sum(diff*diff,axis=1) - weights[i]
+                      
+        bestseed = np.where(dist < bestdist, np.ones(len(undeformed))*(i+1),bestseed)
+        bestdist = np.where(dist < bestdist, dist,bestdist)
+    return bestseed 
+
 
 # --------------------------------------------------------------------
 #                                MAIN
@@ -68,6 +107,9 @@ parser.add_option('-c', '--configuration', dest='config', action='store_true',
                   help='output material configuration [%default]')
 parser.add_option('--secondphase', type='float', dest='secondphase', metavar= 'float',
                   help='volume fraction of randomly distribute second phase [%default]')
+parser.add_option('--laguerre', dest='laguerre', action='store_true',
+                  help='for weighted voronoi (Laguerre) tessellation [%default]')
+
 
 parser.set_defaults(grid   = (0,0,0))
 parser.set_defaults(size   = (0.0,0.0,0.0))
@@ -77,6 +119,7 @@ parser.set_defaults(phase          = 1)
 parser.set_defaults(crystallite    = 1)
 parser.set_defaults(secondphase    = 0.0)
 parser.set_defaults(config = False)
+parser.set_defaults(laguerre = False)
 
 (options,filenames) = parser.parse_args()
 
@@ -105,22 +148,29 @@ for file in files:
   table = damask.ASCIItable(file['input'],file['output'],buffered = False)
   table.head_read()
 
-  coordsCol = table.labels_index('1_coords')
+  labels = []
-  if coordsCol < 0: 
+  if np.any(table.labels_index(['1_coords','2_coords','3_coords'])) == -1:
-    coordsCol = table.labels_index('x')                                                            # try if file is in legacy format
+    parser.error("missing seed coordinate column")
-    if coordsCol < 0: 
+  else:
-      file['croak'].write('column 1_coords/x not found...\n')
+    labels += ['1_coords','2_coords','3_coords']
-      continue
 
-  eulerCol = table.labels_index('phi1')
   hasEulers = np.all(table.labels_index(['phi1','Phi','phi2'])) != -1
-  grainCol = table.labels_index('microstructure')
+  if hasEulers:
-  hasGrains = grainCol != -1
+    labels += ['phi1','Phi','phi2']
     
-  table.data_readArray()
+  hasGrains = table.labels_index('microstructure') != -1
-  coords = table.data[:,coordsCol:coordsCol+3]
+  if hasGrains:
-  eulers = table.data[:,eulerCol:eulerCol+3] if hasEulers else np.zeros(3*len(coords))
+    labels += ['microstructure']
-  grain = table.data[:,grainCol] if hasGrains else 1+np.arange(len(eulers))
+    
+  hasWeight = table.labels_index('weight') != -1
+  if hasWeight:
+    labels += ['weight']
+      
+  table.data_readArray(labels)
+  coords = table.data[:,table.labels_index(['1_coords','2_coords','3_coords'])]
+  eulers = table.data[:,table.labels_index(['phi1','Phi','phi2'])] if hasEulers else np.zeros(3*len(coords))
+  grain = table.data[:,table.labels_index('microstructure')] if hasGrains else 1+np.arange(len(coords))
+  weights = table.data[:,table.labels_index('weight')] if hasWeight else np.zeros(len(coords))
   grainIDs = np.unique(grain).astype('i')
 
 
@@ -179,7 +229,6 @@ for file in files:
     continue
 
 #--- prepare data ---------------------------------------------------------------------------------
-  coords = (coords*info['size']).T
   eulers = eulers.T
 
 #--- switch according to task ---------------------------------------------------------------------
@@ -208,6 +257,9 @@ for file in files:
     x = (np.arange(info['grid'][0])+0.5)*info['size'][0]/info['grid'][0]
     y = (np.arange(info['grid'][1])+0.5)*info['size'][1]/info['grid'][1]
     z = (np.arange(info['grid'][2])+0.5)*info['size'][2]/info['grid'][2]
+    
+    if options.laguerre == False :
+      coords = (coords*info['size']).T
       undeformed = np.vstack(map(np.ravel, meshgrid2(x, y, z)))
   
       file['croak'].write('tessellating...\n')
@@ -216,6 +268,9 @@ for file in files:
                 np.eye(3),\
                 undeformed,coords)//3**3 + 1                                                          # floor division to kill periodic images
       indices = grain[indices-1]
+    else :
+      undeformed = np.vstack(np.meshgrid(x, y, z)).reshape(3,-1).T
+      indices = laguerreTessellation(undeformed, coords)
 
     newInfo['microstructures'] = info['microstructures']
     for i in grainIDs:

processing/pre/seeds_fromRandom.py

@@ -25,14 +25,28 @@ parser.add_option('-g','--grid', dest='grid', type='int', nargs=3, metavar='int
                   help='min a,b,c grid of hexahedral box %default')
 parser.add_option('-r', '--rnd', dest='randomSeed', type='int', metavar='int', \
                   help='seed of random number generator [%default]')
+parser.add_option('-w', '--weights', dest='weights', action='store_true',
+                  help = 'assign random weigts (Gaussian Distribution) to seed points for laguerre tessellation [%default]')
+parser.add_option('--mean', dest='mean', type='float', metavar='float', \
+                  help='mean of Gaussian Distribution for weights [%default]')
+parser.add_option('--sigma', dest='sigma', type='float', metavar='float', \
+                  help='standard deviation of Gaussian Distribution for weights [%default]')
+      
+
 
 parser.set_defaults(randomSeed = None)
 parser.set_defaults(grid = (16,16,16))
 parser.set_defaults(N = 20)
+parser.set_defaults(weights=False)
+parser.set_defaults(mean = 0.0)
+parser.set_defaults(sigma = 1.0)
+
 
 (options,filename) = parser.parse_args()
 options.grid = np.array(options.grid)
 
+labels = "1_coords\t2_coords\t3_coords\tphi1\tPhi\tphi2"
+
 # ------------------------------------------ setup file handle -------------------------------------
 if filename == []:
   file = {'output':sys.stdout, 'croak':sys.stderr}
@@ -48,6 +62,8 @@ if options.N > gridSize:
   options.N = gridSize
 if options.randomSeed == None:
   options.randomSeed = int(os.urandom(4).encode('hex'), 16)
+
+
 np.random.seed(options.randomSeed)                                                      # init random generators
 random.seed(options.randomSeed)
 
@@ -76,14 +92,24 @@ seeds[1,:] = (np.mod(seedpoints//                 options.grid[0] ,options.grid[
 seeds[2,:] = (np.mod(seedpoints//(options.grid[1]*options.grid[0]),options.grid[2])\
              +np.random.random())/options.grid[2]
 
+table = np.transpose(np.concatenate((seeds,grainEuler),axis = 0))
+
+if options.weights :
+   weight = np.random.normal(loc=options.mean, scale=options.sigma, size=options.N)
+   weight /= np.sum(weight)
+   table = np.append(table, weight.reshape(options.N,1), axis=1)
+   labels += "\tweight"
+
+
+
 header = ["5\theader",
           scriptID + " " + " ".join(sys.argv[1:]),
           "grid\ta {}\tb {}\tc {}".format(options.grid[0],options.grid[1],options.grid[2]),
           "microstructures\t{}".format(options.N),
           "randomSeed\t{}".format(options.randomSeed),
-          "1_coords\t2_coords\t3_coords\tphi1\tPhi\tphi2",
+          "%s"%labels,
          ]
 
 for line in header:
   file['output'].write(line+"\n")
-np.savetxt(file['output'],np.transpose(np.concatenate((seeds,grainEuler),axis = 0)),fmt='%10.6f',delimiter='\t')
+np.savetxt(file['output'], table, fmt='%10.6f', delimiter='\t')