%%  Plots orthoganol planar slice of FEAtool result using iso lines/contours
%%  Note: reqires 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, orthoplane, val, lvl)
%%      PROB:       the fea variable exported by FEAtool
%%      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
%%      lvls:       the number of iso levels to plot on the plane.

function isolineslice (PROB, orthoplane, val, lvls)
% ISOLINESLICE isolineslice (PROB, orthoplane, val, lvls)
     
    % Use fea exported from FEAtool file menu
    % Extract relevant data from fea structure
       
    expression  = PROB.gui.post.sliceexpr;
    cmap        = PROB.gui.post.colormap;
    
    %Get slice directives from gui post struct when updated
    %xslice      = PROB.gui.post.slicex
    %yslice      = PROB.gui.post.slicey
    %zslice      = PROB.gui.post.slicez
    
    errmsg = sprintf (...
        "isolineslice: unrecognized orthoplane argument '%s',use only '%s', '%s', or '%s'", ...
            orthoplane, "XY", "XZ", "YZ");      
    if (orthoplane == 'YZ')
        xslice = val;
        yslice = [];
        zslice = [];
    elseif (orthoplane == 'XZ')
        yslice = val;
        xslice = [];
        zslice = [];
    elseif (orthoplane == 'XY')
        zslice = val;
        xslice = [];
        yslice = [];
    else
    	error (errmsg);     
    end
       
    % Plot and extract 3D slice.
    f = figure();
    h = postplot( PROB, 'parent', f, 'sliceexpr', expression, ...
        'slicex', xslice, 'slicey', yslice, 'slicez', zslice );
    h = findall(h,'type','patch');   % Get Matlab graphics handle to patch/slice.

    % Create 2D grid.
    grid.p = get(h,'Vertices')';
        
    if (orthoplane == 'YZ')
        grid.p(1,:) = [];   % Delete x-coordinates.
    elseif (orthoplane == 'XZ')
    	grid.p(2,:) = [];   % Delete y-coordinates.
    elseif (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:
    % postplot( fea2, 'parent', f2, 'surfexpr', 'u', 'isoexpr', 'u' );
    postplot( fea2, 'parent', f2, 'isoexpr', 'u', 'isomap', cmap, 'isolev', lvls );
    