%%  Plots orthogonal, planar slice of FEATool result using iso lines/contours
%%      optionally adds color band contours. 3D plot of slice is also generated.
%%  Note: reqires prior export of fea result from FEAtool File menu.
%%  See:
%%      http://forum.featool.com/Slice-Plot-of-3D-Result-Showing-Iso-Contour-Lines-on-a-Slice-tp151.html
%%  Usage: isolineslice (PROB, evalexpr, orthoplane, val, lvl, cmap, addcolors)
%%      PROB:       fea, the variable exported by FEAtool
%%      evalexpr:	the mathematical expression to be evaluated and plotted
%%                  such as 's_dc*sqrt(Vx^2+Vy^2+Vz^2)', current density A/m-sq
%%      orthoplane:	'XY', 'YZ', or 'XZ' to denote orientation of slice plane
%%      val:        value of the invariant coordinate in the plane 
%%                      (e.g., for 'XY', z is invariant) e.g., 0.5
%%      lvls:       the number of iso levels to plot on the plane, e.g., 20
%%      cmap:       the color map to be used (for both 3D and 2D plots)
%%                  e.g., 'jet'
%%      addcolors:  (==1) Add colored gradient contour bands.
%%
%%      Example:	isolineslice(fea, 's_dc*sqrt(Vx^2+Vy^2+Vz^2)', 'XY', 0.0, 15, 'jet', 1 )
%%   Note:  Use single quotes,'sometext',(not double quotes "sometext")
%%          surrounding ascii/text/string/character data.
%%         

function isolineslice (PROB, evalexpr, orthoplane, val, lvls, cmap, addcolors)

% ISOLINESLICE isolineslice (PROB, evalexpr, orthoplane, val, lvls, cmap, isocolors)

    % Use fea exported from FEAtool file menu
    % Extract relevant data from fea structure
       
    % Check evalexpr input to be Voltage, E field (V/m) or current density
    % if not, set default.
    if isequal(evalexpr, 's_dc*sqrt(Vx^2+Vy^2+Vz^2)')
        quantity_and_units = 'Current Density (A/m)';
    elseif isequal (evalexpr, 'V' )
        quantity_and_units = 'Voltage (Volts)';
    elseif isequal(evalexpr, 'sqrt(Vx^2+Vy^2+Vz^2)')
        quantity_and_units = 'Voltage Gradient (Volts/meter), aka Electric Field';
    else   
       % set to the default
       evalexpr = 's_dc*sqrt(Vx^2+Vy^2+Vz^2)';      % set to current density
       quantity_and_units = 'Current Density (A/m)';
      end;
    
    titlestr2d = sprintf ('%s Iso Contours', quantity_and_units);
    
    % Check cmap input
    if ~isequal (cmap , 'jet') &&  ~isequal (cmap, 'cool') && ~isequal(cmap, 'spring') && ...
       ~isequal (cmap, 'summer')&&  ~isequal (cmap, 'autumn') && ~isequal (cmap, 'winter')  ...
       % set to the defautl, default
       cmap = 'parula';
    end;
    
    errmsg = sprintf (...
        "isolineslice: unrecognized orthoplane argument '%s',use only '%s', '%s', or '%s'", ...
            orthoplane, 'XY', 'XZ', 'YZ');      
    if (isequal (orthoplane, 'YZ'))
        xslice = val;
        yslice = [];
        zslice = [];
        xlabel_2D = 'y (m)';
        ylabel_2D = 'z (m)';
        orthoaxis = 'x';
    elseif (isequal (orthoplane, 'XZ'))
        yslice = val;
        xslice = [];
        zslice = [];
        xlabel_2D = 'x (m)';
        ylabel_2D = 'z (m)';
        orthoaxis = 'y';
    elseif (isequal (orthoplane, 'XY'))
        zslice = val;
        xslice = [];
        yslice = [];
        xlabel_2D = 'x (m)';
        ylabel_2D = 'y (m)';
        orthoaxis = 'z';
    else
    	error (errmsg);  
        return;
    end

    % Plot and extract 3D slice for specified orthoplane.
    % Do this uncondionally in case it was not done in featool.
    
    f = figure();
    h = postplot( PROB, 'parent', f, 'sliceexpr', evalexpr, ...
        'slicex', xslice, 'slicey', yslice, 'slicez', zslice, 'colormap', cmap );
    xlabel ('x (m)');
    ylabel ('y (m)');
    zlabel ('z (m)');

    titlestr3D = sprintf ('%s slice %s at %s coordinate %.2e', orthoplane, ...
                           quantity_and_units, orthoaxis, val);
    title (titlestr3D);

    h = findall(h,'type','patch');   % Get Matlab graphics handle to patch/slice.

    % Create 2D grid.
    grid.p = get(h,'Vertices')';
       
    if (isequal(orthoplane, 'YZ'))
        grid.p(1,:) = [];   % Delete x-coordinates.
    elseif (isequal(orthoplane, 'XZ'))
    	grid.p(2,:) = [];   % Delete y-coordinates.
    elseif (isequal(orthoplane, 'XY'))
    	grid.p(3,:) = [];   % Delete z-coordinates.
    else
    	error (errmsg);
    end
       
    grid.c = get(h,'Faces')';
    grid.a = gridadj(grid.c,2);
    grid.b = gridbdr(grid.p,grid.c,grid.a);
    grid.s = ones(1,size(grid.c,2));

    % Create help fea2 data struct.
    fea2.sdim = {'x', 'y'};
    fea2.grid = grid;
    fea2 = addphys(fea2,@poisson);
    
    
    fea2 = parsephys(fea2);
    fea2 = parseprob(fea2);
    fea2.sol.u = get(h,'FaceVertexCData');

    f2 = figure();
    
    % If want to display iso colors also:
    if (addcolors == 1) 
        postplot( fea2, 'parent', f2, 'surfexpr', 'u', 'isoexpr', 'u' );
    end
    
    postplot( PROB, 'parent', f2, 'boundary', 'on', 'surfexpr', 0, 'facecolor', 'none', 'colorbar', 'off', 'axis', 'off' );
    postplot( fea2, 'parent', f2, 'isoexpr', 'u', 'isomap', cmap, 'isolev',...
        lvls, 'title', 'Slice ISO Lines', 'xlabel', 'x', 'ylabel', 'y', 'zlabel', 'z', 'colorbar', 'off' );
    view(2);
    
    xlabel (xlabel_2D);
    ylabel (ylabel_2D);

    title (titlestr2d);

%   Now again not relying on gui?
       
%     Hmm, actually the whole "gui" field should probably not be present 
%     during fea struct export, so I would not recommend using any 
%     information from there. I think you can reconstruct the slice input 
%     from the actual slice coordinates, something like:
% 
%
%     
%     h = findall(0,'tag','fea_slice');   % Get patch directly from FEATool GUI.
% 
%    if( isempty(h) )   % MAYBE Plot slice manually if not available.
%        % but for now consider it an error/mis-use
%        error ("Must plot slice first using FEATool postprocesing settings");
%        % h = postplot( PROB, 'parent', f, 'sliceexpr', evalexpr, ...
%        % 'slicex', xslice, 'slicey', yslice, 'slicez', zslice );
%    end
% ...
%    p = fea.grid.p
%    TOL specifies a distance tolerance (default eps). UG-p759
%    TOL = sqrt(eps);
%    for i=1:3
%      if( all( abs(p(i,:) - p(i,1)) < TOL ) )
%        slice_pos_xyz{i} = mean(p(i,:));
%      else
%        slice_pos_xyz{i} = [];
%      end	 
%    end
% 
% Randal wrote:
% I would also like to add more control of the iso lines and the legend 
% such as the plotted range (begin and end values) and how to scale the 
% legend to fit on the page.
% 
% If the "isoval" input is an integer scalar then the contour lines are 
% computed as "linspace(min_val, max_val, isoval)", alternatively you 
% can input a vector of specific value yourself, 
% for example "isoval := [0.1, 0.2, 0.5]" or 
%     "isoval := linspace(0.1,0.5,10)".