Changeset 18796


Ignore:
Timestamp:
11/17/14 13:21:48 (10 years ago)
Author:
bdef
Message:

BUG:ReWriting export paraview in python

File:
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/m/contrib/paraview/exportVTK.py

    r18516 r18796  
    11import numpy
    2 def exportVTK(filename,model,options):
    3     '''
    4     vtk export
    5     function exportVTK(filename,model)
    6     creates a directory with the vtk files for displays in paraview
    7     (only work for triangle and wedges based on their number of nodes)
    8    
    9     Give only the results for nw but could be extended to geometry, mask...
    10    
    11     input: filename   destination
    12     (string)
    13     ------------------------------------------------------------------
    14     model      this is md
    15     ------------------------------------------------------------------
    16     By default only the results are exported, you can add whichever
    17     field you need as a string:
    18     add 'geometry' to export md.geometry
    19    
    20     Basile de Fleurian:
    21     '''
     2import os
     3import model
     4import glob
     5def exportVTK(filename,model,*args):
     6        '''
     7        vtk export
     8        function exportVTK(filename,model)
     9        creates a directory with the vtk files for displays in paraview
     10        (only work for triangle and wedges based on their number of nodes)
     11       
     12        Give only the results for nw but could be extended to geometry, mask...
     13       
     14        input: filename   destination
     15        (string)
     16        ------------------------------------------------------------------
     17model      this is md
     18        ------------------------------------------------------------------
     19        By default only the results are exported, you can add whichever
     20        field you need as a string:
     21        add 'geometry' to export md.geometry
    2222
    23 [path,name,ext]=fileparts(filename)
    24 separator=filesep
    25 mkdir(filename)
     23        Basile de Fleurian:
     24        '''
     25        Dir=os.path.basename(filename)
     26        Path=filename[:-len(Dir)]
    2627
    27 #get the element related variables
    28 if model.mesh.shape[0]==2:
    29     points=[model.mesh.x,model.mesh.y,zeros(model.mesh.numberofvertices,1)];
    30 else:
    31     points=[model.mesh.x,model.mesh.y,model.mesh.z]
     28        if os.path.exists(filename):
     29                print ('File {} allready exist'.format(filename))
     30                newname=raw_input('Give a new name or "delete" to replace: ')
     31                if newname=='delete':
     32                        filelist = glob.glob(filename+'/*')
     33                        for oldfile in filelist:
     34                                os.remove(oldfile)
     35                else:
     36                        print ('New file name is {}'.format(newname))
     37                        filename=newname
     38                        os.mkdir(filename)
     39        else:
     40                os.mkdir(filename)
    3241
    33 [num_of_points,dim]=numpy.size(points)
    34 [num_of_elt]=numpy.size(model.mesh.elements,1)
    35 [point_per_elt]=numpy.size(model.mesh.elements,2)
     42        #get the element related variables
     43        if 'z' in dict.keys(model.mesh.__dict__):
     44                points=numpy.column_stack((model.mesh.x,model.mesh.y,model.mesh.z))
     45                dim=3
     46        else:
     47                points=numpy.column_stack((model.mesh.x,model.mesh.y,numpy.zeros(numpy.shape(model.mesh.x))))
     48                dim=2
    3649
    37 #Select the type of element function of the number of nodes per elements
    38 if point_per_elt==3:
    39     celltype=5 #triangles
    40 elif point_per_elt==6:
    41     celltype=13 #wedges
    42 else:
    43     error('Your Element definition is not taken into account \n')
    44    
    45 #this is the result structure
    46 res_struct=model.results
    47 if (len(fields(res_struct))>0):
    48     #Getting all the solutions of the model
    49     solnames=fields(res_struct)
    50     num_of_sols=numpy.length(solnames)
    51     num_of_timesteps=1
    52     #%building solutionstructure
    53     for solution in num_of_sols:
    54         sol_struct[i]=res_struct.(solnames[i]);
    55         #looking for multiple time steps
    56         if(numpy.size(sol_struct[i],2)>num_of_timesteps):
    57             num_of_timesteps=numpy.size(sol_struct[i],2);
     50        num_of_points=numpy.size(model.mesh.x)
     51        num_of_elt=numpy.shape(model.mesh.elements)[0]
     52        point_per_elt=numpy.shape(model.mesh.elements)[1]
     53               
     54        #Select the type of element function of the number of nodes per elements
     55        if point_per_elt==3:
     56                celltype=5 #triangles
     57        elif point_per_elt==6:
     58                celltype=13 #wedges
     59        else:
     60                error('Your Element definition is not taken into account \n')
    5861
    59 else:
    60     num_of_timesteps=1
     62        #this is the result structure
     63        res_struct=model.results
     64        if (len(res_struct.__dict__)>0):
     65                #Getting all the solutions of the model
     66                solnames=(dict.keys(res_struct.__dict__))
     67                num_of_sols=len(solnames)
     68                num_of_timesteps=1
     69                out_freq=model.settings.output_frequency
     70                #%building solutionstructure
     71                for solution in solnames:
     72                        #looking for multiple time steps
     73                        if (len(res_struct.__dict__[solution])>num_of_timesteps):
     74                                num_of_timesteps=len(res_struct.__dict__[solution])
     75                                num_of_timesteps=int(num_of_timesteps/out_freq)+1
     76        else:
     77                num_of_timesteps=1
    6178
    62 for step in num_of_timesteps:
    63    
    64     timestep=step
    65    
    66     fid = open(strcat(path,filesep,name,filesep,name,'.vtk',int2str(timestep),'.vtk'),'w+')
    67     fid.write('# vtk DataFile Version 2.0 \n')
    68     fid.write('Data for run %s \n' % model.miscellaneous.name)
    69     fid.write('ASCII \n')
    70     fid.write('DATASET UNSTRUCTURED_GRID \n')
     79        for step in range(0,num_of_timesteps):
     80                timestep=step
     81                fid=open((filename +'/Timestep.vtk'+str(timestep)+'.vtk'),'w+')
     82                fid.write('# vtk DataFile Version 2.0 \n')
     83                fid.write('Data for run %s \n' % model.miscellaneous.name)
     84                fid.write('ASCII \n')
     85                fid.write('DATASET UNSTRUCTURED_GRID \n')
     86                fid.write('POINTS %d float\n' % num_of_points)
     87                if(dim==3):
     88                        for point in points:
     89                                fid.write('%f %f %f \n'%(point[0], point[1], point[2]))
     90                elif(dim==2):
     91                        for point in points:
     92                                fid.write('%f %f %f \n'%(point[0], point[1], point[2]))
     93                       
     94                fid.write('CELLS %d %d\n' %(num_of_elt, num_of_elt*(point_per_elt+1)))
     95               
     96                if point_per_elt==3:
     97                        for elt in range(0, num_of_elt):
     98                                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))
     99                elif point_per_elt==6:
     100                        for elt in range(0, num_of_elt):
     101                                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))
     102                else:
     103                        print 'Number of nodes per element not supported'
     104
     105                fid.write('CELL_TYPES %d\n' %num_of_elt)
     106                for elt in range(0, num_of_elt):
     107                        fid.write('%d\n' %celltype)
     108
     109                fid.write('POINT_DATA %s \n' %str(num_of_points))
    71110       
    72     fid.write(fid,'POINTS %d float\n',num_of_points)
    73     if(dim==3):
    74         s='%f %f %f \n'
    75     elif(dim==2):
    76         s='%f %f \n'
    77    
    78     P=[points zeros(num_of_points,3-dim)]
    79     fid.write(fid,s,transpose de P)
    80        
    81     fid.write('CELLS %d %d\n' % num_of_elt % num_of_elt*(point_per_elt+1))
    82     s='%d'
    83     for j=1:point_per_elt:
    84         s=horzcat(s,{' %d'})
    85        
    86     s=cell2mat(horzcat(s,{'\n'}))
    87     fid.write(fid,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]transpose)
    88    
    89     fid.write(fid,'CELL_TYPES %d\n',num_of_elt)
    90     s='%d\n'
    91     fid.write(fid,s,celltype*ones(num_of_elt,1))
    92     fid.write(fid,'POINT_DATA %s \n',num2str(num_of_points))
    93    
    94     #loop over the different solution structures
    95     if 'num_of_sols' in locals():
    96         for j=1:num_of_sols:
    97             #dealing with results on different timesteps
    98             if(numpy.size(sol_struct{j},2)>timestep):
    99                 timestep = step
    100             else:
    101                 timestep = numpy.size(sol_struct{j},2)
    102                
    103             #getting the number of fields in the solution
    104             fieldnames=fields(sol_struct{j}(timestep))
    105             num_of_fields=numpy.length(fieldnames)
    106                
    107             #check which field is a real result and print
    108             for k=1:num_of_fields:
    109                 if ((numel(sol_struct{j}(timestep).(fieldnames{k})))==num_of_points):
    110                     #paraview does not like NaN, replacing
    111                     nanval=find(isnan(sol_struct{j}(timestep).(fieldnames{k})))
    112                     sol_struct{j}(timestep).(fieldnames{k})(nanval)=-9999
    113                     #also checking for verry small value that mess up
    114                     smallval=(abs(sol_struct{j}(timestep).(fieldnames{k}))<1.0e-20)
    115                     sol_struct{j}(timestep).(fieldnames{k})(smallval)=0.0
    116                     fid.write('SCALARS %s float 1 \n' % fieldnames{k})
    117                     fid.write('LOOKUP_TABLE default\n')
    118                     s='%e\n'
    119                     fid.write(s % sol_struct{j}(timestep).(fieldnames{k}))
    120                    
    121     #loop on arguments, if something other than result is asked, do
    122     #it now
    123     for j= 1:nargin-2:
    124         res_struct=model.(varargin{j})
    125         fieldnames=fields(res_struct)
    126         num_of_fields=numpy.length(fieldnames)
    127         for k=1:num_of_fields:
    128             if ((numel(res_struct.(fieldnames{k})))==num_of_points):
    129                 #paraview does not like NaN, replacing
    130                 nanval=find(isnan(res_struct.(fieldnames{k})))
    131                 res_struct.(fieldnames{k})(nanval)=-9999
    132                 #also checking for verry small value that mess up
    133                 smallval=(abs(res_struct.(fieldnames{k}))<1.0e-20)
    134                 res_struct.(fieldnames{k})(smallval)=0.0
    135                 fid.write(fid,'SCALARS %s float 1 \n',fieldnames{k})
    136                 fid.write(fid,'LOOKUP_TABLE default\n')
    137                 s='%e\n'
    138                 fid.write(fid,s,res_struct.(fieldnames{k}))
    139     fid.close();
     111                #loop over the different solution structures
     112                if 'solnames' in locals():
     113                        for sol in solnames:
     114                                #dealing with results on different timesteps
     115                                if(len(res_struct.__dict__[sol])>timestep):
     116                                        timestep = step
     117                                else:
     118                                        timestep = len(res_struct.__dict__[sol])
     119                               
     120                                #getting the  fields in the solution
     121                                fieldnames=dict.keys(res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__)
     122                       
     123                                #check which field is a real result and print
     124                                for field in fieldnames:
     125                                        if ((numpy.size(res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__[field]))==num_of_points):
     126                                                fid.write('SCALARS %s float 1 \n' % field)
     127                                                fid.write('LOOKUP_TABLE default\n')
     128                                                for node in range(0,num_of_points):
     129                                                        #paraview does not like NaN, replacing
     130                                                        if numpy.isnan(res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__[field][node]):
     131                                                                fid.write('%e\n' % -9999.9999)
     132                                                        #also checking for verry small value that mess up
     133                                                        elif (abs(res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__[field][node])<1.0e-20):
     134                                                                fid.write('%e\n' % 0.0)
     135                                                        else:
     136                                                                fid.write('%e\n' % res_struct.__dict__[sol].__getitem__(timestep*out_freq-1).__dict__[field][node])
     137                                       
     138                #loop on arguments, if something other than result is asked, do
     139                #it now
     140                for other in args:
     141                        other_struct=model.__dict__[other]
     142                        othernames=(dict.keys(other_struct.__dict__))
     143                        for field in othernames:
     144                                if ((numpy.size(other_struct.__dict__[field]))==num_of_points):
     145                                        fid.write('SCALARS %s float 1 \n' % field)
     146                                        fid.write('LOOKUP_TABLE default\n')
     147                                        for node in range(0,num_of_points):
     148                                                #paraview does not like NaN, replacing
     149                                                if numpy.isnan(other_struct.__dict__[field][node]):
     150                                                        fid.write('%e\n' % -9999.9999)
     151                                                #also checking for verry small value that mess up
     152                                                elif (abs(other_struct.__dict__[field][node])<1.0e-20):
     153                                                        fid.write('%e\n' % 0.0)
     154                                                else:
     155                                                        fid.write('%e\n' % other_struct.__dict__[field][node])
     156        fid.close();
Note: See TracChangeset for help on using the changeset viewer.