From 952ad4f8fec00552c2f100151b88e145368c7fd2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 8 Oct 2020 18:33:40 +0200 Subject: [PATCH] functionality for geom generation in python lib --- processing/pre/geom_fromDREAM3D.py | 88 ++---------------------------- processing/pre/geom_fromTable.py | 57 ++----------------- python/damask/_geom.py | 38 ++++++++++++- 3 files changed, 44 insertions(+), 139 deletions(-) diff --git a/processing/pre/geom_fromDREAM3D.py b/processing/pre/geom_fromDREAM3D.py index daf7d2ab9..54ed4ca0a 100755 --- a/processing/pre/geom_fromDREAM3D.py +++ b/processing/pre/geom_fromDREAM3D.py @@ -1,12 +1,8 @@ #!/usr/bin/env python3 import os -import sys from optparse import OptionParser -import h5py -import numpy as np - import damask @@ -64,88 +60,12 @@ parser.set_defaults(pointwise = 'CellData', if options.basegroup is None: parser.error('No base group selected') -rootDir ='DataContainers' - - if filenames == []: parser.error('no input file specified.') for name in filenames: - damask.util.report(scriptName,name) + damask.util.report(scriptName,name) - errors = [] + geom = damask.Geom.load_DREAM3D(name,options.basegroup,options.pointwise) + damask.util.croak(geom) - inFile = h5py.File(name, 'r') - group_geom = os.path.join(rootDir,options.basegroup,'_SIMPL_GEOMETRY') - try: - size = inFile[os.path.join(group_geom,'DIMENSIONS')][...] \ - * inFile[os.path.join(group_geom,'SPACING')][...] - grid = inFile[os.path.join(group_geom,'DIMENSIONS')][...] - origin = inFile[os.path.join(group_geom,'ORIGIN')][...] - except KeyError: - errors.append('Geometry data ({}) not found'.format(group_geom)) - - - group_pointwise = os.path.join(rootDir,options.basegroup,options.pointwise) - if options.average is None: - label = 'Point' - - dataset = os.path.join(group_pointwise,options.quaternion) - try: - quats = np.reshape(inFile[dataset][...],(np.product(grid),4)) - rot = [damask.Rotation.from_quaternion(q,True,P=+1) for q in quats] - except KeyError: - errors.append('Pointwise orientation (quaternion) data ({}) not readable'.format(dataset)) - - dataset = os.path.join(group_pointwise,options.phase) - try: - phase = np.reshape(inFile[dataset][...],(np.product(grid))) - except KeyError: - errors.append('Pointwise phase data ({}) not readable'.format(dataset)) - - microstructure = np.arange(1,np.product(grid)+1,dtype=int).reshape(grid,order='F') - - - else: - label = 'Grain' - - dataset = os.path.join(group_pointwise,options.microstructure) - try: - microstructure = np.transpose(inFile[dataset][...].reshape(grid[::-1]),(2,1,0)) # convert from C ordering - except KeyError: - errors.append('Link between pointwise and grain average data ({}) not readable'.format(dataset)) - - group_average = os.path.join(rootDir,options.basegroup,options.average) - - dataset = os.path.join(group_average,options.quaternion) - try: - rot = [damask.Rotation.from_quaternion(q,True,P=+1) for q in inFile[dataset][...][1:]] # skip first entry (unindexed) - except KeyError: - errors.append('Average orientation data ({}) not readable'.format(dataset)) - - dataset = os.path.join(group_average,options.phase) - try: - phase = [i[0] for i in inFile[dataset][...]][1:] # skip first entry (unindexed) - except KeyError: - errors.append('Average phase data ({}) not readable'.format(dataset)) - - if errors != []: - damask.util.croak(errors) - continue - - config_header = [''] - for i in range(np.nanmax(microstructure)): - config_header += ['[{}{}]'.format(label,i+1), - '(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*rot[i].as_Eulers(degrees = True)), - ] - config_header += [''] - for i in range(np.nanmax(microstructure)): - config_header += ['[{}{}]'.format(label,i+1), - '(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(phase[i],i+1), - ] - - header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\ - + config_header - geom = damask.Geom(microstructure,size,origin,comments=header) - damask.util.croak(geom) - - geom.save_ASCII(os.path.splitext(name)[0]+'.geom',compress=False) + geom.save_ASCII(os.path.splitext(name)[0]+'.geom',compress=False) diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py index a0de6a4c5..d3cb2d50d 100755 --- a/processing/pre/geom_fromTable.py +++ b/processing/pre/geom_fromTable.py @@ -2,11 +2,8 @@ import os import sys -from io import StringIO from optparse import OptionParser -import numpy as np - import damask @@ -40,64 +37,20 @@ parser.add_option('-q', '--quaternion', dest = 'quaternion', type = 'string', metavar='string', help = 'quaternion label') -parser.add_option('--axes', - dest = 'axes', - type = 'string', nargs = 3, metavar = ' '.join(['string']*3), - help = 'orientation coordinate frame in terms of position coordinate frame [+x +y +z]') - -parser.set_defaults(pos = 'pos', - ) +parser.set_defaults(pos= 'pos') (options,filenames) = parser.parse_args() if filenames == []: filenames = [None] -if np.sum([options.quaternion is not None, - options.microstructure is not None]) != 1: - parser.error('need either microstructure or quaternion (and optionally phase) as input.') -if options.microstructure is not None and options.phase is not None: - parser.error('need either microstructure or phase (and mandatory quaternion) as input.') -if options.axes is not None and not set(options.axes).issubset(set(['x','+x','-x','y','+y','-y','z','+z','-z'])): - parser.error('invalid axes {} {} {}.'.format(*options.axes)) - for name in filenames: damask.util.report(scriptName,name) - table = damask.Table.load(StringIO(''.join(sys.stdin.read())) if name is None else name) - table.sort_by(['{}_{}'.format(i,options.pos) for i in range(3,0,-1)]) # x fast, y slow - grid,size,origin = damask.grid_filters.cell_coord0_gridSizeOrigin(table.get(options.pos)) + labels = [] + for l in [options.quaternion,options.phase,options.microstructure]: + if l is not None: labels.append(l) - config_header = table.comments - - if options.microstructure: - microstructure = table.get(options.microstructure).reshape(grid,order='F') - - elif options.quaternion: - q = table.get(options.quaternion) - phase = table.get(options.phase).astype(int) if options.phase else \ - np.ones((table.data.shape[0],1),dtype=int) - - unique,unique_inverse = np.unique(np.hstack((q,phase)),return_inverse=True,axis=0) - microstructure = unique_inverse.reshape(grid,order='F') + 1 - - config_header = [''] - for i,data in enumerate(unique): - ori = damask.Rotation(data[0:4]) - config_header += ['[Grain{}]'.format(i+1), - '(gauss)\tphi1 {:.2f}\tPhi {:.2f}\tphi2 {:.2f}'.format(*ori.as_Eulers(degrees = True)), - ] - if options.axes is not None: config_header += ['axes\t{} {} {}'.format(*options.axes)] - - config_header += [''] - for i,data in enumerate(unique): - config_header += ['[Grain{}]'.format(i+1), - '(constituent)\tphase {}\ttexture {}\tfraction 1.0'.format(int(data[4]),i+1), - ] - - header = [scriptID + ' ' + ' '.join(sys.argv[1:])]\ - + config_header - geom = damask.Geom(microstructure,size,origin, - comments=header) + geom = damask.Geom.load_table(name,options.pos,labels) damask.util.croak(geom) geom.save_ASCII(sys.stdout if name is None else os.path.splitext(name)[0]+'.geom',compress=False) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 2ccbb1988..ec65154bb 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -1,15 +1,18 @@ import copy import multiprocessing as mp from functools import partial +from os import path import numpy as np +import h5py from scipy import ndimage,spatial from . import environment -from . import Rotation from . import VTK from . import util from . import grid_filters +from . import Table +from . import Rotation class Geom: @@ -353,7 +356,7 @@ class Geom: Number of periods per unit cell. Defaults to 1. materials : (int, int), optional Material IDs. Defaults to (1,2). - + Notes ----- The following triply-periodic minimal surfaces are implemented: @@ -397,6 +400,35 @@ class Geom: ) + @staticmethod + def load_DREAM3D(fname,base_group,point_data=None,material='FeatureIds'): + root_dir ='DataContainers' + f = h5py.File(fname, 'r') + g = path.join(root_dir,base_group,'_SIMPL_GEOMETRY') + size = f[path.join(g,'DIMENSIONS')][()] * f[path.join(g,'SPACING')][()] + grid = f[path.join(g,'DIMENSIONS')][()] + origin = f[path.join(g,'ORIGIN')][()] + group_pointwise = path.join(root_dir,base_group,point_data) + + ma = np.arange(1,np.product(grid)+1,dtype=int) if point_data is None else \ + np.reshape(f[path.join(group_pointwise,material)],grid.prod()) + + return Geom(ma.reshape(grid,order='F'),size,origin,util.execution_stamp('Geom','from_DREAM3D')) + + + @staticmethod + def load_table(fname,coordinates,labels): + table = Table.load(fname).sort_by([f'{i}_{coordinates}' for i in range(3,0,-1)]) + + grid,size,origin = grid_filters.cell_coord0_gridSizeOrigin(table.get(coordinates)) + + labels_ = [labels] if isinstance(labels,str) else labels + _,unique_inverse = np.unique(np.hstack([table.get(l) for l in labels_]),return_inverse=True,axis=0) + ma = unique_inverse.reshape(grid,order='F') + 1 + + return Geom(ma,size,origin,util.execution_stamp('Geom','from_table')) + + def save(self,fname,compress=True): """ Generates vtk rectilinear grid. @@ -536,7 +568,7 @@ class Geom: coords_rot = R.broadcast_to(tuple(self.grid))@coords with np.errstate(all='ignore'): - mask = np.where(np.sum(np.power(coords_rot/r,2.0**exponent),axis=-1) > 1.0,True,False) + mask = np.sum(np.power(coords_rot/r,2.0**np.array(exponent)),axis=-1) > 1.0 if periodic: # translate back to center mask = np.roll(mask,((c-np.ones(3)*.5)*self.grid).astype(int),(0,1,2))