Index: /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.py (revision 18795) +++ /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.py (revision 18796) @@ -1,139 +1,156 @@ import numpy -def exportVTK(filename,model,options): - ''' - vtk export - function exportVTK(filename,model) - creates a directory with the vtk files for displays in paraview - (only work for triangle and wedges based on their number of nodes) - - Give only the results for nw but could be extended to geometry, mask... - - input: filename destination - (string) - ------------------------------------------------------------------ - model this is md - ------------------------------------------------------------------ - By default only the results are exported, you can add whichever - field you need as a string: - add 'geometry' to export md.geometry - - Basile de Fleurian: - ''' +import os +import model +import glob +def exportVTK(filename,model,*args): + ''' + vtk export + function exportVTK(filename,model) + creates a directory with the vtk files for displays in paraview + (only work for triangle and wedges based on their number of nodes) + + Give only the results for nw but could be extended to geometry, mask... + + input: filename destination + (string) + ------------------------------------------------------------------ +model this is md + ------------------------------------------------------------------ + By default only the results are exported, you can add whichever + field you need as a string: + add 'geometry' to export md.geometry -[path,name,ext]=fileparts(filename) -separator=filesep -mkdir(filename) + Basile de Fleurian: + ''' + Dir=os.path.basename(filename) + Path=filename[:-len(Dir)] -#get the element related variables -if model.mesh.shape[0]==2: - points=[model.mesh.x,model.mesh.y,zeros(model.mesh.numberofvertices,1)]; -else: - points=[model.mesh.x,model.mesh.y,model.mesh.z] + if os.path.exists(filename): + print ('File {} allready exist'.format(filename)) + newname=raw_input('Give a new name or "delete" to replace: ') + if newname=='delete': + filelist = glob.glob(filename+'/*') + for oldfile in filelist: + os.remove(oldfile) + else: + print ('New file name is {}'.format(newname)) + filename=newname + os.mkdir(filename) + else: + os.mkdir(filename) -[num_of_points,dim]=numpy.size(points) -[num_of_elt]=numpy.size(model.mesh.elements,1) -[point_per_elt]=numpy.size(model.mesh.elements,2) + #get the element related variables + if 'z' in dict.keys(model.mesh.__dict__): + points=numpy.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) + dim=3 + else: + points=numpy.column_stack((model.mesh.x,model.mesh.y,numpy.zeros(numpy.shape(model.mesh.x)))) + dim=2 -#Select the type of element function of the number of nodes per elements -if point_per_elt==3: - celltype=5 #triangles -elif point_per_elt==6: - celltype=13 #wedges -else: - error('Your Element definition is not taken into account \n') - -#this is the result structure -res_struct=model.results -if (len(fields(res_struct))>0): - #Getting all the solutions of the model - solnames=fields(res_struct) - num_of_sols=numpy.length(solnames) - num_of_timesteps=1 - #%building solutionstructure - for solution in num_of_sols: - sol_struct[i]=res_struct.(solnames[i]); - #looking for multiple time steps - if(numpy.size(sol_struct[i],2)>num_of_timesteps): - num_of_timesteps=numpy.size(sol_struct[i],2); + num_of_points=numpy.size(model.mesh.x) + num_of_elt=numpy.shape(model.mesh.elements)[0] + point_per_elt=numpy.shape(model.mesh.elements)[1] + + #Select the type of element function of the number of nodes per elements + if point_per_elt==3: + celltype=5 #triangles + elif point_per_elt==6: + celltype=13 #wedges + else: + error('Your Element definition is not taken into account \n') -else: - num_of_timesteps=1 + #this is the result structure + res_struct=model.results + if (len(res_struct.__dict__)>0): + #Getting all the solutions of the model + solnames=(dict.keys(res_struct.__dict__)) + num_of_sols=len(solnames) + num_of_timesteps=1 + out_freq=model.settings.output_frequency + #%building solutionstructure + for solution in solnames: + #looking for multiple time steps + if (len(res_struct.__dict__[solution])>num_of_timesteps): + num_of_timesteps=len(res_struct.__dict__[solution]) + num_of_timesteps=int(num_of_timesteps/out_freq)+1 + else: + num_of_timesteps=1 -for step in num_of_timesteps: - - timestep=step - - fid = open(strcat(path,filesep,name,filesep,name,'.vtk',int2str(timestep),'.vtk'),'w+') - fid.write('# vtk DataFile Version 2.0 \n') - fid.write('Data for run %s \n' % model.miscellaneous.name) - fid.write('ASCII \n') - fid.write('DATASET UNSTRUCTURED_GRID \n') + for step in range(0,num_of_timesteps): + timestep=step + fid=open((filename +'/Timestep.vtk'+str(timestep)+'.vtk'),'w+') + fid.write('# vtk DataFile Version 2.0 \n') + fid.write('Data for run %s \n' % model.miscellaneous.name) + fid.write('ASCII \n') + fid.write('DATASET UNSTRUCTURED_GRID \n') + fid.write('POINTS %d float\n' % num_of_points) + if(dim==3): + for point in points: + fid.write('%f %f %f \n'%(point[0], point[1], point[2])) + elif(dim==2): + for point in points: + fid.write('%f %f %f \n'%(point[0], point[1], point[2])) + + fid.write('CELLS %d %d\n' %(num_of_elt, num_of_elt*(point_per_elt+1))) + + if point_per_elt==3: + for elt in range(0, num_of_elt): + fid.write('3 %d %d %d\n' %(model.mesh.elements[elt,0]-1,model.mesh.elements[elt,1]-1,model.mesh.elements[elt,2]-1)) + elif point_per_elt==6: + for elt in range(0, num_of_elt): + fid.write('6 %d %d %d %d %d %d\n' %(model.mesh.elements[elt,0]-1,model.mesh.elements[elt,1]-1,model.mesh.elements[elt,2]-1,model.mesh.elements[elt,3]-1,model.mesh.elements[elt,4]-1,model.mesh.elements[elt,5]-1)) + else: + print 'Number of nodes per element not supported' + + fid.write('CELL_TYPES %d\n' %num_of_elt) + for elt in range(0, num_of_elt): + fid.write('%d\n' %celltype) + + fid.write('POINT_DATA %s \n' %str(num_of_points)) - fid.write(fid,'POINTS %d float\n',num_of_points) - if(dim==3): - s='%f %f %f \n' - elif(dim==2): - s='%f %f \n' - - P=[points zeros(num_of_points,3-dim)] - fid.write(fid,s,transpose de P) - - fid.write('CELLS %d %d\n' % num_of_elt % num_of_elt*(point_per_elt+1)) - s='%d' - for j=1:point_per_elt: - s=horzcat(s,{' %d'}) - - s=cell2mat(horzcat(s,{'\n'})) - fid.write(fid,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]transpose) - - fid.write(fid,'CELL_TYPES %d\n',num_of_elt) - s='%d\n' - fid.write(fid,s,celltype*ones(num_of_elt,1)) - fid.write(fid,'POINT_DATA %s \n',num2str(num_of_points)) - - #loop over the different solution structures - if 'num_of_sols' in locals(): - for j=1:num_of_sols: - #dealing with results on different timesteps - if(numpy.size(sol_struct{j},2)>timestep): - timestep = step - else: - timestep = numpy.size(sol_struct{j},2) - - #getting the number of fields in the solution - fieldnames=fields(sol_struct{j}(timestep)) - num_of_fields=numpy.length(fieldnames) - - #check which field is a real result and print - for k=1:num_of_fields: - if ((numel(sol_struct{j}(timestep).(fieldnames{k})))==num_of_points): - #paraview does not like NaN, replacing - nanval=find(isnan(sol_struct{j}(timestep).(fieldnames{k}))) - sol_struct{j}(timestep).(fieldnames{k})(nanval)=-9999 - #also checking for verry small value that mess up - smallval=(abs(sol_struct{j}(timestep).(fieldnames{k}))<1.0e-20) - sol_struct{j}(timestep).(fieldnames{k})(smallval)=0.0 - fid.write('SCALARS %s float 1 \n' % fieldnames{k}) - fid.write('LOOKUP_TABLE default\n') - s='%e\n' - fid.write(s % sol_struct{j}(timestep).(fieldnames{k})) - - #loop on arguments, if something other than result is asked, do - #it now - for j= 1:nargin-2: - res_struct=model.(varargin{j}) - fieldnames=fields(res_struct) - num_of_fields=numpy.length(fieldnames) - for k=1:num_of_fields: - if ((numel(res_struct.(fieldnames{k})))==num_of_points): - #paraview does not like NaN, replacing - nanval=find(isnan(res_struct.(fieldnames{k}))) - res_struct.(fieldnames{k})(nanval)=-9999 - #also checking for verry small value that mess up - smallval=(abs(res_struct.(fieldnames{k}))<1.0e-20) - res_struct.(fieldnames{k})(smallval)=0.0 - fid.write(fid,'SCALARS %s float 1 \n',fieldnames{k}) - fid.write(fid,'LOOKUP_TABLE default\n') - s='%e\n' - fid.write(fid,s,res_struct.(fieldnames{k})) - fid.close(); + #loop over the different solution structures + if 'solnames' in locals(): + for sol in solnames: + #dealing with results on different timesteps + if(len(res_struct.__dict__[sol])>timestep): + timestep = step + else: + timestep = len(res_struct.__dict__[sol]) + + #getting the fields in the solution + fieldnames=dict.keys(res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__) + + #check which field is a real result and print + for field in fieldnames: + if ((numpy.size(res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__[field]))==num_of_points): + fid.write('SCALARS %s float 1 \n' % field) + fid.write('LOOKUP_TABLE default\n') + for node in range(0,num_of_points): + #paraview does not like NaN, replacing + if numpy.isnan(res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__[field][node]): + fid.write('%e\n' % -9999.9999) + #also checking for verry small value that mess up + elif (abs(res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__[field][node])<1.0e-20): + fid.write('%e\n' % 0.0) + else: + fid.write('%e\n' % res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__[field][node]) + + #loop on arguments, if something other than result is asked, do + #it now + for other in args: + other_struct=model.__dict__[other] + othernames=(dict.keys(other_struct.__dict__)) + for field in othernames: + if ((numpy.size(other_struct.__dict__[field]))==num_of_points): + fid.write('SCALARS %s float 1 \n' % field) + fid.write('LOOKUP_TABLE default\n') + for node in range(0,num_of_points): + #paraview does not like NaN, replacing + if numpy.isnan(other_struct.__dict__[field][node]): + fid.write('%e\n' % -9999.9999) + #also checking for verry small value that mess up + elif (abs(other_struct.__dict__[field][node])<1.0e-20): + fid.write('%e\n' % 0.0) + else: + fid.write('%e\n' % other_struct.__dict__[field][node]) + fid.close();