From c311ac47ca74eddab55ba17d6ee0f9b5a9c258c2 Mon Sep 17 00:00:00 2001 From: chen Date: Wed, 12 Oct 2016 20:20:15 -0400 Subject: [PATCH] force using rectangular mesh Paraview cannot handle large number of polyvertices using xdmf, forcing a meshed structure to by pass the problem. --- lib/damask/DS_HDF5.xml | 177 ++++++++++++++++++---------------- lib/damask/h5table.py | 26 +++-- processing/post/ascii2hdf5.py | 54 +++++++++-- 3 files changed, 152 insertions(+), 105 deletions(-) diff --git a/lib/damask/DS_HDF5.xml b/lib/damask/DS_HDF5.xml index af3523bfb..1277ce8d2 100644 --- a/lib/damask/DS_HDF5.xml +++ b/lib/damask/DS_HDF5.xml @@ -1,5 +1,6 @@ + attr / @@ -12,174 +13,186 @@ / - + + + Scalar + /Geometry/Vx + Geometry + Vector along x define the spectral mesh + + + + Scalar + /Geometry/Vy + Geometry + Vector along y defines the spectral mesh + + + + Scalar + /Geometry/Vz + Geometry + Vector along z defines the spectral mesh + + - scalar - /ip - + Scalar + /Geometry/ip + Geometry - - scalar - /node - + Scalar + /Geometry/node + Geometry - - vector - /grain - + Scalar + /Geometry/grain + Geometry - - vector - /pos - + Vector + /Geometry/pos + Geometry - - integer - /elem - + Scalar + /Geometry/elem + Geometry - + - integer - /phase + Scalar + /Crystallite/phase Crystallite - - integer - /texture + Scalar + /Crystallite/texture Crystallite - - scalar - /volume + Scalar + /Crystallite/volume Crystallite - - vector - /orientation + Vector + /Crystallite/orientation Crystallite - - vector - /eulerangles + Vector + /Crystallite/eulerangles Crystallite Bunnge Euler angles in degrees - - vector - /grainrotation + Vector + /Crystallite/grainrotation Crystallite - - tensor - /f + Tensor + /Crystallite/f Crystallite deformation gradient (F) -

- tensor - /p + Tensor + /Crystallite/p Crystallite Pikola Kirkhoff stress -

+ + Tensor + /Crystallite/Cauchy + Crystallite + Cauchy stress tensor + + + + Tensor + /Crystallite/lnV + Crystallite + + + + + Scalar + /Crystallite/MisesCauchy + Crystallite + von Mises equivalent of Cauchy stress + + + + Scalar + /Crystallite/MiseslnV + Crystallite + left von Mises strain + + + - vector - /resistance_slip + Vector + /Constitutive/resistance_slip Constitutive - vector - /shearrate_slip + Vector + /Constitutive/shearrate_slip Constitutive - vector - /resolvedstress_slip + Vector + /Constitutive/resolvedstress_slip Constitutive - scalar - /totalshear + Scalar + /Constitutive/totalshear Constitutive - vector - /accumulatedshear_slip + Matrix + /Constitutive/accumulatedshear_slip Constitutive vector contains accumulated shear per slip system - - tensor - /Cauchy - Derived - Cauchy stress tensor - - - - tensor - /lnV - Derived - - - - - scalar - /MisesCauchy - Derived - von Mises equivalent of Cauchy stress - - - - scalar - /MiseslnV - Derived - left von Mises strain - + \ No newline at end of file diff --git a/lib/damask/h5table.py b/lib/damask/h5table.py index a393c4805..f0c3fd72f 100644 --- a/lib/damask/h5table.py +++ b/lib/damask/h5table.py @@ -23,12 +23,6 @@ except(NameError): unicode=str -# ------------------------------------------------------- # -# Singleton class for converting feature name to H5F path # -# ------------------------------------------------------- # -# NOTE: -# use simple function to mimic the singleton class in -# C++/Java def lables_to_path(label, dsXMLPath=None): """read the xml definition file and return the path.""" if dsXMLPath is None: @@ -41,23 +35,26 @@ def lables_to_path(label, dsXMLPath=None): "DS_HDF5.xml") # This current implementation requires that all variables # stay under the root node, the nesting is defined through the - # h5path. This could be improved easily with more advanced parsing - # using ET interface, but for now I can not see the benefits in doing - # so. + # h5path. + # Allow new derived data to be put under the root tree = ET.parse(dsXMLPath) - dataType = tree.find('{}/type'.format(label)).text - h5path = tree.find('{}/h5path'.format(label)).text + try: + dataType = tree.find('{}/type'.format(label)).text + h5path = tree.find('{}/h5path'.format(label)).text + except: + dataType = "Scalar" + h5path = "/{}".format(label) # just put it under root return (dataType, h5path) class H5Table(object): - """Interface class for h5py + """light weight interface class for h5py DESCRIPTION ----------- Interface/wrapper class for manipulating data in HDF5 with DAMASK specialized data structure. - -->Minimal API design. + -->try to maintain a minimal API design. PARAMETERS ---------- h5f_path: str @@ -114,7 +111,8 @@ class H5Table(object): dsXMLPath=self.dsXMLFile) with h5py.File(self.h5f_path, 'a') as h5f: h5f_dst = h5f[h5f_path] # get the handle for target dataset(table) - rst_data = h5f_dst.read_direct(np.zeros(h5f_dst.shape)) + rst_data = np.zeros(h5f_dst.shape) + h5f_dst.read_direct(rst_data) return rst_data def add_data(self, feature_name, dataset, cmd_log=None): diff --git a/processing/post/ascii2hdf5.py b/processing/post/ascii2hdf5.py index e2f1622e0..d8c9daaab 100755 --- a/processing/post/ascii2hdf5.py +++ b/processing/post/ascii2hdf5.py @@ -7,6 +7,11 @@ # new one to the system wide definition file # # /lib/damask/DS_HDF5.xml # # or specify your own when initializing HDF5 class # +# 2. Somehow the point cloud structure cannot be properly handled # +# by Xdmf, which is a descriptive wrapper for visualizing HDF5 # +# using Paraview. The current solution is using cell structured # +# HDF5 so that Xdmf can describe the data shape as a rectangular # +# mesh rather than polyvertex. # # ------------------------------------------------------------------- # import os @@ -14,13 +19,24 @@ import damask import numpy as np from optparse import OptionParser + scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -# -------------------------------------------------------------------- -# MAIN -# -------------------------------------------------------------------- +# ----- helper function ----- # +def get_rectMshVectors(xyz_array, posNum): + """take in a xyz array from rectangular mesh and figure out Vx, Vy, Vz""" + # need some improvement, and only works for rectangular grid + v = sorted(list(set(xyz_array[:, posNum]))) + v_interval = (v[2]+v[1])/2.0 - (v[1]+v[0])/2.0 + v_start = (v[1]+v[0])/2.0 - v_interval + v_end = (v[-1]+v[-2])/2.0 + v_interval + V = np.linspace(v_start, v_end, len(v)+1) + return V + + +# ----- MAIN ---- # desp_msg = "Convert DAMASK ascii table to HDF5 file" parser = OptionParser(option_class=damask.extendableOption, usage='%prog options [file[s]]', @@ -31,8 +47,14 @@ parser.add_option('-D', '--DefinitionFile', type = 'string', metavar = 'string', help = 'definition file for H5 data storage') +parser.add_option('-p', + '--pos', '--position', + dest = 'pos', + type = 'string', metavar = 'string', + help = 'label of coordinates [%default]') -parser.set_defaults(DefinitionFile='default') +parser.set_defaults(DefinitionFile='default', + pos='pos') (options,filenames) = parser.parse_args() @@ -54,17 +76,31 @@ labels_idx = [asciiTable.label_index(label) for label in labels] featuresDim = [labels_idx[i+1] - labels_idx[i] for i in xrange(len(labels)-1)] featuresDim.append(fullTable.shape[1] - labels_idx[-1]) +# ----- figure out size and grid ----- # +pos_idx = asciiTable.label_index('pos') +xyz_array = asciiTable.data[:, pos_idx:pos_idx+3] +Vx = get_rectMshVectors(xyz_array, 0) +Vy = get_rectMshVectors(xyz_array, 1) +Vz = get_rectMshVectors(xyz_array, 2) +# use the dimension of the rectangular grid to reshape all other data +mshGridDim = (len(Vx)-1, len(Vy)-1, len(Vz)-1) + # ----- create a new HDF5 file and save the data -----# # Will overwrite existing HDF5 file with the same name - h5f = damask.H5Table(filename.replace(".txt", ".h5"), new_file=True, dsXMLFile=defFile) # adding increment number as root level attributes h5f.add_attr('inc', incNum) +# add the mesh grid data now +h5f.add_data("Vx", Vx) +h5f.add_data("Vy", Vy) +h5f.add_data("Vz", Vz) + +# add the rest of data from table for fi in xrange(len(labels)): featureName = labels[fi] - if featureName == 'inc': continue + if 'inc' in featureName: continue # remove trouble maker "("" and ")" if "(" in featureName: featureName = featureName.replace("(", "") if ")" in featureName: featureName = featureName.replace(")", "") @@ -72,8 +108,8 @@ for fi in xrange(len(labels)): featureDim = featuresDim[fi] # grab the data hook dataset = fullTable[:, featureIdx:featureIdx+featureDim] - # reshape tensor to 3x3 --> just make life slightly easier - if dataset.shape[1] == 9: - dataset = dataset.reshape((dataset.shape[0], 3, 3)) + # mapping 2D data onto a 3D rectangular mesh to get 4D data + dataset = dataset.reshape((mshGridDim[0], mshGridDim[1], mshGridDim[2], + dataset.shape[1])) # write out data h5f.add_data(featureName, dataset)