
function [GOBJ_SLICED, stat] = ...
    geom_slice_model ( GOBJ, orthoplane, val, plot )
% [GOBJ_SLICED, stat] = geom_slice_model ( GOBJ, orthoplane, val, plot )
%
%   slices entire model GOBJ with orthogonal (axis-perpendicular) plane
%
% Slices an FEATool GOBJ with a [cut-]plane effectively dividing any object
% that it interects at the plane of intersection, thus creating two
% distinct components which share an internal boundariy at the
% intersection. These boundaries my be integrated over in the 
% Postprocessing mode. The component pieces generated by the slice may be
% subsequently deleted, transformed, or combined with other component 
% objects to further develop the model geometry in the usual manner.
%
% The first argument is a geometry object such as fea.geom, orthoplane and
% val define the cutting plane, and plot enables the plotting of figures
% including the unmodified model that is input (GEOB) and the output
% sliced/cut model geometry (GEOB_SLICED).  The final plot adds a
% highlighed and enlarged component to illustrate the position of the 
% underlying slice/cut.
%
% Inputs:
%	GOBJ         :  a geometric object struct as in the fea struct member 
%                   geom(fea.geom). 
%	orthoplane   :	'XY', 'YZ', or 'XZ' to denote orientation of the 
%                       slice plane
%	val          :	The value of the invariant coordinate in the plane
%                       (e.g., for 'XY', z is invariant) e.g., 0.5
% 	plot         :	[true] If omitted or true, the geometry model is 
%                   plotted at various points, including before and after 
%                   the slice is done.
%	Example      :	geom_slice_model(fea.geom, 'XY', 0.5, true)
%                   geom_slice_model(mygeom, 'XY', 0.5)
%
% Outputs:
%   GOBJ_SLICED : the resulting sliced geometry model.
%
%   stat        :
%    0 : normal return. GOBJ_SLICED output variable is set to result.
%   -1 : wrong number of arguments.
%   -2 : invalid slice plane designation. (orthoplane)
%   -3 : invalid slice plane location. val places it outside the model GOBJ
%
% Note: Each time a plot is displayed, there is a delay of about 2 seconds 
% in the processing to allow for time to view the plot. This is governed by
% a variable, dwell, that is passed as an argument to the "pause" function.
% It is set near the top of the code and can be changed to a more suitable
% value (including 0).
%
% See also      geom_slice_model_gui   geom_bounds_and_size

% Set true for development
debug_log = false;

% initialize return variables.
stat = 0;
GOBJ_SLICED = struct;

% number of seconds to pause when displaying plots
% (possible future option)
dwell = 2;

if( nargin == 3 )
    do_plot = true;
else
    if (nargin == 4)
        do_plot = plot;
    else
        % indicate bad argument list
        disp ("**********");
        disp(...
       "Err: geom_slice_model, incorrect number of arguments: " + nargin);
        disp ("**********");
        disp ("");
        disp ("See 'help geom_slice_model'");
        stat = -1;
        return;
    end  
end

if debug_log
    % override input plot request
    do_plot = true;
end

if do_plot
    % Create or obtain figure handle f to plot progress during
    f = gcf; %figure();
end

% Determine the size and location of a block that will completely enclose
% the entire geometry - all its objects/components using an external
% function.
[p_min, p_max, size] = geom_bounds_and_size (GOBJ);

% print_bounds_and_size;
update_plot ("Input Model Before Slice");

% The block defined by p_min and p_max is considered to have two halves, 
% blockA and blockB, so that they share a face that is in the slicing 
% plane. blockA and blockB are created and defined below.

% set x, y, and z position variable indexes for convinenence
% and readability, e.g., p_min(x,y,z) would be p_min(
x = 1;  %   p(1) is the p x coordinate
y = 2;  %   p(2) is the p y coordinate
z = 3;  %   p(3) is the p z coordinate

% determine plane (XY, YZ, or XZ) and its invariant value, (x, y, or z)
% e.g., for the 'XY' plane the invariant value is a z coordinate.
slice_plane = orthoplane;
slice_location = val;

% Process input text argument, slice_plane.
% If called by wapper, slice_plane will have already been checked,
if (isequal (slice_plane, 'XY'))
    slice_axis_index	= z;
    slice_axis_name     ='z';
elseif (isequal (slice_plane, 'YZ'))
    slice_axis_index    = x;
    slice_axis_name     ='x';
elseif (isequal (slice_plane, 'XZ'))
    slice_axis_index    = y;
    slice_axis_name     ='y';
else
    % do not continue with invalid slice_plane
    err_msg = sprintf("geom_slice_model: Invalid slice plane, %s", ...
        slice_plane);
    disp (err_msg);
    
    stat = -2;
    return;
end

% check that val is between its axis min and max (before adding 
% the arbitraty 20% pad)

% Note: this cannot be checked in the caller (wrapper) unless p_min/p_max
% are calculated there too... or passed in ... nah.
if ~(val > p_min (slice_axis_index) & ...
     val < p_max (slice_axis_index))
 
    % the proposed cut plane does not intersect the model geometry.
    % at most it shares a face/line/point.  It cold not divide the 
    % model into two pieces.
    
    % do not continue with this invalid val
    err_msg = sprintf( ...
    "Invalid %s slice plane - it does not intersect model at %c = %f", ...
        slice_plane, slice_axis_name, val);
    disp (err_msg);
  
    plot_slice_plane (GOBJ, 'InvalidSlicePlane', "Err: " + err_msg);
    
    error_dlg = errordlg(err_msg, 'Input Invalid');
    waitfor(error_dlg);

    stat = -3;
    return;
end
    
% plot_slice_plane (GOBJ, 'Review selected plane', "Proposed Cut Plane");
    
% Add arbitrary 20 percent to the p_min/p_max enclosing block dimensions
% by reducing p_min by 10% of size and adding 10% of size to p_max.
p_min = p_min - 0.1*size;
p_max = p_max + 0.1*size;

% for consistency, update size also.
size = 1.2*size;

% print_bounds_and_size;

% Create the bounding points (diagonal vertex coordinates) for blockA and
% blockB - haveing a common face being the slice/cut plane.

% blockA's minimum point vertex, blockA_min, is the same as that of the 
% larger enclosing block defined by p_min and p_max.
blockA_min = p_min;

% The coordinates of blockA_max, the diagonal vertex, will be the same as 
% p_max for two of is dimensions and the other dimension is set by the 
% desired slicing plane location.
blockA_max = p_max;
blockA_max(slice_axis_index) = slice_location;


% In a similar but converse manner for half-block blockB:
blockB_max = p_max;
blockB_min = p_min; 
blockB_min(slice_axis_index) = slice_location;

%=========================
% Split the p_min/p_max block by a slice_plane locted at the slice_location
% and define parameters of both half blocks that share that bisecting plane
% as a face.
%=========================

% Create the two block halves of the enclosing block (two half-blocks), 
% blockA and blockB, from the diagonal point pairs created above.

% blockA
blockA = gobj_block( ...
        [blockA_min(x)], [blockA_max(x)], ...
        [blockA_min(y)], [blockA_max(y)], ...
        [blockA_min(z)], [blockA_max(z)], ...
        'SLICEblockA1' );  % tentative tag id, see below

% blockB
blockB = gobj_block( ...
        [blockB_min(x)], [blockB_max(x)], ...
        [blockB_min(y)], [blockB_max(y)], ...
        [blockB_min(z)], [blockB_max(z)], ...
        'SLICEblockB1' );  % tentative tag id, see below
    
% for each GOBJ object (GOBJ.objects), obj
%   duplicate obj in place
%   locate tag or index for duplicate object (dup_obj_tag)
%   add half blockA to GOBJ
%   subtract blockA from object (obj_tag)
%   add half blockB to GOBJ
%   subtract blockB from duplicate object (dup_obj_tag)
%   Note:   if a subtract should fail, revert to the GOBJ state before that
%           iteration (of i) and proceed with the next iteration.

% first, get the number and list of tags for all component objects in model
% to be sliced.
num_model_objs = length (GOBJ.objects);
for i = 1:1:num_model_objs
    model_obj_tags{i} = GOBJ.objects{i}.tag;
end

% initialized blockA and blockB tag numeric suffix to 0..

for i = 1:1:num_model_objs
    % Save the GOBJ state in case a geom_apply_formula error occurs.
    % It is saved in one place (here, at the top of the dup/subt loop.)
    % It is restored to prev_geo if either subtract fails and the loop 
    % restarts at the next iteration of i.

    prev_geo = GOBJ;
  
    % get tag for model object being sliced in this pass from the static
    % list made before the dup/subtract loop began.
    obj_tag = model_obj_tags{i};  
    
    % Duplicate the obj
    [GOBJ, dup_obj_tag] = copy_geometry_object( obj_tag, GOBJ );
    
    % blockA.tag is set above to SLICEblockA1, tentatively. The actual tag 
    % assigned to the added object will be selcted by geom_add_gobj so as 
    % to be unique. No attempt is made to create a unique tag beforehand.
    % After the block is added, its tag is read from the last objects{}.tag
    % field and the blockA is replaced with the actual tag before the 
    % subtract operation.
    
    GOBJ = geom_add_gobj( GOBJ, blockA );
    blockA.tag = GOBJ.objects{end}.tag;
            
    % update_plot("after add blockA.tag:"+blockA.tag);

    % Now proceed with blockA subtraction from the ith model object
    % (the one with the tag: obj_tag).
    
    formula = sprintf ('%s - %s', obj_tag, blockA.tag );
    if debug_log
        fprintf ("slicing object, i=%i, obj_tag=%s\n", i, obj_tag);
    end
    [GOBJ, gstat ] = geom_apply_formula( GOBJ, formula, 1e-5, false);

    % gstat output variable is the number of failed geometry object
    % combinations.
   
    if gstat ~= 0
        % if zero (0) is not returned from geom_apply_formula, IT IS 
        % ASSUMED that the operation could not be done because the object
        % being sliced was completely contained by or outside of the 
        % SLICEblock[A|B]n.  If a test is ever developed for this (complete
        % intersection or no intersection), then this test should be added 
        % before the operation is done and the object skipped to avoid
        % having to reverse course upon a failed subtraction.
       
        % The recovery from this is to restore the geometry object, GOBJ, 
        % to its state before this iteration of 'for i=' and resume with 
        % the next iteration.
        
        if debug_log
            disp("WARN: geom_apply_formula for blockA, return gstat: "+...
                gstat+", for: '"+formula+"'");
        end
        % abandon this 'i' iteration restoring GOBJ to its previous, 
        % top of the loop value.
        GOBJ = prev_geo;     
    else
        % Normal return from geom_apply_formula:
                
        % Subtract blockB from duplicate object with dup_obj_tag.

        % Note: As with blockA (see above) the actual blockB tag is set
        % after the geom_add_gobj operation.
        GOBJ = geom_add_gobj( GOBJ, blockB );
        blockB.tag = GOBJ.objects{end}.tag;
        
        formula = sprintf ('%s - %s', dup_obj_tag{1}, blockB.tag );
        if debug_log
            fprintf ("slicing object, i=%i, dup_obj_tag=%s\n", ...
                i, dup_obj_tag{1});
        end
        
        % gstat output variable is the number of failed geometry object
        % combinations.
        [GOBJ, gstat ] = geom_apply_formula( GOBJ, formula, 1e-5, false );  

        if gstat ~= 0
            % See comment beneath enclosing 'if (~isequal(gstat , 0))'
            if debug_log
                disp( ...
                "WARN: geom_apply_formula for blockB return gstat : "+ ...
                gstat +", for: '"+formula+"'");
            end
            
            % abandon this 'i' iteration restoring GOBJ to its previous, 
            % top of the loop value.
            GOBJ = prev_geo;     
        end    
    end  
end

plot_title = sprintf ( ...
    "Model geometry sliced by %s cut plane at %c = %f", ...
    orthoplane, slice_axis_name, val);
update_plot (plot_title);

% stat should be ==0 here, but if not, just return and let the caller
% check for unanticipated error.

% returm the sliced geometry
GOBJ_SLICED  = GOBJ;

% show highlighted slice plane (does not modify GOBJ)
plot_slice_plane(GOBJ, '      SlicePlane', plot_title);

return;

    function print_object_tags (g)
        for j=1:1:length(g.objects)
            disp ("  "+g.objects{j}.tag);
        end
    end


    function print_bounds_and_size()
        fprintf ( ...
           'geometry obj bounds and size\n       x         y         z\n');
        fprintf ('p_min: %10f  %10f  %10f\n', p_min);
        fprintf ('p_max: %10f  %10f  %10f\n', p_max);
        fprintf (' size: %10f  %10f  %10f\n', size);
    end


    function update_plot (plot_title)
        if do_plot
            clf (f);
            title (plot_title);
            plotgeom (GOBJ,'parent', f, 'fontsize', 12);
            view(3);
            pause(dwell);
        end
    end 

    function plot_slice_plane (geo_obj, slice_plane_tag, plot_title)
        % This function uses global variables 
        % Copy geo_obj so that the following modifications are not made to 
        % global geo_obj.
        
        geo_cpy = geo_obj;
        
        % Nudge the value for p1(slice_axis_index) down a bit, and
        % nudge the value for p2(slice_axis_index) up a bit.
        % Then create a block whose diagonal vertexes are p1 and p2.
        % This creates a thin (based on size of "a bit") block that
        % illustrates the location of the slice plane, which is located in
        % the center 
        
        % controls the thin-ness of the plane illustrative block.
        a_bit = .0001 * max(size);
        
        % make the illustrative plane block extend 20% further than its
        % intersection with the model.
        p1 = p_min - 0.2*size;
        p2 = p_max + 0.2*size;
        
        % corners of pseudo slice block to plot
        p1 (slice_axis_index) = val - a_bit;
        p2 (slice_axis_index) = val + a_bit;
        
        % For readability, create appropriately named variables for the
        % each member of p1 and p2:
        [xmin, ymin, zmin] = deal (p1(1), p1(2), p1(3));
        [xmax, ymax, zmax] = deal (p2(1), p2(2), p2(3));
        
        % create the very thin block to represent the slice
        block_obj = gobj_block( ...
            [xmin], [xmax], ...
            [ymin], [ymax], ...
            [zmin], [zmax], ...
            slice_plane_tag, false);
        
        % add block-plane to copy of sliced geometry and 
        geo_cpy = geom_add_gobj(geo_cpy, block_obj);
        if do_plot
            clf (f);
            title (plot_title);
            % highlight the last component added to GOBJ (the pseudo slice)
            plotgeom (geo_cpy,'parent', f, 'fontsize', 12, ...
                'highlight', length(geo_cpy.objects));
            view(3);
            pause(dwell);
        end
    end
end
    



