1#!/usr/bin/env python 2import sys 3import os 4import vtk 5import numpy 6 7from siconos.io.mechanics_hdf5 import MechanicsHdf5 8 9## the best way to dump all data 10# $ h5dump toto.hdf5 > toto.txt 11 12import h5py 13 14import getopt 15 16def usage(): 17 """ 18 {0} <hdf5 file> 19 """.format(sys.argv[0]) 20 print '{0}: Usage'.format(sys.argv[0]) 21 print """ 22 {0} [--help] [--output_frequency=n] [--output_filename=] <hdf5 file> 23 """ 24 25try: 26 opts, args = getopt.gnu_getopt(sys.argv[1:], '', 27 ['output_frequency=', 28 'output_filename=', 29 'cf-scale=']) 30except getopt.GetoptError, err: 31 sys.stderr.write('{0}\n'.format(str(err))) 32 usage() 33 exit(2) 34 35min_time = None 36max_time = None 37scale_factor = 1 38output_frequency=10 39output_filename=None 40for o, a in opts: 41 42 if o == '--help': 43 usage() 44 exit(0) 45 46 elif o == '--output_frequency': 47 output_frequency = float(a) 48 49 elif o == '--output_filename': 50 out_filename=a 51 52 elif o == '--cf-scale': 53 scale_factor = float(a) 54 55 56if len(args) > 0: 57 io_filename = args[0] 58 if output_filename == None: 59 out_filename=''.join(io_filename.rsplit('.')[:-1])+'-filtered.hdf5' 60else: 61 usage() 62 exit(1) 63 64 65 66with MechanicsHdf5(io_filename=io_filename, mode='r') as io: 67 with MechanicsHdf5(io_filename=out_filename, mode='w') as out: 68 69 70 71 hdf1 = io._out 72 hdf2 = out._out 73 74 # copy /data/input 75 hdf2.__delitem__('data/input') 76 h5py.h5o.copy(hdf1.id, "data/input", hdf2.id, "data/input") 77 # copy /data/nslaws 78 hdf2.__delitem__('data/nslaws') 79 h5py.h5o.copy(hdf1.id, "data/nslaws", hdf2.id, "data/nslaws") 80 # copy /data/ref 81 hdf2.__delitem__('data/ref') 82 h5py.h5o.copy(hdf1.id, "data/ref", hdf2.id, "data/ref") 83 84 print('***************************************************** ') 85 print('************ Parsing simulation data ****************') 86 print('***************************************************** ') 87 88 def load(): 89 90 ispos_data = io.static_data() 91 idpos_data = io.dynamic_data() 92 93 icf_data = io.contact_forces_data()[:] 94 95 isolv_data = io.solver_data() 96 97 return ispos_data, idpos_data, icf_data, isolv_data 98 99 spos_data, dpos_data, cf_data, solv_data = load() 100 101 #print('io._data',io._data) 102 #print('static position data : spos_data',spos_data) 103 #print('spos_data.value',spos_data.value) 104 #print('dynamic position data : dpos_data',dpos_data) 105 print('dpos_data.value',dpos_data.value) 106 print('cf_data',cf_data) 107 #print('solv_data',solv_data) 108 109 times = list(set(dpos_data[:, 0])) 110 times.sort() 111 print('len(times)',len(times)) 112 113 if (len(times) ==0 ): 114 print('no results in the hdf5 file') 115 else: 116 117 118 print('Results for ',len(times),' steps in the hdf5 file') 119 #ndyna = len(numpy.where(dpos_data[:, 0] == times[0])) does not work 120 ndyna = len(dpos_data[:, 0])/len(times) 121 print('ndyna =', ndyna) 122 if len(spos_data) > 0: 123 nstatic = len(numpy.where(spos_data[:, 0] == times[0])) 124 nstatic = spos_data.shape[0] 125 else: 126 nstatic = 0 127 print('nstatic =', nstatic) 128 # instances = set(dpos_data[:, 1]) 129 130 131 132 # filtering 133 134 135 p=0 136 current_line=0 137 for k in range(len(times)): 138 #print(times[k]) 139 if (k+1 < len(times) ): 140 time_step=times[k+1]-times[k] 141 #print('time_step',time_step) 142 143 if (k%output_frequency==0): 144 if k==0 : 145 print('filter for k',k,'at times', times[k], 'p', p) 146 out._dynamic_data.resize((p+1)*ndyna,0) 147 out._dynamic_data[p*ndyna:(p+1)*ndyna,:] = numpy.array(dpos_data[k*ndyna:(k+1)*ndyna,:]) 148 #print('times',dpos_data[k*ndyna:(k+1)*ndyna,0]) 149 150 out._static_data.resize((p+1)*nstatic,0) 151 out._static_data[p*nstatic:(p+1)*nstatic,:] = numpy.array(spos_data[0:nstatic,:]) 152 153 out._solv_data.resize((p+1),0) 154 out._solv_data[p:(p+1),:] = numpy.array(solv_data[0:1,:]) 155 156 157 id_f = numpy.where(abs(cf_data[:, 0] - times[k]) < time_step*1e-5)[0] 158 if len(id_f) == 0: 159 print('no contact data at time',times[k]) 160 else: 161 #print('index of contact :', min(id_f), max(id_f)) 162 out._cf_data.resize(max(id_f)+1,0) 163 out._cf_data[min(id_f):max(id_f),:] = numpy.array(cf_data[min(id_f):max(id_f),:]) 164 current_line = max(id_f) 165 else: 166 print('filter for k',k,'at times', times[k], 'p', p) 167 out._dynamic_data.resize((p+1)*ndyna,0) 168 #print( dpos_data[k*ndyna:(k+1)*ndyna,:]) 169 #print('times',dpos_data[(k+1)*ndyna:(k+2)*ndyna,0]) 170 out._dynamic_data[p*ndyna:(p+1)*ndyna,:] = dpos_data[(k+1)*ndyna:(k+2)*ndyna,:] 171 172 # out._static_data.resize((p+1)*nstatic,0) 173 # #print( dpos_data[k*nstatic:(k+1)*nstatic,:]) 174 # out._static_data[p*nstatic:(p+1)*nstatic,:] = spos_data[k*nstatic:(k+1)*nstatic,:] 175 176 out._solv_data.resize((p+1),0) 177 out._solv_data[p:(p+1),:] = numpy.array(solv_data[k:k+1,:]) 178 179 id_f = numpy.where(abs(cf_data[:, 0] - times[k]) < time_step*1e-5)[0] 180 if len(id_f) == 0: 181 print('no contact data at time',times[k]) 182 else: 183 #print('index of contact :', min(id_f), max(id_f)) 184 new_line = current_line+max(id_f)-min(id_f)+1 185 #print('new_line',new_line) 186 #print('current_line',current_line) 187 #print('size of contact data', max(id_f)-min(id_f)+1) 188 out._cf_data.resize(new_line,0) 189 #print('fill out._cf_data indices', current_line, new_line-1) 190 out._cf_data[current_line:new_line,:] = numpy.array(cf_data[min(id_f):max(id_f)+1,:]) 191 current_line=new_line 192 #print('current_line',current_line) 193 194 p = p+1 195 #print(dpos_data) 196 197 print(out._dynamic_data.shape) 198 print(out._static_data.shape) 199 print(out.static_data().value) 200 print(out._cf_data.shape) 201 print(out._solv_data.shape) 202