![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 23211
- Timestamp:
- 09/04/18 06:16:04 (7 years ago)
- Location:
- issm/trunk-jpl/src/m/contrib/defleurian
- Files:
-
- 3 edited
-
netCDF/export_netCDF.py (modified) (6 diffs)
-
paraview/enveloppeVTK.py (modified) (8 diffs)
-
paraview/exportVTK.py (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py
r21644 r23211 18 18 print ('New file name is {}'.format(newname)) 19 19 filename=newname 20 20 21 21 NCData=Dataset(filename, 'w', format='NETCDF4') 22 22 NCData.description = 'Results for run' + md.miscellaneous.name … … 100 100 Subgroup.__setattr__('classtype',md.__dict__[group].__class__.__name__) 101 101 subfields=dict.keys(md.__dict__[group].__dict__[field].__dict__) 102 102 103 103 for subfield in subfields: 104 104 if str(subfield)!='outlog': 105 105 Var=md.__dict__[group].__dict__[field].__dict__[subfield] 106 106 DimDict=CreateVar(NCData,Var,subfield,Subgroup,DimDict) 107 107 108 108 NCData.close() 109 109 … … 130 130 str:str, 131 131 dict:str} 132 132 133 133 val_dim=np.shape(val_shape)[0] 134 134 135 #Now define and fill up variable 135 136 #treating scalar string or bool as atribute 136 if val_type ==str or val_type==unicode or val_type==bool:137 if val_type in [str,unicode,bool]: 137 138 Group.__setattr__(str(field).swapcase(), str(var)) 138 139 #treating list as string table … … 147 148 if val_shape==0: 148 149 ncvar= [] 149 else: 150 else: 150 151 for elt in range(0,val_shape[0]): 151 152 ncvar[elt] = var[elt] … … 164 165 ncvar[elt,1]=str(dict.values(var)[elt]) #converting to str to avoid potential problems 165 166 #Now dealing with numeric variables 166 el se:167 elif val_type in [float,'float64',np.float64,int,'int64']: 167 168 dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim) 168 169 ncvar = Group.createVariable(str(field),TypeDict[val_type],dimensions,zlib=True) … … 175 176 except TypeError: #type does not accept nan, get vallue of the variable 176 177 ncvar[:] = var 178 else: 179 print('WARNING type "{}" is unknown for "{}.{}"'.format(val_type,Group.name,field)) 177 180 return DimDict 178 181 -
issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py
r21303 r23211 9 9 creates a directory with the vtk files for displays in paraview 10 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 11 12 Give only the results for nw but could be extended to geometry, mask... 13 14 input: filename destination 15 15 (string) 16 16 ------------------------------------------------------------------ 17 model this is md 17 model this is md 18 18 ------------------------------------------------------------------ 19 19 By default only the results are exported, you can add whichever … … 40 40 os.mkdir(filename) 41 41 42 IsEnveloppe=np.where(model.mesh.vertexonbase | model.mesh.vertexonsurface) 43 #get the element related variables 42 # {{{get the element related variables 44 43 if 'z' in dict.keys(model.mesh.__dict__): 45 points=np.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) 46 num_of_elt=np.size(np.isnan(model.mesh.lowerelements))+np.size(np.isnan(model.mesh.upperelements)) 47 low_elt_num=np.size(np.isnan(model.mesh.lowerelements)) 48 top_elt_num=np.size(np.isnan(model.mesh.upperelements)) 44 is_enveloppe=np.logical_or(model.mesh.vertexonbase,model.mesh.vertexonsurface) 45 enveloppe_index=np.where(is_enveloppe)[0] 46 convert_index=np.nan*np.ones(np.shape(model.mesh.x)) 47 convert_index=np.asarray([[i,np.where(enveloppe_index==i)[0][0]] for i,val in enumerate(convert_index) if any(enveloppe_index==i)]) 48 points=np.column_stack((model.mesh.x[enveloppe_index], 49 model.mesh.y[enveloppe_index], 50 model.mesh.z[enveloppe_index])) 51 low_elt_num=np.size(np.where(np.isnan(model.mesh.lowerelements))) 52 top_elt_num=np.size(np.where(np.isnan(model.mesh.upperelements))) 53 num_of_elt=low_elt_num+top_elt_num 54 connect=model.mesh.elements[np.where(is_enveloppe[model.mesh.elements-1])].reshape(int(num_of_elt),3)-1 55 for elt in range(0, num_of_elt): 56 connect[elt,0]=convert_index[np.where(convert_index==connect[elt,0])[0],1][0] 57 connect[elt,1]=convert_index[np.where(convert_index==connect[elt,1])[0],1][0] 58 connect[elt,2]=convert_index[np.where(convert_index==connect[elt,2])[0],1][0] 59 49 60 else: 50 points=np.column_stack((model.mesh.x,model.mesh.y,np.zeros(np.shape(model.mesh.x)))) 61 points=np.column_stack((model.mesh.x, 62 model.mesh.y, 63 np.zeros(np.shape(model.mesh.x)))) 51 64 num_of_elt=np.shape(model.mesh.elements)[0] 52 53 num_of_points=np.size(points)[0] 54 dim=np.size(points)[1] 55 point_per_elt=np.shape(model.mesh.elements)[1] 56 65 connect=model.mesh.elements-1 66 enveloppe_index=np.arange(0,np.size(model.mesh.x)) 67 68 every_nodes=np.size(model.mesh.x) 69 num_of_points=np.size(enveloppe_index) 70 dim=3 71 point_per_elt=3 57 72 celltype=5 #triangles 58 59 #this is the result structure 73 74 # }}} 75 # {{{this is the result structure 60 76 res_struct=model.results 61 77 if (len(res_struct.__dict__)>0): … … 64 80 num_of_sols=len(solnames) 65 81 num_of_timesteps=1 66 #%building solutionstructure 82 #%building solutionstructure 67 83 for solution in solnames: 68 84 #looking for multiple time steps … … 72 88 else: 73 89 num_of_timesteps=1 74 90 # }}} 91 # {{{write header and mesh 75 92 for step in range(0,num_of_timesteps): 76 93 timestep=step … … 83 100 for point in points: 84 101 fid.write('%f %f %f \n'%(point[0], point[1], point[2])) 85 102 86 103 fid.write('CELLS %d %d\n' %(num_of_elt, num_of_elt*(point_per_elt+1))) 87 104 88 # if exist('low_elt_num')89 # triaconnect=zeros(num_of_elt,3);90 # triaconnect(1:low_elt_num,:)=model.mesh.elements(find(isnan(model.mesh.lowerelements)),1:3);91 # upshift=-min(min(model.mesh.elements(find(isnan(model.mesh.upperelements)),4:6)))+1+max(max(model.mesh.elements(find(isnan(model.mesh.lowerelements)),1:3)));92 # triaconnect(1+low_elt_num:num_of_elt,:)=model.mesh.elements(find(isnan(model.mesh.upperelements)),4:6)+upshift;93 # fprintf(fid,s,[(3)*ones(num_of_elt,1) triaconnect-1]');94 105 for elt in range(0, num_of_elt): 95 96 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)) 97 106 fid.write('3 %d %d %d\n' %(connect[elt,0], 107 connect[elt,1], 108 connect[elt,2])) 109 98 110 fid.write('CELL_TYPES %d\n' %num_of_elt) 99 111 for elt in range(0, num_of_elt): … … 101 113 102 114 fid.write('POINT_DATA %s \n' %str(num_of_points)) 103 104 # loop over the different solution structures115 # }}} 116 # {{{loop over the different solution structures 105 117 if 'solnames' in locals(): 106 118 for sol in solnames: … … 110 122 else: 111 123 timestep = np.size(res_struct.__dict__[sol]) 112 124 113 125 #getting the fields in the solution 114 126 if(np.size(res_struct.__dict__[sol])>1): … … 123 135 fieldstruct=res_struct.__dict__[sol].__dict__[field] 124 136 125 if ((np.size(fieldstruct))== num_of_points):137 if ((np.size(fieldstruct))==every_nodes): 126 138 fid.write('SCALARS %s float 1 \n' % field) 127 139 fid.write('LOOKUP_TABLE default\n') 128 140 for node in range(0,num_of_points): 129 141 #paraview does not like NaN, replacing 130 if np.isnan(fieldstruct[ node]):142 if np.isnan(fieldstruct[enveloppe_index[node]]): 131 143 fid.write('%e\n' % -9999.9999) 132 144 #also checking for verry small value that mess up 133 elif (abs(fieldstruct[ node])<1.0e-20):145 elif (abs(fieldstruct[enveloppe_index[node]])<1.0e-20): 134 146 fid.write('%e\n' % 0.0) 135 147 else: 136 fid.write('%e\n' % fieldstruct[ node])137 138 # loop on arguments, if something other than result is asked, do139 #it now 148 fid.write('%e\n' % fieldstruct[enveloppe_index[node]]) 149 # }}} 150 # {{{loop on arguments, if something other than result is asked, do it now 151 140 152 for other in args: 141 153 other_struct=model.__dict__[other] 142 154 othernames=(dict.keys(other_struct.__dict__)) 143 155 for field in othernames: 144 # fprintf(fid,s,res_struct.(fieldnames{k})(IsEnveloppe)); 145 if ((np.size(other_struct.__dict__[field]))==num_of_points): 156 if ((np.size(other_struct.__dict__[field]))==every_nodes): 146 157 fid.write('SCALARS %s float 1 \n' % field) 147 158 fid.write('LOOKUP_TABLE default\n') 148 159 for node in range(0,num_of_points): 149 160 #paraview does not like NaN, replacing 150 if np.isnan(other_struct.__dict__[field][ node]):161 if np.isnan(other_struct.__dict__[field][enveloppe_index[node]]): 151 162 fid.write('%e\n' % -9999.9999) 152 163 #also checking for verry small value that mess up 153 elif (abs(other_struct.__dict__[field][ node])<1.0e-20):164 elif (abs(other_struct.__dict__[field][enveloppe_index[node]])<1.0e-20): 154 165 fid.write('%e\n' % 0.0) 155 166 else: 156 fid.write('%e\n' % other_struct.__dict__[field][node]) 167 fid.write('%e\n' % other_struct.__dict__[field][enveloppe_index[node]]) 168 169 # }}} 157 170 fid.close(); -
issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py
r21568 r23211 9 9 creates a directory with the vtk files for displays in paraview 10 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 11 12 Give only the results for nw but could be extended to geometry, mask... 13 14 input: filename destination 15 15 (string) 16 16 ------------------------------------------------------------------ 17 model this is md 17 model this is md 18 18 ------------------------------------------------------------------ 19 19 By default only the results are exported, you can add whichever … … 67 67 num_of_sols=len(solnames) 68 68 num_of_timesteps=1 69 #%building solutionstructure 69 #%building solutionstructure 70 70 for solution in solnames: 71 71 #looking for multiple time steps … … 91 91 for point in points: 92 92 fid.write('%f %f %f \n'%(point[0], point[1], point[2])) 93 93 94 94 fid.write('CELLS %d %d\n' %(num_of_elt, num_of_elt*(point_per_elt+1))) 95 95 96 96 if point_per_elt==3: 97 97 for elt in range(0, num_of_elt): … … 117 117 else: 118 118 timestep = np.size(res_struct.__dict__[sol]) 119 119 120 120 #getting the fields in the solution 121 121 if(np.size(res_struct.__dict__[sol])>1): … … 187 187 # reloaded variable are generally of dim 2 188 188 elif np.shape(other_struct.__dict__[field])[0]==num_of_points: 189 # we want only vector 190 if np.shape(other_struct.__dict__[field])[1]==1: 189 191 fid.write('SCALARS %s float 1 \n' % field) 190 192 fid.write('LOOKUP_TABLE default\n') 191 193 for node in range(0,num_of_points): 192 194 #paraview does not like NaN, replacing 195 print other_struct.__dict__[field][node] 193 196 if np.isnan(other_struct.__dict__[field][node]): 194 197 fid.write('%e\n' % -9999.9999) 195 #also checking for verry small value that mess up198 #also checking for verry small value that mess up 196 199 elif (abs(other_struct.__dict__[field][node])<1.0e-20): 197 200 fid.write('%e\n' % 0.0)
Note:
See TracChangeset
for help on using the changeset viewer.
