Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 23210) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 23211) @@ -18,5 +18,5 @@ print ('New file name is {}'.format(newname)) filename=newname - + NCData=Dataset(filename, 'w', format='NETCDF4') NCData.description = 'Results for run' + md.miscellaneous.name @@ -100,10 +100,10 @@ Subgroup.__setattr__('classtype',md.__dict__[group].__class__.__name__) subfields=dict.keys(md.__dict__[group].__dict__[field].__dict__) - + for subfield in subfields: if str(subfield)!='outlog': Var=md.__dict__[group].__dict__[field].__dict__[subfield] DimDict=CreateVar(NCData,Var,subfield,Subgroup,DimDict) - + NCData.close() @@ -130,9 +130,10 @@ str:str, dict:str} - + val_dim=np.shape(val_shape)[0] + #Now define and fill up variable #treating scalar string or bool as atribute - if val_type==str or val_type==unicode or val_type==bool: + if val_type in [str,unicode,bool]: Group.__setattr__(str(field).swapcase(), str(var)) #treating list as string table @@ -147,5 +148,5 @@ if val_shape==0: ncvar= [] - else: + else: for elt in range(0,val_shape[0]): ncvar[elt] = var[elt] @@ -164,5 +165,5 @@ ncvar[elt,1]=str(dict.values(var)[elt]) #converting to str to avoid potential problems #Now dealing with numeric variables - else: + elif val_type in [float,'float64',np.float64,int,'int64']: dimensions,DimDict=GetDim(NCData,var,val_shape,DimDict,val_dim) ncvar = Group.createVariable(str(field),TypeDict[val_type],dimensions,zlib=True) @@ -175,4 +176,6 @@ except TypeError: #type does not accept nan, get vallue of the variable ncvar[:] = var + else: + print('WARNING type "{}" is unknown for "{}.{}"'.format(val_type,Group.name,field)) return DimDict Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py (revision 23210) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/enveloppeVTK.py (revision 23211) @@ -9,11 +9,11 @@ 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 + + Give only the results for nw but could be extended to geometry, mask... + + input: filename destination (string) ------------------------------------------------------------------ -model this is md +model this is md ------------------------------------------------------------------ By default only the results are exported, you can add whichever @@ -40,22 +40,38 @@ os.mkdir(filename) - IsEnveloppe=np.where(model.mesh.vertexonbase | model.mesh.vertexonsurface) - #get the element related variables + # {{{get the element related variables if 'z' in dict.keys(model.mesh.__dict__): - points=np.column_stack((model.mesh.x,model.mesh.y,model.mesh.z)) - num_of_elt=np.size(np.isnan(model.mesh.lowerelements))+np.size(np.isnan(model.mesh.upperelements)) - low_elt_num=np.size(np.isnan(model.mesh.lowerelements)) - top_elt_num=np.size(np.isnan(model.mesh.upperelements)) + is_enveloppe=np.logical_or(model.mesh.vertexonbase,model.mesh.vertexonsurface) + enveloppe_index=np.where(is_enveloppe)[0] + convert_index=np.nan*np.ones(np.shape(model.mesh.x)) + convert_index=np.asarray([[i,np.where(enveloppe_index==i)[0][0]] for i,val in enumerate(convert_index) if any(enveloppe_index==i)]) + points=np.column_stack((model.mesh.x[enveloppe_index], + model.mesh.y[enveloppe_index], + model.mesh.z[enveloppe_index])) + low_elt_num=np.size(np.where(np.isnan(model.mesh.lowerelements))) + top_elt_num=np.size(np.where(np.isnan(model.mesh.upperelements))) + num_of_elt=low_elt_num+top_elt_num + connect=model.mesh.elements[np.where(is_enveloppe[model.mesh.elements-1])].reshape(int(num_of_elt),3)-1 + for elt in range(0, num_of_elt): + connect[elt,0]=convert_index[np.where(convert_index==connect[elt,0])[0],1][0] + connect[elt,1]=convert_index[np.where(convert_index==connect[elt,1])[0],1][0] + connect[elt,2]=convert_index[np.where(convert_index==connect[elt,2])[0],1][0] + else: - points=np.column_stack((model.mesh.x,model.mesh.y,np.zeros(np.shape(model.mesh.x)))) + points=np.column_stack((model.mesh.x, + model.mesh.y, + np.zeros(np.shape(model.mesh.x)))) num_of_elt=np.shape(model.mesh.elements)[0] - - num_of_points=np.size(points)[0] - dim=np.size(points)[1] - point_per_elt=np.shape(model.mesh.elements)[1] - + connect=model.mesh.elements-1 + enveloppe_index=np.arange(0,np.size(model.mesh.x)) + + every_nodes=np.size(model.mesh.x) + num_of_points=np.size(enveloppe_index) + dim=3 + point_per_elt=3 celltype=5 #triangles - - #this is the result structure + + # }}} + # {{{this is the result structure res_struct=model.results if (len(res_struct.__dict__)>0): @@ -64,5 +80,5 @@ num_of_sols=len(solnames) num_of_timesteps=1 - #%building solutionstructure + #%building solutionstructure for solution in solnames: #looking for multiple time steps @@ -72,5 +88,6 @@ else: num_of_timesteps=1 - + # }}} + # {{{write header and mesh for step in range(0,num_of_timesteps): timestep=step @@ -83,17 +100,12 @@ 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 exist('low_elt_num') - # triaconnect=zeros(num_of_elt,3); - # triaconnect(1:low_elt_num,:)=model.mesh.elements(find(isnan(model.mesh.lowerelements)),1:3); - # 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))); - # triaconnect(1+low_elt_num:num_of_elt,:)=model.mesh.elements(find(isnan(model.mesh.upperelements)),4:6)+upshift; - # fprintf(fid,s,[(3)*ones(num_of_elt,1) triaconnect-1]'); 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)) - + fid.write('3 %d %d %d\n' %(connect[elt,0], + connect[elt,1], + connect[elt,2])) + fid.write('CELL_TYPES %d\n' %num_of_elt) for elt in range(0, num_of_elt): @@ -101,6 +113,6 @@ fid.write('POINT_DATA %s \n' %str(num_of_points)) - - #loop over the different solution structures + # }}} + # {{{loop over the different solution structures if 'solnames' in locals(): for sol in solnames: @@ -110,5 +122,5 @@ else: timestep = np.size(res_struct.__dict__[sol]) - + #getting the fields in the solution if(np.size(res_struct.__dict__[sol])>1): @@ -123,35 +135,36 @@ fieldstruct=res_struct.__dict__[sol].__dict__[field] - if ((np.size(fieldstruct))==num_of_points): + if ((np.size(fieldstruct))==every_nodes): 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 np.isnan(fieldstruct[node]): + if np.isnan(fieldstruct[enveloppe_index[node]]): fid.write('%e\n' % -9999.9999) #also checking for verry small value that mess up - elif (abs(fieldstruct[node])<1.0e-20): + elif (abs(fieldstruct[enveloppe_index[node]])<1.0e-20): fid.write('%e\n' % 0.0) else: - fid.write('%e\n' % fieldstruct[node]) - - #loop on arguments, if something other than result is asked, do - #it now + fid.write('%e\n' % fieldstruct[enveloppe_index[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: -# fprintf(fid,s,res_struct.(fieldnames{k})(IsEnveloppe)); - if ((np.size(other_struct.__dict__[field]))==num_of_points): + if ((np.size(other_struct.__dict__[field]))==every_nodes): 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 np.isnan(other_struct.__dict__[field][node]): + if np.isnan(other_struct.__dict__[field][enveloppe_index[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): + elif (abs(other_struct.__dict__[field][enveloppe_index[node]])<1.0e-20): fid.write('%e\n' % 0.0) else: - fid.write('%e\n' % other_struct.__dict__[field][node]) + fid.write('%e\n' % other_struct.__dict__[field][enveloppe_index[node]]) + + # }}} fid.close(); Index: /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 23210) +++ /issm/trunk-jpl/src/m/contrib/defleurian/paraview/exportVTK.py (revision 23211) @@ -9,11 +9,11 @@ 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 + + Give only the results for nw but could be extended to geometry, mask... + + input: filename destination (string) ------------------------------------------------------------------ -model this is md +model this is md ------------------------------------------------------------------ By default only the results are exported, you can add whichever @@ -67,5 +67,5 @@ num_of_sols=len(solnames) num_of_timesteps=1 - #%building solutionstructure + #%building solutionstructure for solution in solnames: #looking for multiple time steps @@ -91,7 +91,7 @@ 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): @@ -117,5 +117,5 @@ else: timestep = np.size(res_struct.__dict__[sol]) - + #getting the fields in the solution if(np.size(res_struct.__dict__[sol])>1): @@ -187,11 +187,14 @@ # reloaded variable are generally of dim 2 elif np.shape(other_struct.__dict__[field])[0]==num_of_points: + # we want only vector + if np.shape(other_struct.__dict__[field])[1]==1: 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 + print other_struct.__dict__[field][node] if np.isnan(other_struct.__dict__[field][node]): fid.write('%e\n' % -9999.9999) - #also checking for verry small value that mess up + #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)