function [PROB,GEO] = buildAndSliceTest
%% Build a model consisting of concentric spheres centered in a block
%% that has been sliced to add a boundary perpendicular to x-axis and
%% in the center.
%%

% Prepare new 3D fea object for function return
fea.sdim = { 'x', 'y', 'z' };
fea = addphys( fea, @conductivemediadc, { 'V' } );

%   Create empty geometry object
g = struct;

%   Create a figure on which to plot
f = figure();

%	Create Block, B1,  containing a spheres, S1, S2 
%   then add a cut plane in the middle perpendicular to x-axis 
%   all centered at (.5, .5, .5)

%   Create B1
gobj = gobj_block( [0], [1], [0], [1], [0], [1], 'B1' );
g = geom_add_gobj( g, gobj );

disp ("...paused...<Enter> to continue");
update_plot ("Create B1... (<Enter> to continue)");
    
%   Create S1
gobj = gobj_sphere( [.5 .5 .5], [.2], [1], 'S1' );
g = geom_add_gobj( g, gobj );
update_plot ("Create S1 at the center of B1");

%   Create S2 outside B1
gobj = gobj_sphere( [.15 .5 .5], [.1], [1], 'S2' );
g = geom_add_gobj( g, gobj );
update_plot ("Create small S2 outside S1 near x=0 face of B1");

%   Duplicate B1, S1, and S2 in-place, creating B2, S3, S4
duplicate_objects ({'B1', 'S1', 'S2'});
update_plot ("B1, S1, and S2 all duplicated in place ==> B2,S3,S4");

%   Create enclosing blocks that share a face, an X-perpendicular plane, 
%   at x=0.5
gobj = gobj_block( [.5], [1.5], [-.5], [1.5], [-.5], [1.5], 'B3' );
g = geom_add_gobj( g, gobj );
gobj = gobj_block( [-.5], [.5], [-.5], [1.5], [-.5], [1.5], 'B4' );
g = geom_add_gobj( g, gobj );
update_plot (...
    "B3 & B4 created with a common face (x=.5) enclosing all other objects.");

%    Subtract B3 & B4 from B1, B2 respectively
g = geom_apply_formula( g, 'B1-B3' );
g = geom_apply_formula( g, 'B2-B4' );
update_plot (...
    "B1 'sliced' by subtracting Blocks B3 from it and B4 from B2, its copy");

%   Re-create enclosing blocks that share X-perpendicular plane, x=0.5
gobj = gobj_block( [.5], [1.5], [-.5], [1.5], [-.5], [1.5], 'B5' );
g = geom_add_gobj( g, gobj );
gobj = gobj_block( [-.5], [.5], [-.5], [1.5], [-.5], [1.5], 'B6' );
g = geom_add_gobj( g, gobj );
update_plot ("Enclosing Blocks, again created ...(B5 and B6)");

%    Subtract B5 & B6 from S1 and  its dup S3 respectively
g = geom_apply_formula( g, 'S1-B5' );
update_plot("after 'S1-B5' but before complementary 'S3-B6'");
g = geom_apply_formula( g, 'S3-B6' );
update_plot("S1 'slicing' completed by 'S3-B6'");

%   Again, re-create enclosing blocks
gobj = gobj_block( [.5], [1.5], [-.5], [1.5], [-.5], [1.5], 'B7' );
g = geom_add_gobj( g, gobj );
gobj = gobj_block( [-.5], [.5], [-.5], [1.5], [-.5], [1.5], 'B8' );
g = geom_add_gobj( g, gobj );
update_plot ("Enclosing Blocks created a third time...(B7 & B8)");

%    Subtract B7 & B8 from S2 and its dup S4 respectively
g = geom_apply_formula( g, 'S2-B7' );
update_plot(".... after 'S2-B7'... no change");
disp ("...about to call 'geom_apply_formula' for non overlapping objects.");
disp ("Hit <enter> to continue...");
%pause;
g = geom_apply_formula( g, 'S4-B8' );
update_plot(...
    "Finally, after 'S4-B8' (crash of GEOMTool) ... again, no change");


%	returm geometry for import by FEATool GUI
set_return_vars;

    function duplicate_objects (gobj)
        for o = gobj
            gobj_tag = {o};
            [g, gobj_tag] = copy_geometry_object( gobj_tag, g );
        end
    end

    function set_return_vars
        fea.geom = g;      
        PROB = fea;
        GEO  = g;
    end

    function update_plot (plot_title)   
        clf (f);
        title (plot_title);
        plotgeom (g,'parent', f);
        view(3);
        %pause;
    end 



end
