Changeset 19197

Timestamp:

03/16/15 10:23:06 (10 years ago)

Author:

bdef

Message:

NEW: adding support for forcings

Location:

issm/trunk-jpl/src/m/contrib/paraview

Files:

• enveloppeVTK.m (modified) (5 diffs)
• exportVTK.m (modified) (2 diffs)

issm/trunk-jpl/src/m/contrib/paraview/enveloppeVTK.m

@@ -51 +51 @@
                 if(size(sol_struct{i},2)>num_of_timesteps);
                         num_of_timesteps=size(sol_struct{i},2);
+                        outstep=model.timestepping.time_step*model.settings.output_frequency
+    end
   end
- end
 else
         num_of_timesteps=1;
@@ -59 +60 @@
         
         timestep=step;
-
         fid = fopen(strcat(path,filesep,name,filesep,'timestep.vtk',int2str(timestep),'.vtk'),'w+');
         fprintf(fid,'# vtk DataFile Version 2.0 \n');
@@ -85 +85 @@
         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,:)=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]');
@@ -109 +107 @@
                                 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);
-                        
+                        resfields=fields(sol_struct{j}(timestep));
+                        num_of_fields=length(resfields);
                         %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);
+                                if ((numel(sol_struct{j}(timestep).(resfields{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;
+                                        nanval=find(isnan(sol_struct{j}(timestep).(resfields{k})));
+                                        sol_struct{j}(timestep).(resfields{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});
+                                        smallval=(abs(sol_struct{j}(timestep).(resfields{k}))<1.0e-20);
+                                        sol_struct{j}(timestep).(resfields{k})(smallval)=0.0;
+                                        fprintf(fid,'SCALARS %s float 1 \n',resfields{k});
                                         fprintf(fid,'LOOKUP_TABLE default\n');
                                         s='%e\n';
-                                        fprintf(fid,s,sol_struct{j}(timestep).(fieldnames{k})(IsEnveloppe));
+                                        fprintf(fid,s,sol_struct{j}(timestep).(resfields{k})(IsEnveloppe));
                     end         
             end 
@@ -149 +145 @@
                                 s='%e\n';
                                 fprintf(fid,s,res_struct.(fieldnames{k})(IsEnveloppe));
-            end         
+                                %check for forcings     
+                        elseif (size(res_struct.(fieldnames{k}),1)==tot_points+1);
+                                %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;
+                                if (size(res_struct.(fieldnames{k}),2)==num_of_timesteps),
+                                        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,timestep));
+                                else,
+                                        %forcing and results not on the same timestep,need some treatment
+                                        fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k});
+                                        fprintf(fid,'LOOKUP_TABLE default\n');
+                                        index=1
+                                        currenttime=((timestep-1)*outstep)+model.timestepping.start_time+model.timestepping.time_step
+                                        while (res_struct.(fieldnames{k})(end,index)<=currenttime);
+                                                index=index+1
+                      end
+                                        uptime=res_struct.(fieldnames{k})(end,index);
+                                        uplim=res_struct.(fieldnames{k})(IsEnveloppe,index);
+                                        uptime
+                                        while (res_struct.(fieldnames{k})(end,index)>=currenttime);
+                                                index=index-1
+                      end
+                                        lowtime=res_struct.(fieldnames{k})(end,index);
+                                        lowlim=res_struct.(fieldnames{k})(IsEnveloppe,index);
+                                        lowtime
+                                        interp=lowlim+(uplim-lowlim)*((currenttime-lowtime)/(uptime-lowtime))
+                                        s='%e\n';
+                                        fprintf(fid,s,interp);
+                                end     
+                  end           
                 end 
         end

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

@@ -59 +59 @@
                 if(size(sol_struct{i},2)>num_of_timesteps);
                         num_of_timesteps=size(sol_struct{i},2);
-    end
+          end
+                outstep=model.timestepping.time_step*model.settings.output_frequency;
   end
 else
@@ -145 +146 @@
                                 s='%e\n';
                                 fprintf(fid,s,res_struct.(fieldnames{k}));
-            end         
+                                %check for forcings     
+                        elseif (size(res_struct.(fieldnames{k}),1)==num_of_points+1);
+                                %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;
+                                if (size(res_struct.(fieldnames{k}),2)==num_of_timesteps),
+                                        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})(1:end-1,timestep));
+                                else,
+                                        %forcing and results not on the same timestep,need some treatment
+                                        fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k});
+                                        fprintf(fid,'LOOKUP_TABLE default\n');
+                                        index=1;
+                                        currenttime=((timestep-1)*outstep)+model.timestepping.start_time;
+                                        while (res_struct.(fieldnames{k})(end,index)<=currenttime);
+                                                if index==size(res_struct.(fieldnames{k}),2)
+                                                        break
+                                                end     
+                                                index=index+1;
+                      end
+                                        uptime=res_struct.(fieldnames{k})(end,index);
+                                        uplim=res_struct.(fieldnames{k})(1:end-1,index);
+                                        while (res_struct.(fieldnames{k})(end,index)>=currenttime);
+                                                if index==1
+                                                        break
+                              end
+                                                index=index-1;
+                      end
+                                        lowtime=res_struct.(fieldnames{k})(end,index);
+                                        lowlim=res_struct.(fieldnames{k})(1:end-1,index);
+                                        if uptime==currenttime,
+                                                interp=uplim;
+                                        elseif lowtime==currenttime,
+                                                interp=lowlim;
+                                        else
+                                                interp=lowlim+(uplim-lowlim)*((currenttime-lowtime)/(uptime-lowtime));
+                                        end
+                                        s='%e\n';
+                                        fprintf(fid,s,interp);
+                                end
+                  end           
                 end 
         end