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