diff --git a/src/matlab/add_patch2slice.m b/src/matlab/add_patch2slice.m index 0b28fbe0..d4233d95 100644 --- a/src/matlab/add_patch2slice.m +++ b/src/matlab/add_patch2slice.m @@ -1,6 +1,6 @@ function add_patch2slice(sdir,sval,snum,xc,yc,zc, xe,ye,ze, q,level,... contourlevels,mappedgrid,manifold,maskflag,grid_number,blockno,... - colormapping); + colormapping) % Internal matlab routine for Clawpack graphics. @@ -13,7 +13,7 @@ function add_patch2slice(sdir,sval,snum,xc,yc,zc, xe,ye,ze, q,level,... % Mask any patches that might be underneath the new patch we want to add. if (maskflag == 1) mask_patches(slice,sdir,level,xe(1),xe(end),ye(1),ye(end),ze(1),ze(end)); -end; +end % Get the new patch new_patch = create_patch(xc,yc,zc,xe,ye,ze,q,sdir,sval,contourlevels,... diff --git a/src/matlab/create_clines.m b/src/matlab/create_clines.m index 4d156a87..9c1e714e 100644 --- a/src/matlab/create_clines.m +++ b/src/matlab/create_clines.m @@ -9,52 +9,52 @@ global creatingclines if (isempty(c)) - % No contour lines on this patch - h = []; - return; -end; + % No contour lines on this patch + h = []; + return; +end st_idx = 1; line_cnt = 0; while (1) - cval = c(1,st_idx); - next_length = c(2,st_idx); - line_cnt = line_cnt + 1; - x_like = zeros(1,next_length) + sval; - y_like = c(1,st_idx+1:st_idx+next_length); - z_like = c(2,st_idx+1:st_idx+next_length); - - [xdata,ydata,zdata] = get_xyz(x_like,y_like,z_like,sdir); - - udata.cartCoords = [xdata', ydata', zdata']; - if (mappedgrid == 1 | manifold == 1) - if (mappedgrid == 1) - if (nargin('mapc2p') == 2) - [xdata,ydata] = mapc2p(xdata,ydata); - else - [xdata, ydata, zdata] = mapc2p(xdata,ydata,zdata); - end; - end; - if (manifold == 1) - creatingclines = 1; % flag for mapc2m to indicate contour lines - % are being generated: may want to shift off - % manifold slightly to avoid hidden line removal - % causing apparent gaps in contours. - set_blocknumber(blockno); - [xdata, ydata, zdata] = mapc2m(xdata,ydata); - creatingclines = 0; - end; - end; - udata.phys_vertices = [xdata', ydata', zdata']; - - h(line_cnt) = line('XData',xdata,'YData',ydata,'ZData',zdata,'Color','k'); - set(h(line_cnt),'Tag','on'); - udata.cval = cval; - set(h(line_cnt),'UserData',udata); - st_idx = st_idx + next_length + 1; - if (st_idx > length(c)) - return; - end; -end; + cval = c(1,st_idx); + next_length = c(2,st_idx); + line_cnt = line_cnt + 1; + x_like = zeros(1,next_length) + sval; + y_like = c(1,st_idx+1:st_idx+next_length); + z_like = c(2,st_idx+1:st_idx+next_length); + + [xdata,ydata,zdata] = get_xyz(x_like,y_like,z_like,sdir); + + udata.cartCoords = [xdata', ydata', zdata']; + if (mappedgrid == 1 || manifold == 1) + if (mappedgrid == 1) + if (nargin('mapc2p') == 2) + [xdata,ydata] = mapc2p(xdata,ydata); + else + [xdata, ydata, zdata] = mapc2p(xdata,ydata,zdata); + end + end + if (manifold == 1) + creatingclines = 1; % flag for mapc2m to indicate contour lines + % are being generated: may want to shift off + % manifold slightly to avoid hidden line removal + % causing apparent gaps in contours. + set_blocknumber(blockno); + [xdata, ydata, zdata] = mapc2m(xdata,ydata); + creatingclines = 0; + end + end + udata.phys_vertices = [xdata', ydata', zdata']; + + h(line_cnt) = line('XData',xdata,'YData',ydata,'ZData',zdata,'Color','k'); + set(h(line_cnt),'Tag','on'); + udata.cval = cval; + set(h(line_cnt),'UserData',udata); + st_idx = st_idx + next_length + 1; + if (st_idx > length(c)) + break; + end +end set_clines(h); diff --git a/src/matlab/create_patch.m b/src/matlab/create_patch.m index 66c8fc26..f336b803 100644 --- a/src/matlab/create_patch.m +++ b/src/matlab/create_patch.m @@ -4,8 +4,8 @@ % Internal matlab routine for Clawpack graphics. -[xe_like, ye_like, ze_like] = get_xyzlike(xe,ye,ze,sdir); -[xc_like, yc_like, zc_like] = get_xyzlike(xc,yc,zc,sdir); +[~, ye_like, ze_like] = get_xyzlike(xe,ye,ze,sdir); +[~, yc_like, zc_like] = get_xyzlike(xc,yc,zc,sdir); % -----------x--------------------------------------- % Create patch with q data @@ -49,43 +49,51 @@ [mv_names{1:3}] = get_xyzlike('mx','my','mz',sdir); [dv_names{1:3}] = get_xyzlike('dx','dy','dz',sdir); -for i = 1:3, - userdata = setfield(userdata,(vc_names{i}),vc{i}); - userdata = setfield(userdata,(ve_names{i}), ve{i}); - userdata = setfield(userdata,(mv_names{i}),length(vc{i})); +for i = 1:3 + % userdata = setfield(userdata,(vc_names{i}),vc{i}); + % userdata = setfield(userdata,(ve_names{i}), ve{i}); + % userdata = setfield(userdata,(mv_names{i}),length(vc{i})); + userdata.(vc_names{i}) = vc{i}; + userdata.(ve_names{i}) = ve{i}; + userdata.(mv_names{i}) = numel(vc{i}); % If we are doing a manifold, then dz==1. dv = ve{i}(2) - ve{i}(1); - userdata = setfield(userdata,(dv_names{i}),(dv == 0) + dv*(dv ~= 0)); -end; + % userdata = setfield(userdata,(dv_names{i}),(dv == 0) + dv*(dv ~= 0)); + userdata.(dv_names{i}) = (dv == 0) + dv*(dv ~= 0); +end [vmin_names{1:3}] = get_xyzlike('xmin','ymin','zmin',sdir); [vmax_names{1:3}] = get_xyzlike('xmax','ymax','zmax',sdir); -userdata = setfield(userdata,vmin_names{1},sval); -userdata = setfield(userdata,vmax_names{1},sval); -for i = 2:3, - userdata = setfield(userdata,(vmin_names{i}), ve{i}(1)); - userdata = setfield(userdata,(vmax_names{i}), ve{i}(end)); -end; +% userdata = setfield(userdata,vmin_names{1},sval); +% userdata = setfield(userdata,vmax_names{1},sval); +userdata.(vmin_names{1}) = sval; +userdata.(vmax_names{1}) = sval; +for i = 2:3 + %userdata = setfield(userdata,(vmin_names{i}), ve{i}(1)); + %userdata = setfield(userdata,(vmax_names{i}), ve{i}(end)); + userdata.(vmin_names{i}) = ve{i}(1); + userdata.(vmax_names{i}) = ve{i}(end); +end % Store Cartesian coodinates for use in mask_patches, and convert patch % vertices to physical coordinates. v = get(p,'Vertices'); userdata.cartCoords = v; -if (mappedgrid == 1 | manifold == 1) +if (mappedgrid == 1 || manifold == 1) if (mappedgrid == 1) if (nargin('mapc2p') == 2) [v(:,1),v(:,2)] = mapc2p(v(:,1),v(:,2)); else [v(:,1),v(:,2),v(:,3)] = mapc2p(v(:,1),v(:,2),v(:,3)); - end; - end; + end + end if (manifold == 1) [v(:,1),v(:,2), v(:,3)] = mapc2m(v(:,1),v(:,2)); - end; + end set(p,'Vertices',v); -end; +end userdata.phys_vertices = v; % ------------------------------------------------------- @@ -97,7 +105,7 @@ if (~isempty(contourlevels)) c = contourc(yc_like,zc_like,qcm2,contourlevels); userdata.contourLines = create_clines(c,sval,sdir,mappedgrid,manifold,blockno); -end; +end % Mesh data for showing coarsened mesh later... % userdata.mesh = create_mesh(sdir,sval,xe,ye,ze,mappedgrid,manifold); diff --git a/src/matlab/create_sliceintersections.m b/src/matlab/create_sliceintersections.m index a42142d9..d4179877 100644 --- a/src/matlab/create_sliceintersections.m +++ b/src/matlab/create_sliceintersections.m @@ -7,80 +7,101 @@ function create_sliceintersections(mappedgrid) ve_names = {'xe','ye','ze'}; isect_names = {'yzIntersect','xzIntersect','xyIntersect'}; -for idir = 1:3, - for jdir = (idir+1):3, - slices_idir = get_slices(sdirs{idir}); - slices_jdir = get_slices(sdirs{jdir}); - - % Now consider all pair-wise combinations of these two sets of slices - for n = 1:length(slices_idir); - for m = 1:length(slices_jdir); - - % cell2mat doesn't seem to exist on some versions of Matlab - % so I have copied my matlab version to ../matlab2 files. - % pvec = cell2mat_claw(slices_idir{n}); - % qvec = cell2mat_claw(slices_jdir{m}); - - pvec = []; - qvec = []; - % Need a big list of all patches we can intersect - for ll = 1:length(slices_idir{n}), - pvec = [pvec slices_idir{n}{ll}]; - end; - - for ll = 1:length(slices_jdir{m}), - qvec = [qvec slices_jdir{m}{ll}]; - end; - - for k = 1:length(pvec), - for l = 1:length(qvec), - p = pvec(k); - q = qvec(l); - - udata_p = get(p,'UserData'); - udata_q = get(q,'UserData'); - - np = udata_p.grid_number; - nq = udata_q.grid_number; - if (np == nq) - % This sets two of the coordinates for the line we are going - % to draw - sval_p = udata_p.sval; - sval_q = udata_q.sval; - - kdir = (1:3); - kdir([idir; jdir]) = []; - - % Third direction is all that is left - ke = getfield(udata_p,ve_names{kdir}); - vdata{idir} = sval_p + 0*ke; - vdata{jdir} = sval_q + 0*ke; - vdata{kdir} = ke; - - if (mappedgrid == 1) - [xdata, ydata, zdata] = mapc2p(vdata{:}); - else - [xdata, ydata, zdata] = deal(vdata{:}); - end; - - h = line('XData',xdata, 'YData',ydata,'ZData',zdata); - - sp = getfield(udata_p,isect_names{kdir}); - sp(end+1).line = h; - sp(end+1).sharedPatch = q; - - sq = getfield(udata_q,isect_names{kdir}); - sq(end+1).line = h; - sq(end+1).sharedPatch = p; - - udata_p = setfield(udata_p,isect_names{kdir},sp); - udata_q = setfield(udata_q,isect_names{kdir},sq); - set(p,'UserData',udata_p); - set(q,'UserData',udata_q); - end; % np == nq - end; % (l) loop over all patches - end; % (k) loop over all patches - end; % (m) loop over all slices - end; % (n) loop over all slices - end; % jdir -end; % idir +vdata = cell(3,1); + +for idir = 1:3 + for jdir = (idir+1):3 + slices_idir = get_slices(sdirs{idir}); + slices_jdir = get_slices(sdirs{jdir}); + + % Now consider all pair-wise combinations of these two sets of slices + for n = 1:length(slices_idir) + for m = 1:length(slices_jdir) + + l1 = 0; + for ll = 1:length(slices_idir{n}) + l1 = l1 + numel(slices_idir{n}{ll}); + end + + pvec = zeros(l1,1); + % Need a big list of all patches we can intersect + k = 1; + for ll = 1:length(slices_idir{n}) + ps = slices_idir{n}{ll}; + for jj = 1:length(ps) + pvec(k) = ps(jj); + k = k + 1; + end + end + + l2 = 0; + for ll = 1:length(slices_jdir{m}) + l2 = l2 + numel(slices_jdir{m}{ll}); + end + + qvec = zeros(l2,1); + % Need a big list of all patches we can intersect + k = 1; + for ll = 1:length(slices_jdir{m}) + ps = slices_jdir{m}{ll}; + for jj = 1:length(ps) + qvec(k) = ps(jj); + k = k + 1; + end + end + + for k = 1:length(pvec) + for l = 1:length(qvec) + p = pvec(k); + q = qvec(l); + + udata_p = get(p,'UserData'); + udata_q = get(q,'UserData'); + + np = udata_p.grid_number; + nq = udata_q.grid_number; + if (np == nq) + % This sets two of the coordinates for the line we are going + % to draw + sval_p = udata_p.sval; + sval_q = udata_q.sval; + + kdir = (1:3); + kdir([idir; jdir]) = []; + + % Third direction is all that is left + % ke = getfield(udata_p,ve_names{kdir}); + ke = udata_p.(ve_names{kdir}); + vdata{idir} = sval_p + 0*ke; + vdata{jdir} = sval_q + 0*ke; + vdata{kdir} = ke; + + if (mappedgrid == 1) + [xdata, ydata, zdata] = mapc2p(vdata{:}); + else + [xdata, ydata, zdata] = deal(vdata{:}); + end + + h = line('XData',xdata, 'YData',ydata,'ZData',zdata); + + % sp = getfield(udata_p,isect_names{kdir}); + sp = udata_p.(isect_names{kdir}); + sp(end+1).line = h; + sp(end+1).sharedPatch = q; + + %sq = getfield(udata_q,isect_names{kdir}); + sq = udata_q.(isect_names{kdir}); + sq(end+1).line = h; + sq(end+1).sharedPatch = p; + + udata_p.(isect_names{kdir}) = sp; + udata_q.(isect_names{kdir}) = sq; + set(p,'UserData',udata_p); + set(q,'UserData',udata_q); + end % np == nq + end % (l) loop over all patches + end % (k) loop over all patches + end % (m) loop over all slices + end % (n) loop over all slices + end % jdir +end % idir diff --git a/src/matlab/drawcontourlines.m b/src/matlab/drawcontourlines.m index c44ff7a2..a6b7c360 100644 --- a/src/matlab/drawcontourlines.m +++ b/src/matlab/drawcontourlines.m @@ -28,38 +28,38 @@ function drawcontourlines(v, sdir,snum) % See also SHOWCONTOURLINES, HIDECONTOURLINES. if (nargin < 2) - sdirs = {'x', 'y','z'}; + sdirs = {'x', 'y','z'}; else - sdirs = {sdir}; -end; + sdirs = {sdir}; +end -for idir = 1:length(sdirs), - sdir = sdirs{idir}; - slices = get_slices(sdir); - if (nargin < 3) - snum = 1:length(slices); - end; - for ns = 1:length(snum), - n = snum(ns); - if (n < 1 | n > length(slices)) - continue; - end; - slice = slices{n}; - for level = 1:length(slice), - pvec = slice{level}; - for k = 1:length(pvec); - p = pvec(k); - udata = get(p,'UserData'); - delete(udata.contourLines); - [xc_like, yc_like, zc_like] = ... - get_xyzlike(udata.xc,udata.yc,udata.zc,sdir); - c = contourc(yc_like,zc_like,udata.q,v); - udata.contourLines = create_clines(c,udata.sval,... - sdir,udata.mappedgrid,udata.manifold,udata.blockno); - set(p,'UserData',udata); - set_cline_visibility(p); - end; - end; - mask_patches_all(slice); % To mask any contourlines. - end; -end; +for idir = 1:length(sdirs) + sdir = sdirs{idir}; + slices = get_slices(sdir); + if (nargin < 3) + snum = 1:length(slices); + end + for ns = 1:length(snum) + n = snum(ns); + if (n < 1 || n > length(slices)) + continue; + end + slice = slices{n}; + for level = 1:length(slice) + pvec = slice{level}; + for k = 1:length(pvec) + p = pvec(k); + udata = get(p,'UserData'); + delete(udata.contourLines); + [~, yc_like, zc_like] = ... + get_xyzlike(udata.xc,udata.yc,udata.zc,sdir); + c = contourc(yc_like,zc_like,udata.q,v); + udata.contourLines = create_clines(c,udata.sval,... + sdir,udata.mappedgrid,udata.manifold,udata.blockno); + set(p,'UserData',udata); + set_cline_visibility(p); + end + end + mask_patches_all(slice); % To mask any contourlines. + end +end diff --git a/src/matlab/mask_patch.m b/src/matlab/mask_patch.m index 15909476..5c00b16d 100644 --- a/src/matlab/mask_patch.m +++ b/src/matlab/mask_patch.m @@ -19,13 +19,13 @@ function mask_patch(p,sdir,xlow,xhigh,ylow,yhigh,zlow,zhigh) s = ds/2; if (s > min([udata.dx, udata.dy, udata.dz])) - error('mask_patch : s is too big for dx, dy, dz\n'); -end; + error('mask_patch : s is too big for dx, dy, dz\n'); +end vc = udata.cartCoords; -[xdata_like,ydata_like, zdata_like] = get_xyzlike(vc(:,1),vc(:,2), vc(:,3),sdir); -[xlow_like, ylow_like, zlow_like] = get_xyzlike(xlow,ylow,zlow,sdir); -[xhigh_like, yhigh_like, zhigh_like] = get_xyzlike(xhigh,yhigh,zhigh,sdir); +[~,ydata_like, zdata_like] = get_xyzlike(vc(:,1),vc(:,2), vc(:,3),sdir); +[~, ylow_like, zlow_like] = get_xyzlike(xlow,ylow,zlow,sdir); +[~, yhigh_like, zhigh_like] = get_xyzlike(xhigh,yhigh,zhigh,sdir); nan_mask = ydata_like > ylow_like+s & ydata_like < yhigh_like-s & ... zdata_like > zlow_like+s & zdata_like < zhigh_like-s; diff --git a/src/matlab/plotclaw1.m b/src/matlab/plotclaw1.m index 2c635f89..355879a1 100644 --- a/src/matlab/plotclaw1.m +++ b/src/matlab/plotclaw1.m @@ -27,8 +27,8 @@ clawdim = 1; -disp(' ') -disp('plotclaw1 plots 1d results from clawpack') +fprintf('\n') +fprintf('plotclaw1 plots 1d results from clawpack\n') set_value('NoQuery','NoQuery',0); @@ -36,36 +36,36 @@ % set plotting parameters: whichfile = which('setplot1'); if strcmp(whichfile,'') - disp('*** No setplot1 file found') + fprintf('*** No setplot1 file found.\n') else if (NoQuery == 0) - inp = input(['Execute setplot1 (default = yes)? '],'s'); + inp = input('Execute setplot1 (default = yes)? ','s'); if (isempty(inp)) inp = 'y'; end else inp = 'y'; end - - inpd = findstr('y',lower(inp)); + + inpd = strfind('y',lower(inp)); if (inpd == 1) setplot1 - disp(' ') - disp(['Executing m-script ' whichfile]) + fprintf('\n') + fprintf(['Executing m-script %s\n' whichfile]) end end -disp(' ') +fprintf('\n') % the file setprob.m can be used to set up any necessary physical parameters % or desired values of plotting parameters for this particular problem. whichfile = which('setprob'); if strcmp(whichfile,'') - %disp('*** No setprob file found') + fprintf('*** No setprob file found\n') else - disp(['Executing m-script ' whichfile]) - disp(' ') - setprob + fprintf(['Executing m-script %s\n' whichfile]) + fprintf('\n') + setprob(); end @@ -75,35 +75,36 @@ Frame = -1; % initialize frame counter -if ~exist('MaxFrames') - disp('MaxFrames parameter not set... you may need to execute setplot1') - return +if ~exist('MaxFrames','var') + fprintf('MaxFrames parameter not set... you may need to execute setplot1\n') + return end set_value('frameinc','plot_interval',1); set_value('outputdir','OutputDir','./'); set_value('outputflag','OutputFlag','ascii'); +set_value('readblocknumber','ReadBlockNumber',0); amrdata = []; while Frame <= MaxFrames - - % pause for input from user to determine if we go to next frame, - % look at data, or skip around. This may reset Frame counter. - - Frame_old = Frame; - queryframe; % changes value of Frame - - if (query_quit) - break; - end - - if (Frame ~= Frame_old | isempty(amrdata)) - [amrdata,t] = readamrdata(clawdim,Frame,outputdir,outputflag); - end; - - % produce the plot: - - plotframe1 % routine claw/matlab/plotframe1.m - % does the plotting for this frame - + + % pause for input from user to determine if we go to next frame, + % look at data, or skip around. This may reset Frame counter. + + Frame_old = Frame; + queryframe % changes value of Frame + + if (query_quit) + break + end + + if (Frame ~= Frame_old || isempty(amrdata)) + [amrdata,t] = readamrdata(clawdim,Frame,outputdir,outputflag); + end + + % produce the plot: + + plotframe1 % routine claw/matlab/plotframe1.m + % does the plotting for this frame + end % main loop on frames diff --git a/src/matlab/plotclaw2.m b/src/matlab/plotclaw2.m index 8ad0ace2..02206c51 100644 --- a/src/matlab/plotclaw2.m +++ b/src/matlab/plotclaw2.m @@ -28,32 +28,33 @@ clawdim = 2; -disp(' ') -disp('plotclaw2 plots 2d results from clawpack or amrclaw') +fprintf('\n') +fprintf('plotclaw2 plots 2d results from clawpack or amrclaw\n') set_value('NoQuery','NoQuery',0); +set_value('maxlevels','MaxLevels',6); % set plotting parameters: whichfile = which('setplot2'); if strcmp(whichfile,'') - disp('*** No setplot2 file found') + fprintf('*** No setplot2 file found\n') else if (NoQuery == 0) - inp = input(['Execute setplot2 (default = yes)? '],'s'); + inp = input('Execute setplot2 (default = yes)? ','s'); if (isempty(inp)) inp = 'y'; end else inp = 'y'; - end; - inpd = findstr('y',lower(inp)); + end + inpd = strfind('y',lower(inp)); if (inpd == 1) setplot2; - disp(['Executing m-script ' whichfile]) - disp(' ') + fprintf('Executing m-script %s\n',whichfile) + fprintf('\n') end end -disp(' ') +fprintf('\n'); % the file setprob.m can be used to set up any necessary physical parameters % or desired values of plotting parameters for this particular problem. @@ -62,17 +63,17 @@ if strcmp(whichfile,'') %disp('*** No setprob file found') else - disp(['Executing m-script ' whichfile]) - disp(' ') - setprob + fprintf('Executing m-script %s\n',whichfile) + fprintf('\n') + setprob(); end %============================================= % MAIN LOOP ON FRAMES: %============================================= -if ~exist('MaxFrames') - disp('MaxFrames parameter not set... you may need to execute setplot2') +if ~exist('MaxFrames','var') + fprintf('MaxFrames parameter not set... you may need to execute setplot2\n') return end @@ -94,9 +95,9 @@ if (query_quit) break; end - if (Frame ~= Frame_old | isempty(amrdata)) + if (Frame ~= Frame_old || isempty(amrdata)) [amrdata,t] = readamrdata(clawdim,Frame,outputdir,outputflag,... - outputprefix,readblocknumber); + outputprefix,readblocknumber); end plotframe2; diff --git a/src/matlab/plotclaw3.m b/src/matlab/plotclaw3.m index 15d62433..1b8cf186 100644 --- a/src/matlab/plotclaw3.m +++ b/src/matlab/plotclaw3.m @@ -32,24 +32,24 @@ % set plotting parameters: whichfile = which('setplot3'); if strcmp(whichfile,'') - disp('*** No setplot3 file found') + fprintf('*** No setplot3 file found\n'); else if (NoQuery == 0) - inp = input(['Execute setplot3 (default = yes)? '],'s'); + inp = input('Execute setplot3 (default = yes)? ','s'); if (isempty(inp)) inp = 'y'; end else inp = 'y'; - end; - inpd = findstr('y',lower(inp)); + end + inpd = strfind('y',lower(inp)); if (inpd == 1) setplot3; - disp(['Executing m-script ' whichfile]) - disp(' ') + fprintf('Executing m-script %s\n',whichfile) + fprintf('\n') end end -disp(' ') +fprintf('\n'); % the file setprob.m can be used to set up any necessary physical parameters % or desired values of plotting parameters for this particular problem. @@ -57,22 +57,23 @@ whichfile = which('setprob'); if strcmp(whichfile,'') %disp('*** No setprob file found') - else - disp(['Executing m-script ' whichfile]) - disp(' ') - setprob - end +else + fprintf('Executing m-script %s\n',whichfile); + fprintf('\n') + setprob; +end %============================================= % MAIN LOOP ON FRAMES: %============================================= -if ~exist('MaxFrames') - disp('MaxFrames parameter not set... you may need to execute setplot3') - return +if ~exist('MaxFrames','var') + fprintf('MaxFrames parameter not set... you may need to execute setplot3\n'); + return end set_value('frameinc','plot_interval',1); +set_value('maxlevels','MaxLevels',6); set_value('outputdir','OutputDir','./'); set_value('outputflag','OutputFlag','ascii'); set_value('outputprefix','plot_prefix','pltstate'); @@ -81,18 +82,22 @@ clear amrdata; Frame = -frameinc; % initialize frame counter while Frame <= MaxFrames - - % pause for input from user to determine if we go to next frame, - % look at data, or skip around. This may reset Frame counter. - - old_Frame = Frame; - queryframe % this changes value of Frame - - if (old_Frame ~= Frame | isempty(amrdata)) - [amrdata,t] = readamrdata(clawdim,Frame,outputdir,outputflag,... - outputprefix,readblocknumber); - end; - - plotframe3; - + + % pause for input from user to determine if we go to next frame, + % look at data, or skip around. This may reset Frame counter. + + old_Frame = Frame; + queryframe % this changes value of Frame + + if (query_quit) + break; + end + + if (old_Frame ~= Frame || isempty(amrdata)) + [amrdata,t] = readamrdata(clawdim,Frame,outputdir,outputflag,... + outputprefix,readblocknumber); + end + + plotframe3; + end % main loop on frames diff --git a/src/matlab/plotframe2.m b/src/matlab/plotframe2.m index a5e92106..4f24b320 100644 --- a/src/matlab/plotframe2.m +++ b/src/matlab/plotframe2.m @@ -21,339 +21,338 @@ % % See SETPLOT, PLOTCLAW2, MANIFOLD. -if (~exist('amrdata')) - error('*** plotframe2 : ''amrdata'' does not exist; Call readamrdata or plotclaw2'); -end; +if (~exist('amrdata','var')) + error('*** plotframe2 : ''amrdata'' does not exist; Call readamrdata or plotclaw2'); +end -if length(amrdata) == 0 - disp(' ') - disp(['Frame ',num2str(Frame),'(',outputflag,') does not exist ***']); - disp(' ') - return; +if isempty(amrdata) + fprintf('\n'); + fprintf('Frame %d (%s) does not exist\n',Frame,outputflag); + fprintf('\n'); + return; end -disp(' '); -disp(['Frame ',num2str(Frame),' at time t = ',num2str(t)]); +fprintf('\n'); +fprintf('Frame %d at time = %g\n',Frame,t); -if exist('beforeframe')==2 - % make an m-file with this name for any other commands you - % want executed before drawing each frame, for example - % if you want to use axes to specify exactly where the - % plot will be in the window, aspect ratio, etc. - beforeframe; +if exist('beforeframe','file') + % make an m-file with this name for any other commands you + % want executed before drawing each frame, for example + % if you want to use axes to specify exactly where the + % plot will be in the window, aspect ratio, etc. + beforeframe; end set_value('forestclaw','ForestClaw',0); % Do some checking to make sure input is right.. if (PlotType <= 3) - - set_value('mappedgrid','MappedGrid',0); - if (mappedgrid == 1 & ~exist('mapc2p')) - error('*** MappedGrid = 1 but no ''mapc2p'' function was found.'); - end; - - set_value('manifold','Manifold',0); - if (manifold == 1 & ~exist('mapc2m')) - error('*** Manifold = 1, but no ''mapc2m'' function was found.'); - end; - - if (manifold == 1) - set_value('view_arg','UserView',3); - else - set_value('view_arg','UserView',2); - end; - - set_value('underoverflag','ShowUnderOverShoots',0); - if (underoverflag == 1 & ~exist('underover')) - error(['*** ShowUnderOverShoots = 1, but no ''underover'' ',... - 'function was found.']); - end - - set_value('plotparallelpartitions','PlotParallelPartitions',0); - - set_value('usercolormapping','UserColorMapping',0); - if (usercolormapping == 1 & ~exist('setcolors')) - error(['*** UserColorMapping = 1, but no ''setcolors'' ',... - 'function was found.']); - end - - if (underoverflag == 1) - colormapping = 'underover'; - elseif (plotparallelpartitions == 1) - colormapping = 'parallelpartitions'; - elseif (usercolormapping == 1) - colormapping = 'usercolormapping'; - else - colormapping = 'default'; - end; - - set_value('cvalues','ContourValues',[]); - if (length(cvalues) == 1 & length(amrdata) > 1) - disp(' '); - disp(' *** Warning: Contour values will be chosen independently'); - disp(' *** on each mesh patch. Set ContourValues to a vector of values'); - disp(' *** to insure contour lines match across mesh boundaries.'); - end; - - - if (PlotType == 2) - if (isempty(cvalues)) - disp(' '); - disp(' *** Warning : You have specified PlotType == 2, but have'); - disp(' *** not set any contour values. Set ContourValues to a '); - disp(' *** vector of values or integer (number of contour lines). '); - end; - end; - - - if (PlotType == 3) - % Change color map for schlieren plot - % colormap(flipud(gray(2048)).^5) - colormap(flipud(gray(2048)).^10) - - if (~isempty(cvalues)) - disp(' '); - disp(' *** Warning : ContourValues is set to a non-empty matrix.'); - disp(' *** Schlieren data will be contoured.'); - end; - end; - - xscoords = []; - yscoords = []; - zscoords = 0; - - sliceCoords = {xscoords, yscoords, zscoords}; - sdir = {'x','y','z'}; - - % Stored in UserData of current figure (i.e. gcf); - clear_amrplot; - create_amrplot(MaxLevels,xscoords, yscoords, zscoords); - newplot; % Use this to deal with any hold on/ hold off issues. - -end; % end of input checking for PlotType <= 3 + + set_value('mappedgrid','MappedGrid',0); + if (mappedgrid == 1 && ~exist('mapc2p','file')) + error('*** MappedGrid = 1 but no ''mapc2p'' function was found.'); + end + + set_value('manifold','Manifold',0); + if (manifold == 1 && ~exist('mapc2m','file')) + error('*** Manifold = 1, but no ''mapc2m'' function was found.'); + end + + if (manifold == 1) + set_value('view_arg','UserView',3); + else + set_value('view_arg','UserView',2); + end + + set_value('underoverflag','ShowUnderOverShoots',0); + if (underoverflag == 1 && ~exist('underover','file')) + error(['*** ShowUnderOverShoots = 1, but no ''underover'' ',... + 'function was found.']); + end + + set_value('plotparallelpartitions','PlotParallelPartitions',0); + + set_value('usercolormapping','UserColorMapping',0); + if (usercolormapping == 1 && ~exist('setcolors','file')) + error(['*** UserColorMapping = 1, but no ''setcolors'' ',... + 'function was found.']); + end + + if (underoverflag == 1) + colormapping = 'underover'; + elseif (plotparallelpartitions == 1) + colormapping = 'parallelpartitions'; + elseif (usercolormapping == 1) + colormapping = 'usercolormapping'; + else + colormapping = 'default'; + end + + set_value('cvalues','ContourValues',[]); + if (length(cvalues) == 1 && length(amrdata) > 1) + fprintf('\n'); + fprintf(' *** Warning: Contour values will be chosen independently\n'); + disp(' *** on each mesh patch. Set ContourValues to a vector of values'); + disp(' *** to insure contour lines match across mesh boundaries.'); + end + + + if (PlotType == 2) + if (isempty(cvalues)) + disp(' '); + disp(' *** Warning : You have specified PlotType == 2, but have'); + disp(' *** not set any contour values. Set ContourValues to a '); + disp(' *** vector of values or integer (number of contour lines). '); + end + end + + + if (PlotType == 3) + % Change color map for schlieren plot + % colormap(flipud(gray(2048)).^5) + colormap(flipud(gray(2048)).^10) + + if (~isempty(cvalues)) + fprintf('\n'); + fprintf(' *** Warning : ContourValues is set to a non-empty matrix.\n'); + fprintf(' *** Schlieren data will be contoured.\n'); + end + end + + xscoords = []; + yscoords = []; + zscoords = 0; + + sliceCoords = {xscoords, yscoords, zscoords}; + sdir = {'x','y','z'}; + + % Stored in UserData of current figure (i.e. gcf); + clear_amrplot; + create_amrplot(MaxLevels,xscoords, yscoords, zscoords); + newplot; % Use this to deal with any hold on/ hold off issues. + +end % end of input checking for PlotType <= 3 if PlotType == 4 - % ------------- - % Scatter plots - % ------------- - - view_arg = 2; - - set_value('usermap1d','UserMap1d',0); - if (usermap1d == 1) - if (~exist('map1d')) - error('*** You have set UserMap1d=1, but no ''map1d'' function was found'); - end; - else - if ((~exist('x0') | ~exist('y0'))) - str = sprintf(['*** plotframe2 : (x0,y0) must be defined before you\n',... - '*** can do a line or scatter plot, or set UserMap1d and create\n',... - '*** a function ''map1d''.']); - error(str); - end; - end; - - set_value('mappedgrid','MappedGrid',0); - - if (exist('ScatterStyle')) - pstyle = ScatterStyle; - elseif (exist('LineStyle')) - pstyle = LineStyle; - else - error([' *** plotframe2 : Set either ''ScatterStyle'' or ',... - '''LineStyle''.']); - end; - - if (~iscell(pstyle)) - error(['*** plotframe2 : ScatterStyle or LineStyle must be',... - 'cell matrices. Use ''setplotstyle'' to set either of these ',... - 'variables']); - end; - - [linestyle,linecolors,markerstyle] = get_plotstyle(pstyle,MaxLevels); - - clear_amrplot; - create_amrplot(MaxLevels); - newplot; % deals with some hold on/off issues... - -end; % end of input checking for PlotType == 4 + % ------------- + % Scatter plots + % ------------- + + view_arg = 2; + + set_value('usermap1d','UserMap1d',0); + if (usermap1d == 1) + if (~exist('map1d','file')) + error('*** You have set UserMap1d=1, but no ''map1d'' function was found'); + end + else + if ((~exist('x0','var') || ~exist('y0','var'))) + error(['*** plotframe2 : (x0,y0) must be defined before you\n',... + '*** can do a line or scatter plot, or set UserMap1d and create\n',... + '*** a function ''map1d''.']); + end + end + + set_value('mappedgrid','MappedGrid',0); + + if (exist('ScatterStyle','var')) + pstyle = ScatterStyle; + elseif (exist('LineStyle','var')) + pstyle = LineStyle; + else + error([' *** plotframe2 : Set either ''ScatterStyle'' or ',... + '''LineStyle''.']); + end + + if (~iscell(pstyle)) + error(['*** plotframe2 : ScatterStyle or LineStyle must be',... + 'cell matrices. Use ''setplotstyle'' to set either of these ',... + 'variables']); + end + + [linestyle,linecolors,markerstyle] = get_plotstyle(pstyle,MaxLevels); + + clear_amrplot; + create_amrplot(MaxLevels); + newplot; % deals with some hold on/off issues... + +end % end of input checking for PlotType == 4 if (PlotType == 5) - set_value('mappedgrid','MappedGrid',0); - if (mappedgrid == 1 & ~exist('mapc2p')) - error('*** MappedGrid = 1 but no ''mapc2p'' function was found.'); - end; - - set_value('manifold','Manifold',0); - if (manifold == 1 & ~exist('mapc2m')) - error('*** Manifold = 1, but no ''mapc2m'' function was found.'); - end; - set_value('view_arg','UserView',3); -end; + set_value('mappedgrid','MappedGrid',0); + if (mappedgrid == 1 && ~exist('mapc2p','file')) + error('*** MappedGrid = 1 but no ''mapc2p'' function was found.'); + end + + set_value('manifold','Manifold',0); + if (manifold == 1 && ~exist('mapc2m','file')) + error('*** Manifold = 1, but no ''mapc2m'' function was found.'); + end + set_value('view_arg','UserView',3); +end qmin = []; qmax = []; -ncells = []; +ncells = zeros(maxlevels,1); %============================================= % MAIN LOOP ON GRIDS FOR THIS FRAME: %============================================= ngrids = length(amrdata); % length of structure array -for ng = 1:ngrids, - - gridno = amrdata(ng).gridno; - blockno = amrdata(ng).blockno; % == 0 if there is only one block - level = amrdata(ng).level; - if (forestclaw) - level = level + 1; % ForestClaw levels start at 0 - end - - % if we're not plotting data at this level, skip to next grid - if (PlotData(level) == 0) - continue; - end - - % Set block number for multi-block calculations. - set_blocknumber(blockno); - if (forestclaw) - mpirank = amrdata(ng).mpirank; - set_mpirank(mpirank); - else - set_mpirank(0); - end - - mx = amrdata(ng).mx; - my = amrdata(ng).my; - - xlow = amrdata(ng).xlow; - ylow = amrdata(ng).ylow; - - dx = amrdata(ng).dx; - dy = amrdata(ng).dy; - - xedge = xlow + (0:mx)*dx; - yedge = ylow + (0:my)*dy; - - xcenter = xedge(1:mx) + dx/2; - ycenter = yedge(1:my) + dy/2; - - % for compatibility with old matlab41/plotframe2 convention: - x = xcenter; - y = ycenter; - - % read q data: - data = amrdata(ng).data; - - - data = data'; - - if (UserVariable == 1) - % User has supplied a function to convert original q variables to - % the variable which is to be plotted, e.g. Mach number, entropy. - qdata = feval(UserVariableFile,data); - q = reshape(qdata,mx,my); - else - q = reshape(data(:,mq),mx,my); - end - - amrdata(ng).q = q; - - % q must be permuted so that it matches the dimensions of xcm,ycm,zcm - % created by meshgrid. - if (PlotType == 3) - % Scheieren plot; we plot the gradient, not the values. - [qx,qy] = gradient(q,dx,dy); - qs = sqrt(qx.^2 + qy.^2); - qmesh = qs'; - else - qmesh = q'; - end; - - % minimum over all grids at this time, but not necessarily on slice - % shown. - qmin = min([qmin,min(q(:))]); - qmax = max([qmax,max(q(:))]); - - % keep count of how many cells at this refinement level: - if length(ncells) < level - ncells(level) = 0; - end; - ncells(level) = ncells(level) + mx*my; - - % ----------------------------------------------- - % plot commands go here - if PlotType <= 3 - - % Add amr patch of manifold into current plot. +for ng = 1:ngrids + + gridno = amrdata(ng).gridno; + blockno = amrdata(ng).blockno; % == 0 if there is only one block + level = amrdata(ng).level; if (forestclaw) - sval = -level; % Subtract one so that levels correspond to levels - else - sval = level; + level = level + 1; % ForestClaw levels start at 0 end - zedge = [sval sval]; - zcenter = [sval sval]; - sdir = 'z'; - snum = 1; % only one slice in 2d plot - % only mask patches underneath if we are plotting a Manifold, or - % using ForestClaw + + % if we're not plotting data at this level, skip to next grid + if (PlotData(level) == 0) + continue; + end + + % Set block number for multi-block calculations. + set_blocknumber(blockno); if (forestclaw) - maskflag = 0; + mpirank = amrdata(ng).mpirank; + set_mpirank(mpirank); else - maskflag = (manifold == 1); + set_mpirank(0); end - add_patch2slice(sdir,sval,snum,xcenter,ycenter,zcenter, ... - xedge,yedge,zedge,qmesh,level,cvalues,mappedgrid,manifold,... - maskflag,ng,blockno,colormapping); - end; % end of plotting for PlotType == 3 - - if (PlotType == 4) - % 1d Line plots - - [xcm,ycm] = meshgrid(xcenter,ycenter); - - if (usermap1d == 1) - % Users should call mapc2p from inside of map1d. - [rvec,qvec] = map1d(xcm,ycm,qmesh); - [rs,cs] = size(rvec); - [rq,cq] = size(qvec); - if (cs > 1 | cq > 1) - error(['plotframe2 : map1d can only return single columns vectors ',... - 'for s or q']); - end; + + mx = amrdata(ng).mx; + my = amrdata(ng).my; + + xlow = amrdata(ng).xlow; + ylow = amrdata(ng).ylow; + + dx = amrdata(ng).dx; + dy = amrdata(ng).dy; + + xedge = xlow + (0:mx)*dx; + yedge = ylow + (0:my)*dy; + + xcenter = xedge(1:mx) + dx/2; + ycenter = yedge(1:my) + dy/2; + + % for compatibility with old matlab41/plotframe2 convention: + x = xcenter; + y = ycenter; + + % read q data: + data = amrdata(ng).data; + + + data = data'; + + if (UserVariable == 1) + % User has supplied a function to convert original q variables to + % the variable which is to be plotted, e.g. Mach number, entropy. + qdata = feval(UserVariableFile,data); + q = reshape(qdata,mx,my); else - if (mappedgrid == 1) - [xpm,ypm] = mapc2p(xcm,ycm); - r = sqrt((xpm - x0).^2 + (ypm - y0).^2); - else - r = sqrt((xcm-x0).^2 + (ycm - y0).^2); - end - rvec = reshape(r,prod(size(r)),1); - qvec = reshape(qmesh,prod(size(qmesh)),1); - end; - - add_line2plot(rvec,qvec,level,markerstyle{level},... - linecolors{level},linestyle{level}); - - end; % end of plotting for PlotType == 4 - - if (PlotType == 5) - if (manifold == 1) - [xpcenter,ypcenter] = mapc2m(xcenter,ycenter); + q = reshape(data(:,mq),mx,my); + end + + amrdata(ng).q = q; + + % q must be permuted so that it matches the dimensions of xcm,ycm,zcm + % created by meshgrid. + if (PlotType == 3) + % Scheieren plot; we plot the gradient, not the values. + [qx,qy] = gradient(q,dx,dy); + qs = sqrt(qx.^2 + qy.^2); + qmesh = qs'; else - xpcenter = xcenter; - ypcenter = ycenter; - end; - h = surf(xpcenter,ypcenter,q); - end; - - if exist('aftergrid') == 2 - % make an m-file with this name for any other commands you - % want executed at the end of drawing each grid - aftergrid; - end; - + qmesh = q'; + end + + % minimum over all grids at this time, but not necessarily on slice + % shown. + qmin = min([qmin,min(q(:))]); + qmax = max([qmax,max(q(:))]); + + % keep count of how many cells at this refinement level: + if length(ncells) < level + ncells(level) = 0; + end + ncells(level) = ncells(level) + mx*my; + + % ----------------------------------------------- + % plot commands go here + if PlotType <= 3 + + % Add amr patch of manifold into current plot. + if (forestclaw) + sval = -level + 1; % Subtract one so that levels correspond to levels + else + sval = level; + end + zedge = [sval sval]; + zcenter = [sval sval]; + sdir = 'z'; + snum = 1; % only one slice in 2d plot + % only mask patches underneath if we are plotting a Manifold, or + % using ForestClaw + if (forestclaw) + maskflag = 0; + else + maskflag = (manifold == 1); + end + add_patch2slice(sdir,sval,snum,xcenter,ycenter,zcenter, ... + xedge,yedge,zedge,qmesh,level,cvalues,mappedgrid,manifold,... + maskflag,ng,blockno,colormapping); + end % end of plotting for PlotType == 3 + + if (PlotType == 4) + % 1d Line plots + + [xcm,ycm] = meshgrid(xcenter,ycenter); + + if (usermap1d == 1) + % Users should call mapc2p from inside of map1d. + [rvec,qvec] = map1d(xcm,ycm,qmesh); + [rs,cs] = size(rvec); + [rq,cq] = size(qvec); + if (cs > 1 || cq > 1) + error(['plotframe2 : map1d can only return single columns vectors ',... + 'for s or q']); + end + else + if (mappedgrid == 1) + [xpm,ypm] = mapc2p(xcm,ycm); + r = sqrt((xpm - x0).^2 + (ypm - y0).^2); + else + r = sqrt((xcm-x0).^2 + (ycm - y0).^2); + end + rvec = reshape(r,numel(r),1); + qvec = reshape(qmesh,nuel(r),1); + end + + add_line2plot(rvec,qvec,level,markerstyle{level},... + linecolors{level},linestyle{level}); + + end % end of plotting for PlotType == 4 + + if (PlotType == 5) + if (manifold == 1) + [xpcenter,ypcenter] = mapc2m(xcenter,ycenter); + else + xpcenter = xcenter; + ypcenter = ycenter; + end + h = surf(xpcenter,ypcenter,q); + end + + if exist('aftergrid','file') + % make an m-file with this name for any other commands you + % want executed at the end of drawing each grid + aftergrid; + end + end % loop on ng (plot commands for each grid) %============================================= @@ -361,42 +360,44 @@ % add title and labels: if UserVariable == 1 - str = sprintf('%s at time %8.4f',UserVariableFile,t); - title(str,'fontsize',15); + str = sprintf('%s at time %8.4f',UserVariableFile,t); + title(str,'fontsize',15); else - str = sprintf('q(%d) at time %8.4f',mq,t); - title(str,'fontsize',15); + str = sprintf('q(%d) at time %8.4f',mq,t); + title(str,'fontsize',15); end if (PlotType <= 3) - setPlotGrid(PlotGrid); - setPlotGridEdges(PlotGridEdges); - if (PlotType == 2) - setslicecolor('w'); - end; - if (manifold == 0 & view_arg == 2) - % To really make sure we get a 2d view and can see all levels. - set(gca,'ZLimMode','auto'); - zlim = get(gca,'zlim'); - % Increase it sligtly to make sure that we can see all the patches - % plotted at the different z levels - zlim = [min(zlim)-0.001 max(zlim)+0.001]; - set(gca,'zlim',zlim); - end; - xlabel('x') - ylabel('y') -end; + setPlotGrid(PlotGrid); + setPlotGridEdges(PlotGridEdges); + if (PlotType == 2) + setslicecolor('w'); + end + if (manifold == 0 && view_arg == 2) + % To really make sure we get a 2d view and can see all levels. + set(gca,'ZLimMode','auto'); + zlim = get(gca,'zlim'); + % Increase it sligtly to make sure that we can see all the patches + % plotted at the different z levels. + % Also, make sure the zero level is visible (so user can add + % curves, etc to plot). + zlim = [min(zlim)-0.001 0+0.001]; + set(gca,'zlim',[-maxlevels-1,0]); + end + xlabel('x') + ylabel('y') +end % Set view point view(view_arg); afterframe_default; -if exist('afterframe')==2 - % make an m-file with this name for any other commands you - % want executed at the end of drawing each frame - % for example to change the axes, or add a curve for a - % boundary - - afterframe; -end; +if exist('afterframe','file') + % make an m-file with this name for any other commands you + % want executed at the end of drawing each frame + % for example to change the axes, or add a curve for a + % boundary + + afterframe; +end diff --git a/src/matlab/plotframe3.m b/src/matlab/plotframe3.m index 6ad17add..84c062b7 100644 --- a/src/matlab/plotframe3.m +++ b/src/matlab/plotframe3.m @@ -24,208 +24,209 @@ % % See SETPLOT, PLOTCLAW3, MAPPEDGRID. -if (~exist('amrdata')) - str = sprintf(['*** Plotframe3 : ''amrdata'' does not exist. call\n',... - '*** readamrdata or plotclaw3']); - error(str); -end; +if ~exist('amrdata','var') + error(['*** Plotframe3 : ''amrdata'' does not exist. call\n',... + '*** readamrdata or plotclaw3']); +end -if length(amrdata) == 0 - disp(' ') - disp(['Frame ',num2str(Frame),'(',outputflag,') does not exist ***']); - disp(' ') - return; -end; +if isempty(amrdata) + fprintf('\n'); + fprintf('Frame %d (%s) does not exist\n',Frame,outputflag); + fprint('\n'); + return; +end -disp(' '); -disp(['Frame ',num2str(Frame),' at time t = ',num2str(t)]); +fprintf('\n'); +fprintf('Frame %d at time t %g\n', Frame, t); -if exist('beforeframe')==2 - % make an m-file with this name for any other commands you - % want executed before drawing each frame, for example - % if you want to use axes to specify exactly where the - % plot will be in the window, aspect ratio, etc. - beforeframe; +if exist('beforeframe','file') + % make an m-file with this name for any other commands you + % want executed before drawing each frame, for example + % if you want to use axes to specify exactly where the + % plot will be in the window, aspect ratio, etc. + beforeframe; end + + % ------------------------------------------------------------------------ % Check PlotType <= 3 inputs (slices) % ------------------------------------------------------------------------ if (PlotType <= 3) - - % The manifold option isn't used in 3d - manifold = 0; - - set_value('mappedgrid','MappedGrid',0); - set_value('view_arg','UserView',3); - set_value('cvalues','ContourValues',[]); - set_value('interpmethod','InterpMethod','*linear'); - - if (mappedgrid == 1) - if (~exist('mapc2p')) - error('*** MappedGrid == 1 but no ''mapc2p'' function was found'); - end; - end; - - set_value('underoverflag','ShowUnderOverShoots',0); - if (underoverflag == 1 & ~exist('underover')) - error(['*** ShowUnderOverShoots = 1, but no ''underover'' ',... - 'function was found.']); - end - - set_value('usercolormapping','UserColorMapping',0); - if (usercolormapping == 1 & ~exist('setcolors')) - error(['*** UserColorMapping = 1, but no ''setcolors'' ',... - 'function was found.']); - end - - if (underoverflag == 1) - colormapping = 'underover'; - elseif (usercolormapping == 1) - colormapping = 'usercolormapping'; - else - colormapping = 'default'; - end; - - if (PlotType == 2) - if (isempty(cvalues)) - disp(' '); - disp(' *** Warning : You have specified PlotType == 2, but have'); - disp(' *** not set any contour values. Set ContourValues to a '); - disp(' *** vector of values or integer (number of contour lines). '); - end; - end; - - % Change color map for schlieren plot - if (PlotType == 3) - - % colormap(flipud(gray(2048)).^5) - colormap(flipud(gray(2048)).^10) - - if (~isempty(cvalues)) - disp(' '); - disp(' *** Warning : ContourValues is set to a non-empty matrix.'); - disp(' *** Schlieren data will be contoured.'); - end; - end; - - set_value('xscoords','xSliceCoords',[]); - set_value('yscoords','ySliceCoords',[]); - set_value('zscoords','zSliceCoords',[]); - - if (length(cvalues) == 1 & (length(amrdata) > 1 | ... - (length(xscoords) + length(yscoords) + length(zscoords) > 1))) - disp(' '); - disp(' *** Warning: Contour values will be chosen independently'); - disp(' *** on each amr patch or slice. Set ContourValues to a '); - disp(' *** vector of values to insure that contour lines match '); - disp(' *** across mesh boundaries or slices.'); - end; - - set_value('isosurfvalues','IsosurfValues',[]); - - m = length(isosurfvalues); - if (m > 0) - % now check input for isosurfaces - - set_value('isoalphas','IsosurfAlphas',1); - isosurfalphas = set_length(isoalphas,m); - - set_value('isc','IsosurfColors','q'); - if (iscell(isc)) - % In case SETPLOTSTYLE was used to to set the isosurfcolors. - isocolors = strvcat(isc); + + % The manifold option isn't used in 3d + manifold = 0; + + set_value('mappedgrid','MappedGrid',0); + set_value('view_arg','UserView',3); + set_value('cvalues','ContourValues',[]); + set_value('interpmethod','InterpMethod','*linear'); + + if (mappedgrid == 1) + if ~exist('mapc2p','file') + error('*** MappedGrid == 1 but no ''mapc2p'' function was found'); + end + end + + set_value('underoverflag','ShowUnderOverShoots',0); + if underoverflag == 1 && ~exist('underover','file') + error(['*** ShowUnderOverShoots = 1, but no ''underover'' ',... + 'function was found.']); + end + + set_value('usercolormapping','UserColorMapping',0); + if (usercolormapping == 1 && ~exist('setcolors','file')) + error(['*** UserColorMapping = 1, but no ''setcolors'' ',... + 'function was found.']); + end + + if (underoverflag == 1) + colormapping = 'underover'; + elseif (usercolormapping == 1) + colormapping = 'usercolormapping'; else - isocolors = isc; - end; - isosurfcolors = set_length(isocolors,m); - - versionstring = version; - version_number = str2num(versionstring(1)); - - for i = 1:length(isosurfalphas) - a = isosurfalphas(i); - has_opengl = strcmp(get(gcf,'Renderer'),'OpenGL'); - if (a < 1 & version_number < 6) - disp(' '); - disp(' *** plotframe3: Transparency for isosurfaces supported only in') - disp(' *** Matlab versions 6 and higher.'); - isosurfalphas(i) = 1; - elseif (a < 1 & has_opengl == 0) - disp(' *** Warning : The graphics renderer OpenGL is not set on your'); - disp(' *** system. Use setopengl to set this renderer. Without it,'); - disp(' *** isosurfaces will not appear transparent'); - isosurfalphas(i) = 1; - end; - end; - else - % These won't be used... - isosurfcolors = ''; - isosurfalphas = []; - end; % end of length(isosurfvalues > 0) - - clear_amrplot; - create_amrplot(MaxLevels,xscoords,yscoords,zscoords,isosurfvalues); - newplot; - -end; % end of PlotType <= 3 + colormapping = 'default'; + end + + if (PlotType == 2) + if (isempty(cvalues)) + disp(' '); + disp(' *** Warning : You have specified PlotType == 2, but have'); + disp(' *** not set any contour values. Set ContourValues to a '); + disp(' *** vector of values or integer (number of contour lines). '); + end + end + + % Change color map for schlieren plot + if (PlotType == 3) + + % colormap(flipud(gray(2048)).^5) + colormap(flipud(gray(2048)).^10) + + if (~isempty(cvalues)) + disp(' '); + disp(' *** Warning : ContourValues is set to a non-empty matrix.'); + disp(' *** Schlieren data will be contoured.'); + end + end + + set_value('xscoords','xSliceCoords',[]); + set_value('yscoords','ySliceCoords',[]); + set_value('zscoords','zSliceCoords',[]); + + if (length(cvalues) == 1 && (length(amrdata) > 1 || ... + (length(xscoords) + length(yscoords) + length(zscoords) > 1))) + disp(' '); + disp(' *** Warning: Contour values will be chosen independently'); + disp(' *** on each amr patch or slice. Set ContourValues to a '); + disp(' *** vector of values to insure that contour lines match '); + disp(' *** across mesh boundaries or slices.'); + end + + set_value('isosurfvalues','IsosurfValues',[]); + + m = length(isosurfvalues); + if (m > 0) + % now check input for isosurfaces + + set_value('isosurfalphas','IsosurfAlphas',1); + isosurfalphas = set_length(isosurfalphas,m); + + set_value('isc','IsosurfColors','q'); + if (iscell(isc)) + % In case SETPLOTSTYLE was used to to set the isosurfcolors. + isocolors = char(isc); + else + isocolors = isc; + end + isosurfcolors = set_length(isocolors,m); + + versionstring = version; + version_number = str2double(versionstring(1)); + + for i = 1:length(isosurfalphas) + a = isosurfalphas(i); + has_opengl = setopengl(); + if (a < 1 && version_number < 6) + fprintf('\n'); + fprintf(' *** plotframe3: Transparency for isosurfaces supported only in\n') + fprintf(' *** Matlab versions 6 and higher.\n'); + isosurfalphas(i) = 1; + elseif (a < 1 && has_opengl) + setopengl(); + elseif (has_opengl == 0) + fprintf(' *** Warning : The graphics renderer OpenGL is not set on your\n'); + fprintf(' *** system. Use setopengl to set this renderer. Without it,\n'); + fprintf(' *** isosurfaces will not appear transparent\n'); + isosurfalphas(i) = 1; + end + end + else + % These won't be used... + isosurfcolors = ''; + isosurfalphas = []; + end % end of length(isosurfvalues > 0) + + clear_amrplot; + create_amrplot(maxlevels,xscoords,yscoords,zscoords,isosurfvalues); + newplot; + +end % end of PlotType <= 3 % ------------------------------------------------------------------------ % Check PlotType == 4 inputs (1d line plots) % ------------------------------------------------------------------------ if (PlotType == 4) - - view_arg = 2; - - set_value('usermap1d','UserMap1d',0); - if (usermap1d == 1) - if (~exist('map1d')) - error('*** You have set UserMap1d=1, but no ''map1d'' function was found'); - end; - else - if ((~exist('x0') | ~exist('y0') | ~exist('z0'))) - str = sprintf(['*** plotframe3 : (x0,y0,z0) must be defined before you\n',... - '*** can do a line or scatter plot. Or, create a ''map1d''',... - ' function and set UserMap1d = 1']); - error(str); - end; - end; - - set_value('mappedgrid','MappedGrid',0); - - if (exist('ScatterStyle')) - pstyle = ScatterStyle; - elseif (exist('LineStyle')) - pstyle = LineStyle; - else - error([' *** plotframe3 : Set either ''ScatterStyle'' or ',... - '''LineStyle''.']); - end; - - if (~iscell(pstyle)) - str = sprintf(['*** plotframe3 : ScatterStyle or LineStyle must be cell \n',... - 'matrices. Use SETPLOTSTYLE to set either of these variables.']); - error(str); - end; - - [linestyle,linecolors,markerstyle] = get_plotstyle(pstyle,MaxLevels); - - clear_amrplot; - create_amrplot(MaxLevels); - newplot; - -end; % end of PlotType == 4 + + view_arg = 2; + + set_value('usermap1d','UserMap1d',0); + if (usermap1d == 1) + if (~exist('map1d','file')) + error('*** You have set UserMap1d=1, but no ''map1d'' function was found'); + end + else + if (~exist('x0','var') || ~exist('y0','var') || ~exist('z0','var')) + error(['*** plotframe3 : (x0,y0,z0) must be defined before you\n',... + '*** can do a line or scatter plot. Or, create a ''map1d''',... + ' function and set UserMap1d = 1']); + end + end + + set_value('mappedgrid','MappedGrid',0); + + if (exist('ScatterStyle','var')) + pstyle = ScatterStyle; + elseif (exist('LineStyle','var')) + pstyle = LineStyle; + else + error([' *** plotframe3 : Set either ''ScatterStyle'' or ',... + '''LineStyle''.']); + end + + if (~iscell(pstyle)) + error(['*** plotframe3 : ScatterStyle or LineStyle must be cell \n',... + 'matrices. Use SETPLOTSTYLE to set either of these variables.']); + end + + [linestyle,linecolors,markerstyle] = get_plotstyle(pstyle,maxlevels); + + clear_amrplot; + create_amrplot(maxlevels); + newplot; + +end % end of PlotType == 4 % ------------------------------------------------------------------------ % Check PlotType == 5 inputs (deprecated isosurface plots) % ------------------------------------------------------------------------ if (PlotType == 5) - disp(' ') - disp(' *** PlotType = 5 is no longer supported as of Clawpack 4.2') - disp(' *** Choose different setting for PlotType in setplot3.m') - disp(' *** Isosurface plots are avaiable with PlotType = 1, 2 or 3') - disp(' *** by setting IsosurfValues and IsosurfColors') - disp(' ') -end; % end of PlotType == 5 + fprintf('\n') + fprintf(' *** PlotType = 5 is no longer supported as of Clawpack 4.2\n') + fprintf(' *** Choose different setting for PlotType in setplot3.m\n') + fprintf(' *** Isosurface plots are avaiable with PlotType = 1, 2 or 3\n') + fprintf(' *** by setting IsosurfValues and IsosurfColors\n') + fprintf('\n') +end % end of PlotType == 5 % End of input check. % ------------------------------------------------------------------ @@ -241,198 +242,198 @@ qmin = []; qmax = []; -ncells = []; +ncells = zeros(maxlevels,1); ngrids = length(amrdata); % length of structure array -for ng = 1:ngrids, - - gridno = amrdata(ng).gridno; - blockno = amrdata(ng).blockno; - level = amrdata(ng).level; - - set_blocknumber(blockno); - - if (PlotData(level) == 0) - % Continue to next patch - continue; - end; - - mx = amrdata(ng).mx; - my = amrdata(ng).my; - mz = amrdata(ng).mz; - - xlow = amrdata(ng).xlow; - ylow = amrdata(ng).ylow; - zlow = amrdata(ng).zlow; - - dx = amrdata(ng).dx; - dy = amrdata(ng).dy; - dz = amrdata(ng).dz; - - data = amrdata(ng).data; - - % Get grid info for plotting grid boxes. - xhigh = xlow + mx*dx; - yhigh = ylow + my*dy; - zhigh = zlow + mz*dz; - - xedge = xlow + (0:mx)*dx; - yedge = ylow + (0:my)*dy; - zedge = zlow + (0:mz)*dz; - - xcenter = xedge(1:mx) + dx/2; - ycenter = yedge(1:my) + dy/2; - zcenter = zedge(1:mz) + dz/2; - - % for compatibility with old matlab41/plotframe3 convention: - x = xcenter; - y = ycenter; - z = zcenter; - - data = data'; - - if (UserVariable == 1) - % User has supplied a function to convert original q variables to - % the variable which is to be plotted, e.g. Mach number, entropy. - qdata = feval(UserVariableFile,data); - q = reshape(qdata,mx,my,mz); - else - q = reshape(data(:,mq),mx,my,mz); - end - - - % q must be permuted so that it matches the dimensions of xcm,ycm,zcm - % created by meshgrid. - if (PlotType == 3) - % Scheieren plot; we plot the gradient, not the values. - [qx,qy,qz] = gradient(q,dx,dy,dz); - qs = sqrt(qx.^2 + qy.^2 + qz.^2); - qmesh = permute(qs,[2,1,3]); - else - qmesh = permute(q,[2,1,3]); - end; - amrdata(ng).q = qmesh; - - % minimum over all grids at this time, but not necessarily on slice - % shown. - qmin = min([qmin,min(min(min(q)))]); - qmax = max([qmax,max(max(max(q)))]); - - % keep count of how many cells at this refinement level: - if length(ncells) < level - ncells(level) = 0; - end; - ncells(level) = ncells(level) + mx*my*mz; - - % ----------------------------------------------- - % plot commands go here - if PlotType <= 3 - - sliceCoords = {xscoords, yscoords, zscoords}; - sdirs = {'x','y','z'}; - maskflag = 1; % Mask out patch areas - - for idir = 1:3, % Loop over directions - slicevals = sliceCoords{idir}; % user specified slice constants - sdir = sdirs{idir}; - - % Loop over all slices that cube of data qmesh might intersect - for n = 1:length(slicevals), - sval = slicevals(n); - [isect,qcm2] = interp_data_3d(xcenter,ycenter,zcenter,... - xedge,yedge,zedge,qmesh,sdir,sval,interpmethod); - if (isect == 1) % qmesh intersects a user specified slice. - % This command adds a patch to the appropriate slice - add_patch2slice(sdir,sval, n,xcenter,ycenter,zcenter,... - xedge,yedge,zedge,qcm2,level,... - cvalues,mappedgrid, manifold,maskflag,ng,blockno,... - colormapping); - end; - end; - end; - hidelevels(find(~PlotData)); - - % Add cube to plot - whether it is visible or not depends on how user - % set PlotCubeEdges. - add_cube2plot(xedge,yedge,zedge,level,mappedgrid); - - % Add isosurfaces to plot - add_isosurface2plot(xcenter,ycenter,zcenter,q,level,... - isosurfvalues,isosurfcolors,isosurfalphas,mappedgrid); - - end; % End plotting for PlotType == 3 - - if (PlotType == 4) - % 1d Line plots - - [xcm,ycm,zcm] = meshgrid(xcenter,ycenter,zcenter); - - if (usermap1d == 1) - % User should call mapc2p from inside of map1d - [rvec,qvec] = map1d(xcm,ycm,zcm,qmesh); - [rs,cs] = size(rvec); - [rq,cq] = size(qvec); - if (cs > 1 | cq > 1) - error(['plotframe3 : map1d can only return single columns vectors ',... - 'for s or q']); - end; +for ng = 1:ngrids + + gridno = amrdata(ng).gridno; + blockno = amrdata(ng).blockno; + level = amrdata(ng).level; + + set_blocknumber(blockno); + + if (PlotData(level) == 0) + % Continue to next patch + continue; + end + + mx = amrdata(ng).mx; + my = amrdata(ng).my; + mz = amrdata(ng).mz; + + xlow = amrdata(ng).xlow; + ylow = amrdata(ng).ylow; + zlow = amrdata(ng).zlow; + + dx = amrdata(ng).dx; + dy = amrdata(ng).dy; + dz = amrdata(ng).dz; + + data = amrdata(ng).data; + + % Get grid info for plotting grid boxes. + xhigh = xlow + mx*dx; + yhigh = ylow + my*dy; + zhigh = zlow + mz*dz; + + xedge = xlow + (0:mx)*dx; + yedge = ylow + (0:my)*dy; + zedge = zlow + (0:mz)*dz; + + xcenter = xedge(1:mx) + dx/2; + ycenter = yedge(1:my) + dy/2; + zcenter = zedge(1:mz) + dz/2; + + % for compatibility with old matlab41/plotframe3 convention: + x = xcenter; + y = ycenter; + z = zcenter; + + data = data'; + + if (UserVariable == 1) + % User has supplied a function to convert original q variables to + % the variable which is to be plotted, e.g. Mach number, entropy. + qdata = feval(UserVariableFile,data); + q = reshape(qdata,mx,my,mz); else - if (mappedgrid == 1) - [xpm,ypm,zpm] = mapc2p(xcm,ycm,zcm); - r = sqrt((xpm - x0).^2 + (ypm - y0).^2 + (zpm - z0).^2); - else - r = sqrt((xcm - x0).^2 + (ycm - y0).^2 + (zcm - z0).^2); - end; - rvec = reshape(r,prod(size(r)),1); - qvec = reshape(qmesh,prod(size(qmesh)),1); - end; - - add_line2plot(rvec,qvec,level,markerstyle{level},... - linecolors{level},linestyle{level}); - - end; % end of plotting for PlotType == 4 - - if exist('aftergrid') == 2 - % make an m-file with this name for any other commands you - % want executed at the end of drawing each grid - aftergrid; - end; - + q = reshape(data(:,mq),mx,my,mz); + end + + + % q must be permuted so that it matches the dimensions of xcm,ycm,zcm + % created by meshgrid. + if (PlotType == 3) + % Scheieren plot; we plot the gradient, not the values. + [qx,qy,qz] = gradient(q,dx,dy,dz); + qs = sqrt(qx.^2 + qy.^2 + qz.^2); + qmesh = permute(qs,[2,1,3]); + else + qmesh = permute(q,[2,1,3]); + end + amrdata(ng).q = qmesh; % Not sure about warning. + + % minimum over all grids at this time, but not necessarily on slice + % shown. + qmin = min([qmin,min(min(min(q)))]); + qmax = max([qmax,max(max(max(q)))]); + + % keep count of how many cells at this refinement level: + if length(ncells) < level + ncells(level) = 0; + end + ncells(level) = ncells(level) + mx*my*mz; + + % ----------------------------------------------- + % plot commands go here + if PlotType <= 3 + + sliceCoords = {xscoords, yscoords, zscoords}; + sdirs = {'x','y','z'}; + maskflag = 1; % Mask out patch areas + + for idir = 1:3 % Loop over directions + slicevals = sliceCoords{idir}; % user specified slice constants + sdir = sdirs{idir}; + + % Loop over all slices that cube of data qmesh might intersect + for n = 1:length(slicevals) + sval = slicevals(n); + [isect,qcm2] = interp_data_3d(xcenter,ycenter,zcenter,... + xedge,yedge,zedge,qmesh,sdir,sval,interpmethod); + if (isect == 1) % qmesh intersects a user specified slice. + % This command adds a patch to the appropriate slice + add_patch2slice(sdir,sval, n,xcenter,ycenter,zcenter,... + xedge,yedge,zedge,qcm2,level,... + cvalues,mappedgrid, manifold,maskflag,ng,blockno,... + colormapping); + end + end + end + hidelevels(find(~PlotData)); + + % Add cube to plot - whether it is visible or not depends on how user + % set PlotCubeEdges. + add_cube2plot(xedge,yedge,zedge,level,mappedgrid); + + % Add isosurfaces to plot + add_isosurface2plot(xcenter,ycenter,zcenter,q,level,... + isosurfvalues,isosurfcolors,isosurfalphas,mappedgrid); + + end % End plotting for PlotType == 3 + + if (PlotType == 4) + % 1d Line plots + + [xcm,ycm,zcm] = meshgrid(xcenter,ycenter,zcenter); + + if (usermap1d == 1) + % User should call mapc2p from inside of map1d + [rvec,qvec] = map1d(xcm,ycm,zcm,qmesh); + [rs,cs] = size(rvec); + [rq,cq] = size(qvec); + if (cs > 1 || cq > 1) + error(['plotframe3 : map1d can only return single columns vectors ',... + 'for s or q']); + end + else + if (mappedgrid == 1) + [xpm,ypm,zpm] = mapc2p(xcm,ycm,zcm); + r = sqrt((xpm - x0).^2 + (ypm - y0).^2 + (zpm - z0).^2); + else + r = sqrt((xcm - x0).^2 + (ycm - y0).^2 + (zcm - z0).^2); + end + rvec = reshape(r,numel(r),1); + qvec = reshape(qmesh,numel(qmesh),1); + end + + add_line2plot(rvec,qvec,level,markerstyle{level},... + linecolors{level},linestyle{level}); + + end % end of plotting for PlotType == 4 + + if exist('aftergrid','file') == 2 + % make an m-file with this name for any other commands you + % want executed at the end of drawing each grid + aftergrid; + end + end % loop on ng (plot commands for each grid) % ----------------------------------------------------- set(gca,'FontName','helvetica'); if UserVariable == 1 - str = sprintf('%s at time %8.4f',UserVariableFile,t); - title(str,'fontsize',15); + str = sprintf('%s at time %8.4f',UserVariableFile,t); + title(str,'fontsize',15); else - str = sprintf('q(%d) at time %8.4f',mq,t); - title(str,'fontsize',15); + str = sprintf('q(%d) at time %8.4f',mq,t); + title(str,'fontsize',15); end % Set user defined grid details, from setplot3.m if (PlotType <= 3) - setPlotGrid(PlotGrid); - setPlotGridEdges(PlotGridEdges); - setPlotCubeEdges(PlotCubeEdges); - if (PlotType == 2) - setslicecolor('w'); - end; - - % Now plot intersections of x-y planes, x-z planes and y-z planes. - create_sliceintersections(mappedgrid); - xlabel('x') - ylabel('y') - zlabel('z') -end; + setPlotGrid(PlotGrid); + setPlotGridEdges(PlotGridEdges); + setPlotCubeEdges(PlotCubeEdges); + if (PlotType == 2) + setslicecolor('w'); + end + + % Now plot intersections of x-y planes, x-z planes and y-z planes. + create_sliceintersections(mappedgrid); + xlabel('x') + ylabel('y') + zlabel('z') +end view(view_arg); afterframe_default; -if exist('afterframe')==2 - afterframe % make an m-file with this name for any other commands you - % want executed at the end of drawing each frame - % for example to change the axes, or add a curve for a - % boundary -end; +if exist('afterframe','file') + afterframe % make an m-file with this name for any other commands you + % want executed at the end of drawing each frame + % for example to change the axes, or add a curve for a + % boundary +end diff --git a/src/matlab/readamrdata.m b/src/matlab/readamrdata.m index 796ce640..8b7ff454 100644 --- a/src/matlab/readamrdata.m +++ b/src/matlab/readamrdata.m @@ -1,4 +1,5 @@ -function [amr,t] = readamrdata(dim,Frame,dir,outputflag,outputprefix,readblocknumber); +function [amr,t] = readamrdata(dim,Frame,dir,outputflag, ... + outputprefix,readblocknumber) % READAMRDATA reads amr data produced by Clawpack output routines. % @@ -15,9 +16,15 @@ % % [...] = READAMRDATA(...,FLAG) uses FLAG to determine what kind of data % to read. FLAG = 'ascii' means read the standard ascii output files -% fort.q00NN produced by Clawpack. FLAG = 'hdf' will read HDF4 files -% produced by Clawpack. FLAG='chombo' will read the HDF5 files produced -% by ChomboClaw. The default option is FLAG='ascii'. +% fort.q00NN produced by Clawpack. +% +% Possible values for FLAG : +% +% 'ascii' : Standard Clawpack output (fort.qXXXX files) +% 'binary' : Binary output (fort.qXXXX and fort.bXXXX files) +% 'forestclaw' : ForestClaw output (reads MPI rank, for example). +% 'hdf' : HDF 4 (deprecated and may not work). +% 'chombo' : Chombo output (deprecated and may not work). % % The output argument AMRDATA is an array of structures whose fields are % some or all of : @@ -50,37 +57,41 @@ % fprintf('Cells at level % d : %d\n',level,lvec(level)); % end; % -% See CLAWGRAPH, CHOMBOCLAW. +% See CLAWGRAPHICS. % if (nargin < 6) - readblocknumber = 0; - if (nargin < 5) - outputprefix = ''; - if (nargin < 4) - outputflag = 'ascii'; - if (nargin < 3) - dir = './'; - end; - end; - end; + readblocknumber = 0; + if (nargin < 5) + outputprefix = ''; + if (nargin < 4) + outputflag = 'ascii'; + if (nargin < 3) + dir = './'; + end + end + end end if (strcmp(dir(end),'/') == 0) dir = [dir, '/']; -end; +end + +outflag = lower(outputflag); +switch outflag + case 'ascii' + [amr,t] = readamrdata_ascii(dim,Frame,dir,readblocknumber); + case 'binary' + [amr,t] = readamrdata_binary(dim,Frame,dir,readblocknumber); + case 'forestclaw' + [amr,t] = readamrdata_forestclaw(dim,Frame,dir); + case 'hdf' + [amr,t] = readamrdata_hdf(dim,Frame,dir); + case 'chombo' + [amr,t] = readamrdata_chombo(dim,Frame,dir,outputprefix); + otherwise + error('readamrdata : ''%s'' is not a valid OutputFlag\n',flag); +end -if (strcmp(lower(outputflag),'ascii') == 1) - [amr,t] = readamrdata_ascii(dim,Frame,dir,readblocknumber); -elseif (strcmp(lower(outputflag),'hdf') == 1) - [amr,t] = readamrdata_hdf(dim,Frame,dir); -elseif (strcmp(lower(outputflag),'chombo') == 1) - [amr,t] = readamrdata_chombo(dim,Frame,dir,outputprefix); -elseif (strcmp(lower(outputflag),'forestclaw') == 1) - [amr,t] = readamrdata_forestclaw(dim,Frame,dir); -else - str = sprintf(['readamrdata : ''%s'' is not a valid OutputFlag. Use ',... - 'flag = ''ascii'' or flag = ''hdf''.'],flag); - error(str); end diff --git a/src/matlab/readamrdata_ascii.m b/src/matlab/readamrdata_ascii.m index b241d138..e7d56a44 100644 --- a/src/matlab/readamrdata_ascii.m +++ b/src/matlab/readamrdata_ascii.m @@ -1,4 +1,4 @@ -function [amr,t] = readamrdata_ascii(dim,Frame,dir,readblocknumber); +function [amr,t] = readamrdata_ascii(dim,Frame,dir,readblocknumber) % Internal routine for Clawpack graphics. @@ -6,19 +6,18 @@ fname = [dir, 'fort.',num2str(n1)]; fname(length(dir)+6) = 't'; -if ~exist(fname) - amr = {}; - t = []; - disp(' '); - disp(['Frame ',num2str(Frame),' (',fname,') does not exist ***']); - disp(' '); - return; +if ~exist(fname,'file') + amr = {}; + t = []; + fprintf('\n'); + fprintf('Frame %d (%s) does not exist ***\n',Frame,fname); + fprintf('\n'); + return; end % Read data from fname = 'fort.tXXXX' fid = fopen(fname); - t = fscanf(fid,'%g',1); fscanf(fid,'%s',1); meqn = fscanf(fid,'%d',1); fscanf(fid,'%s',1); ngrids = fscanf(fid,'%d',1); fscanf(fid,'%s',1); @@ -26,63 +25,85 @@ % change the file name to read the q data: fname(length(dir) + 6) = 'q'; -if ~exist(fname) - amr = {}; - t = []; - return; +if ~exist(fname,'file') + amr = {}; + t = []; + return; end fid = fopen(fname); -disp(['Reading data from ',fname]); - -for ng = 1:ngrids, - - % read parameters for this grid: - - amrdata.gridno = fscanf(fid,'%d',1); fscanf(fid,'%s',1); - amrdata.level = fscanf(fid,'%d',1); fscanf(fid,'%s',1); - if (dim > 1) - if (readblocknumber) - amrdata.blockno = fscanf(fid,'%d',1); fscanf(fid,'%s',1); +fprintf('Reading data from %s ',fname); + +if dim == 1 + amr_s = struct('gridno',[],'level',[],'mx',[],'xlow',[],... + 'dx',[],'data',[]); +elseif dim == 2 + amr_s = struct('gridno',[],'level',[],'blockno',[],... + 'mx',[],'my',[],'xlow',[],'ylow',[],'dx',[],'dy',[],'data',[]); +else + amr_s = struct('gridno',[],'level',[],'blockno',[],... + 'mx',[],'my',[],'mz',[],'xlow',[],'ylow',[],'zlow',[],... + 'dx',[],'dy',[],'dz',[],'data',[],'q',[]); +end +amr = repmat(amr_s,ngrids,1); + +M = floor(ngrids/10); + + +for ng = 1:ngrids + + % read parameters for this grid: + + amrdata.gridno = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + amrdata.level = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + if (dim > 1) + if (readblocknumber) + amrdata.blockno = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + else + amrdata.blockno = 1; + end + end + amrdata.mx = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + if (dim > 1) + amrdata.my = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + if (dim > 2) + amrdata.mz = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + end + end + + amrdata.xlow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + if (dim > 1) + amrdata.ylow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + if (dim > 2) + amrdata.zlow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + end + end + + amrdata.dx = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + if (dim > 1) + amrdata.dy = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + if (dim > 2) + amrdata.dz = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + end + end + + % read q data: + if (dim == 1) + amrdata.data = fscanf(fid,'%g',[meqn,amrdata.mx]); + elseif (dim == 2) + amrdata.data = fscanf(fid,'%g',[meqn,amrdata.mx*amrdata.my]); else - amrdata.blockno = 1; + amrdata.data = fscanf(fid,'%g',[meqn,amrdata.mx*amrdata.my*amrdata.mz]); end - end - amrdata.mx = fscanf(fid,'%d',1); fscanf(fid,'%s',1); - if (dim > 1) - amrdata.my = fscanf(fid,'%d',1); fscanf(fid,'%s',1); - if (dim > 2) - amrdata.mz = fscanf(fid,'%d',1); fscanf(fid,'%s',1); - end; - end; - - amrdata.xlow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); - if (dim > 1) - amrdata.ylow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); - if (dim > 2) - amrdata.zlow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); - end; - end; - - amrdata.dx = fscanf(fid,'%g',1); fscanf(fid,'%s',1); - if (dim > 1) - amrdata.dy = fscanf(fid,'%g',1); fscanf(fid,'%s',1); - if (dim > 2) - amrdata.dz = fscanf(fid,'%g',1); fscanf(fid,'%s',1); - end; - end; - - % read q data: - if (dim == 1) - amrdata.data = fscanf(fid,'%g',[meqn,amrdata.mx]); - elseif (dim == 2) - amrdata.data = fscanf(fid,'%g',[meqn,amrdata.mx*amrdata.my]); - else - amrdata.data = fscanf(fid,'%g',[meqn,amrdata.mx*amrdata.my*amrdata.mz]); - end; - - amr(ng) = amrdata; - + + amr(ng) = amrdata; + + if mod(ng,M) == 0 + fprintf('.'); + end + end +fprintf(' Done\n'); + fclose(fid); diff --git a/src/matlab/readamrdata_binary.m b/src/matlab/readamrdata_binary.m new file mode 100644 index 00000000..b21dc928 --- /dev/null +++ b/src/matlab/readamrdata_binary.m @@ -0,0 +1,180 @@ +function [amr,t] = readamrdata_binary(dim,Frame,dir,readblocknumber) + +% Internal routine for Clawpack graphics. + +n1 = Frame+10000; +fname = [dir, 'fort.',num2str(n1)]; +fname(length(dir)+6) = 't'; + +if ~exist(fname,'file') + amr = {}; + t = []; + fprintf('\n'); + fprintf('Frame %d (%s) does not exist ***\n',Frame,fname); + fprintf('\n'); + return +end + +% Read data from fname = 'fort.tXXXX' +fid = fopen(fname); + +t = fscanf(fid,'%g',1); fscanf(fid,'%s',1); +meqn = fscanf(fid,'%d',1); fscanf(fid,'%s',1); +ngrids = fscanf(fid,'%d',1); fscanf(fid,'%s',1); +fscanf(fid,'%d',1); fscanf(fid,'%s',1); % Read maux +ndim = fscanf(fid,'%d',1); fscanf(fid,'%s',1); +mbc = fscanf(fid,'%d',1); fscanf(fid,'%s',1); +fclose(fid); + +if ndim ~= dim + error('readamrdata_binary : Dimensions do not agree'); +end + +if ndim <= 2 + strip_ghost = true; +else + strip_ghost = true; +end + +% change the file name to read the binary data: +fname(length(dir) + 6) = 'q'; +if ~exist(fname,'file') + amr = {}; + t = []; + return; +end + +fid = fopen(fname); + +fname_bin = fname; +fname_bin(length(dir) + 6) = 'b'; + +fprintf('Reading data from %s ',fname_bin); + +if ~exist(fname_bin,'file') + amr = {}; + t = []; + return +end + + +if dim == 1 + amr_s = struct('gridno',[],'level',[],'blockno',[],'mx',[],'xlow',[],... + 'dx',[],'data',[]); +elseif dim == 2 + amr_s = struct('gridno',[],'level',[],'blockno',[],... + 'mx',[],'my',[],'xlow',[],'ylow',[],'dx',[],'dy',[],'data',[]); +else + amr_s = struct('gridno',[],'level',[],'blockno',[],... + 'mx',[],'my',[],'mz',[],'xlow',[],'ylow',[],'zlow',[],... + 'dx',[],'dy',[],'dz',[],'data',[],'q',[]); +end +amr = repmat(amr_s,ngrids,1); + +fid_bin = fopen(fname_bin); + +M = floor(ngrids/10); + +for ng = 1:ngrids + + % read parameters for this grid: + + amrdata.gridno = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + amrdata.level = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + if (dim > 1) + if (readblocknumber) + amrdata.blockno = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + else + amrdata.blockno = 1; + end + end + amrdata.mx = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + if (dim > 1) + amrdata.my = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + if (dim > 2) + amrdata.mz = fscanf(fid,'%d',1); fscanf(fid,'%s',1); + end + end + + amrdata.xlow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + if (dim > 1) + amrdata.ylow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + if (dim > 2) + amrdata.zlow = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + end + end + + amrdata.dx = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + if (dim > 1) + amrdata.dy = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + if (dim > 2) + amrdata.dz = fscanf(fid,'%g',1); fscanf(fid,'%s',1); + end + end + + mx = amrdata.mx; + if (dim > 1) + my = amrdata.my; + if (dim > 2) + mz = amrdata.mz; + end + end + + if strip_ghost + % Ghost cells written out + if (dim == 1) + ntot = (mx+2*mbc)*meqn; + q1d = fread(fid_bin,ntot,'double'); + q = reshape(q1d,meqn,ntot); + q1 = zeros(meqn,mx); + for m = 1:meqn + q1(m,:) = q(m,mbc+1:mx+mbc); + end + amrdata.data = q1; + elseif dim == 2 + ntot = (mx+2*mbc)*(my+2*mbc)*meqn; + q1d = fread(fid_bin,ntot,'double'); + q = reshape(q1d,meqn,ntot); + q1 = zeros(meqn,mx*my); + for m = 1:meqn + qt = reshape(q(m,:),(mx+2*mbc),(my+2*mbc)); + q1(m,:) = reshape(qt(mbc+1:mx+mbc,mbc+1:my+mbc),1,mx*my); + end + amrdata.data = q1; + elseif dim == 3 + ntot = (mx+2*mbc)*(my+2*mbc)*(mz+2*mbc)*meqn; + q1d = fread(fid_bin,ntot,'double'); + q = reshape(q1d,meqn,ntot); + q1 = zeros(meqn,mx*my*mz); + for m = 1:meqn + qt = reshape(q(m,:),(mx + 2*mbc),(my + 2*mbc),(mz + 2*mbc)); + q1(m,:) = reshape(qt(mbc+1:mx+mbc,mbc+1:my+mbc,mbc+1:mz+mbc),1,mx*my*mz); + end + amrdata.data = q1; + end + else + % Ghost cells not written out + if (dim == 1) + ntot = mx*meqn; + elseif dim == 2 + ntot = mx*my*meqn; + elseif (dim == 3) + ntot = mx*my*mz*meqn; + end + q1d = fread(fid_bin,ntot,'double'); + q = reshape(q1d,meqn,ntot); + amrdata.data = q; + end + + amr(ng) = amrdata; + + if mod(ng,M) == 0 + fprintf('.'); + end + +end +fprintf(' Done\n'); + +fclose(fid); + +end diff --git a/src/matlab/set_colors.m b/src/matlab/set_colors.m index e300aaf9..d53129cb 100644 --- a/src/matlab/set_colors.m +++ b/src/matlab/set_colors.m @@ -30,6 +30,12 @@ function set_colors(p,x,y,z,q,colormapping) cm_buff = 3; % Number of buffer under/over color entries uo = underover(); % User-defined function + if ~isfield(uo,'tol_lower') + uo.tol_lower = uo.tol; + end + if ~isfield(uo,'tol_upper') + uo.tol_upper = uo.tol; + end nmax = length(uo.colormap); cm_extended = [kron(ones(cm_buff,1),uo.color_under(:)'); ... @@ -42,7 +48,7 @@ function set_colors(p,x,y,z,q,colormapping) mfix = (uo.value_lower-uo.tol_lower) <= q & q <= uo.value_lower; q(mfix) = uo.value_lower; mfix = uo.value_upper <= q & q <= (uo.value_upper + uo.tol_upper); - q(mfix) = uo.value_upper-1e-8; % So floor works + q(mfix) = uo.value_upper-0.0*uo.tol_upper; % So floor works % ---------------------------------------------------- % Create index values in range [qlo,qhi] into @@ -50,7 +56,7 @@ function set_colors(p,x,y,z,q,colormapping) % ---------------------------------------------------- idx = 0*q + nan; % This will replace the 'cdata' property - % in the patch handle 'p'. + % in the patch handle 'p'. % Mask for values in the range [qlo, qhi]. m0 = uo.value_lower <= q & q <= uo.value_upper; @@ -60,13 +66,13 @@ function set_colors(p,x,y,z,q,colormapping) idx(m0) = cm_buff + floor(1 + slope*(nmax-1)); % Set under shoots to 1 and over shoots to nmax+2 - m_under = q <= uo.value_lower-uo.tol_lower; + m_under = q < (uo.value_lower-uo.tol_lower); idx(m_under) = 1; % first index in cm_extended - m_over = q >= (uo.value_upper + uo.tol_upper); + m_over = q > (uo.value_upper + uo.tol_upper); idx(m_over) = nmax + 2*cm_buff; % last index of cm_extended - m_leftover = ~(m0 | m_under | m_over); + m_leftover = ~m0 & ~m_under & ~m_over; idx(m_leftover) = -1; % ----------------------- @@ -87,7 +93,7 @@ function set_colors(p,x,y,z,q,colormapping) cm_extended = [cm_extended; color_nan]; % Set index to last entry in the colormap. - fv_idx(fv_idx < 0) = length(cm_extended); + fv_idx(fv_idx < 0) = length(cm_extended); % Hardwire colors for the patch set(p,'FaceVertexCData',cm_extended(fv_idx,:)); diff --git a/src/matlab/setopengl.m b/src/matlab/setopengl.m index e9de5a14..53841c35 100644 --- a/src/matlab/setopengl.m +++ b/src/matlab/setopengl.m @@ -1,4 +1,4 @@ -function setopengl() +function has_opengl = setopengl() % SETOPENGL sets the graphics renderer to OpenGL. % @@ -9,7 +9,10 @@ function setopengl() % it is suggested that it be used whenever possible. % % SETOPENGL sets the 'Renderer' property of the current figure to -% 'OpenGL'. +% 'opengl'. +% +% HAS_OPENGL = SETOPENGL() returns TRUE if the system has OpenGL and +% FALSE otherwise; % % Some graphics hardware may not support the OpenGL renderer, or may % cause certain Matlab routines to crash. For this reason, the choice @@ -19,14 +22,18 @@ function setopengl() rset = set(gcf,'Renderer'); found_opengl = 0; -for i = 1:length(rset), +for i = 1:length(rset) if (strcmp(lower(rset(i)),'opengl')) found_opengl = 1; end end if (~found_opengl) - disp('*** Warning : The OpenGL renderer is not available on your system.'); + disp('*** Warning : The OpenGL renderer is not available on your system.'); else - set(gcf,'Renderer','opengl'); -end; + if (nargout == 0) + set(gcf,'Renderer','opengl'); + else + has_opengl = found_opengl; + end +end diff --git a/src/matlab/sliceloop.m b/src/matlab/sliceloop.m index 618a1845..a32ee127 100644 --- a/src/matlab/sliceloop.m +++ b/src/matlab/sliceloop.m @@ -17,23 +17,23 @@ function sliceloop(sdir) % % While looping over frames, the user can enter 'x', 'y' or 'z' to loop % over slices while in the frame loop. -% +% % See also SHOWSLICES, HIDESLICES, SETPLOT3, PLOTCLAW3. if (nargin == 0) - sdir = input(['Enter x, y, or z to specify direction over which to',... - ' loop : '],'s'); -end; + sdir = input(['Enter x, y, or z to specify direction over which to',... + ' loop : '],'s'); +end setviews; slices = get_slices(sdir); -if (length(slices) == 0) - fprintf('sliceLoop : %sSliceCoords == []\n',sdir); - return; % Nothing to loop over -end; +if isempty(slices) + fprintf('sliceLoop : %sSliceCoords == []\n',sdir); + return; % Nothing to loop over +end % First hide all slices in direction dir. hideslices(sdir); @@ -42,29 +42,29 @@ function sliceloop(sdir) next_slice = 0; last_slice = 0; while (notdone) - s = input('SLICELOOP : Hit for next slice, or f to return to Frame loop : ','s'); - - if (isempty(s)) - next_slice = mod(next_slice+1,length(slices)+1); -% elseif (strcmp(s,'j')) -% next_slice = input('Input slice number : '); -% elseif (strcmp(s,'k')) -% keyboard; - elseif (strcmp(s,'f')) - return; - end; - hideslices(sdir,last_slice); - fprintf('\n'); - fprintf('Showing slice %d\n',next_slice); - showslices(sdir,next_slice); - - last_slice = next_slice; - - if exist('afterslice') == 2 - % make an m-file with this name for any other commands you - % want executed at the end of drawing each slice, for example - % to print a gif file for use in making an animation of a sweep-through - afterslice; - end - -end; + s = input('SLICELOOP : Hit for next slice, or f to return to Frame loop : ','s'); + + if (isempty(s)) + next_slice = mod(next_slice+1,length(slices)+1); + % elseif (strcmp(s,'j')) + % next_slice = input('Input slice number : '); + % elseif (strcmp(s,'k')) + % keyboard; + elseif (strcmp(s,'f')) + return; + end + hideslices(sdir,last_slice); + fprintf('\n'); + fprintf('Showing slice %d\n',next_slice); + showslices(sdir,next_slice); + + last_slice = next_slice; + + if exist('afterslice', 'file') + % make an m-file with this name for any other commands you + % want executed at the end of drawing each slice, for example + % to print a gif file for use in making an animation of a sweep-through + afterslice; + end + +end diff --git a/src/matlab/underover.m b/src/matlab/underover.m index 1e6f23dd..446d7706 100644 --- a/src/matlab/underover.m +++ b/src/matlab/underover.m @@ -57,4 +57,15 @@ % See also SETCOLORS, COLORBAR_UNDEROVER. % -error('underover : No user supplied ''underover''; see UNDEROVER'); +yrbcolormap; +cm = colormap; +uo = struct('color_under',[0 1 1],... % cyan + 'color_over',[1 0 1],... % magenta + 'color_nan',[1 1 1],... % white + 'value_lower',0, ... % theoretical minimum + 'value_upper',1,... % theoretical maximum + 'tol_lower', 1e-2,... + 'tol_upper',1e-2,... + 'colormap',cm); % Colors for non under/over shoot values + +