Index: /issm/trunk-jpl/src/m/contrib/paraview/enveloppeVTK.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/paraview/enveloppeVTK.m (revision 18513) +++ /issm/trunk-jpl/src/m/contrib/paraview/enveloppeVTK.m (revision 18513) @@ -0,0 +1,155 @@ +function enveloppeVTK(filename,model,varargin) +% vtk export +% function enveloppeVTK(filename,model) +% creates a directory with the vtk files for displays in paraview +% only export the enveloppe result (surface and base) on trias +% +% 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: + +[path,name,ext]=fileparts(filename); +separator=filesep; +mkdir(filename); +IsEnveloppe=find(model.mesh.vertexonbase | model.mesh.vertexonsurface); + +%get the element related variables +if dimension(model.mesh)==2, + points=[model.mesh.x model.mesh.y zeros(model.mesh.numberofvertices,1)]; + [num_of_elt]=size(model.mesh.elements,1); + [point_per_elt]=size(model.mesh.elements,2); +else + points=[model.mesh.x(IsEnveloppe) model.mesh.y(IsEnveloppe) model.mesh.z(IsEnveloppe)]; + [num_of_elt]=size(find(isnan(model.mesh.lowerelements)),1)+size(find(isnan(model.mesh.upperelements)),1); + [low_elt_num]=size(find(isnan(model.mesh.lowerelements)),1); + [top_elt_num]=size(find(isnan(model.mesh.upperelements)),1); +end + +celltype=5; %triangles +[num_of_points,dim]=size(points); +tot_points=model.mesh.numberofvertices; + +%this is the result structure +res_struct=model.results; +%checking for results +if (length(fields(res_struct))>0); + %Getting all the solutions of the model + solnames=fields(res_struct); + num_of_sols=length(solnames); + num_of_timesteps=1; + %building solution structure + for i=1:num_of_sols + sol_struct{i}=res_struct.(solnames{i}); + %looking for multiple time steps + if(size(sol_struct{i},2)>num_of_timesteps); + num_of_timesteps=size(sol_struct{i},2); + end + end +else + num_of_timesteps=1; +end +for step=1:num_of_timesteps; + + timestep=step; + + fid = fopen(strcat(path,filesep,name,filesep,name,'.vtk',int2str(timestep),'.vtk'),'w+'); + fprintf(fid,'# vtk DataFile Version 2.0 \n'); + fprintf(fid,'Data for run %s \n',model.miscellaneous.name); + fprintf(fid,'ASCII \n'); + fprintf(fid,'DATASET UNSTRUCTURED_GRID \n'); + + fprintf(fid,'POINTS %d float\n',num_of_points); + if(dim==3); + s='%f %f %f \n'; + elseif(dim==2); + s='%f %f \n'; + end + P=[points zeros(num_of_points,3-dim)]; + fprintf(fid,s,P'); + + fprintf(fid,'CELLS %d %d\n',num_of_elt,num_of_elt*(3+1)); + s='%d'; + for j=1:3 + s=horzcat(s,{' %d'}); + end + s=cell2mat(horzcat(s,{'\n'})); + + %build the connection matrix for the top and bottom elements + 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]'); + else + fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]'); + end + + fprintf(fid,'CELL_TYPES %d\n',num_of_elt); + s='%d\n'; + fprintf(fid,s,celltype*ones(num_of_elt,1)); + fprintf(fid,'POINT_DATA %s \n',num2str(num_of_points)); + + %loop over the different solution structures + if (exist('num_of_sols')); + for j=1:num_of_sols + %dealing with results on different timesteps + if(size(sol_struct{j},2)>timestep); + timestep = step; + else + timestep = size(sol_struct{j},2); + end + + %getting the number of fields in the solution + fieldnames=fields(sol_struct{j}(timestep)); + num_of_fields=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})))==tot_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; + fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k}); + fprintf(fid,'LOOKUP_TABLE default\n'); + s='%e\n'; + fprintf(fid,s,sol_struct{j}(timestep).(fieldnames{k})(IsEnveloppe)); + end + end + end + end + %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=length(fieldnames); + for k=1:num_of_fields + if ((numel(res_struct.(fieldnames{k})))==tot_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; + fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k}); + fprintf(fid,'LOOKUP_TABLE default\n'); + s='%e\n'; + fprintf(fid,s,res_struct.(fieldnames{k})(IsEnveloppe)); + end + end + end + fclose(fid); +end Index: /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.m (revision 18512) +++ /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.m (revision 18513) @@ -25,11 +25,14 @@ if dimension(model.mesh)==2, points=[model.mesh.x model.mesh.y zeros(model.mesh.numberofvertices,1)]; + [num_of_points,dim]=size(points); + [num_of_elt]=size(model.mesh.elements,1); + [point_per_elt]=size(model.mesh.elements,2); else points=[model.mesh.x model.mesh.y model.mesh.z]; + [num_of_points,dim]=size(points); + [num_of_elt]=size(model.mesh.elements,1); + [point_per_elt]=size(model.mesh.elements,2); end -[num_of_points,dim]=size(points); -[num_of_elt]=size(model.mesh.elements,1); -[point_per_elt]=size(model.mesh.elements,2); %Select the type of element function of the number of nodes per elements @@ -86,5 +89,5 @@ end s=cell2mat(horzcat(s,{'\n'})); - fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]'); + fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]'); fprintf(fid,'CELL_TYPES %d\n',num_of_elt); Index: /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.py (revision 18513) +++ /issm/trunk-jpl/src/m/contrib/paraview/exportVTK.py (revision 18513) @@ -0,0 +1,131 @@ +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: + ''' + +[path,name,ext]=fileparts(filename) +separator=filesep +mkdir(filename) + +#get the element related variables +points=[model.mesh.x model.mesh.y model.mesh.z] +[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) + +#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 +#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); + +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') + + 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 + 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();