function [ thickness, fea ] = fzero_featool_example()
% FZERO_FEATOOL_EXAMPLE Example for plane stress for hole in plate.
%
%   Example m-script using the built-in MATLAB function FZERO to find
%   for which thickness the maximum stress in the x-direction for the
%   FEATool hole in plate plane stress model is equal to 3e7.
%
%   See also fzero, ex_planestress1.

% Copyright 2013-2019 Precise Simulation, Ltd.


sx_target = 3e7;
thickness_range = [1e-4 1e-2];   % Initial search range.

% Problem function definition.
pfun = @(thickness) sx_target - compute_max_stressx(thickness);

[thickness,diff,exitflag,output] = ...
    fzero( pfun, thickness_range, optimset('Display','iter','TolX',1e-6) );

fprintf( '\n   Thickness found = %g for which the max target x-stress differ by %4.1g%%.\n\n', thickness, abs(diff/sx_target)*100 );

if( ~nargout )
  clear thickness
elseif( nargout>=2 )
  fea = get_fea_struct( thickness );
end


%------------------------------------------------------------------------------%
function [ sx_max ] = compute_max_stressx( thickness )

fea = get_fea_struct( thickness );   % Finite element fea problem struct.

fea.sol.u = solvestat( fea, 'fid', [] );   % Call solver without terminal output.

% Define an expression for the stress in the x-direction.
E11  = 'E_pss/(1-nu_pss^2)';
E12  = ['nu_pss*',E11];
s_sx = [E11,'*ux + ',E12,'*vy'];

[~,sx_max] = minmaxsubd( s_sx, fea );   % Calculate the maximum stress in the x-direction.

%------------------------------------------------------------------------------%
function [ fea_out ] = get_fea_struct( thickness )

persistent fea

if( isempty(fea) )   % Setup new fea struct (from the ex_planestress1 example).

  force = 1000;
  nu    = 0.3;
  E     = 210e9;

  % Geometry and grid.
  fea.sdim = { 'x' 'y' };
  gobj1 = gobj_rectangle( 0, 0.05, 0, 0.05, 'R1' );
  gobj2 = gobj_circle( [0, 0], 0.005, 'C1' );
  fea.geom.objects = { gobj1 gobj2 };
  fea = geom_apply_formula( fea, 'R1-C1' );
  try
    fea.grid = gridgen( fea, 'hmax', 0.001, 'meshgen', 'gridgen2d', 'fid', [] );
  catch
    fea.grid = gridgen( fea, 'hmax', 0.001 );
  end

  % Physics mode and equation coefficients.
  fea = addphys(fea,@planestress);
  fea.phys.pss.eqn.coef{1,end} = { nu };
  fea.phys.pss.eqn.coef{2,end} = { E  };
  fea.phys.pss.sfun            = { 'sflag1', 'sflag1' };

  % Boundary conditions.
  n_bdr    = max(fea.grid.b(3,:));
  bctype   = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
  bc_left  = findbdr( fea, ['x<=',num2str(0.001)] );
  bc_right = findbdr( fea, ['x>=',num2str(0.05-0.001)] );
  bc_lower = findbdr( fea, ['y<=',num2str(0.001)] );
  bctype{1,bc_left}  = 1;
  bctype{2,bc_lower} = 1;
  fea.phys.pss.bdr.coef{1,5} = bctype;
  bccoef = mat2cell( zeros(2,n_bdr), [1 1], ones(1,n_bdr) );
  bccoef{1,bc_right} = 'force/area';
  fea.phys.pss.bdr.coef{1,end} = bccoef;

  fea.expr = { 'force', force ;
               'area', '2*0.05*thickness' ;
               'thickness', thickness };

  fea = parsephys( fea );
  fea = parseprob( fea );

else  % Update thickness in stored fea struct.

  fea.expr{3,end} = thickness;

end

fea = parseprob(fea);
fea_out = fea;
