function [ grid, stats ] = gridgen_scale( sin, varargin )
%@file gridgen_scale.m @brief GRIDGEN_SCALE Grid generation with scaling.
%
%   [ GRID, STATS ] = GRIDGEN_SCALE( SIN, VARARGIN ) Generates a
%   scaled grid/mesh for the geometry defined by the objects in SIN by
%   calling an external grid generation algorithm. SIN can be either a
%   valid fea problem struct with geometry defined in
%   SIN.geom.objects, a cell array of multiple geometry objects, or a
%   single geometry object. The geometry is first scaled according to
%   the input arguments, the scaled geometry is then meshed, and
%   finally the mesh the re-scaled to fit the original
%   geometry. Accepts the following property/value pairs are available
%   with the default grid generation algorithm.
%
%       Property    Value/{Default}           Description
%       -----------------------------------------------------------------------------------
%       axis        int {1}                   Scaling axis (1=x, 2=y, 3=z)
%       p_scal      scalar/array {midpoint}   Scaling points (along axis)
%       s_scal      scalar/array {[]}         Scaling factors (per scaling point)
%       varargin    pv pairs {defaults}       Property/value pairs to pass to gridgen
%
%   P_SCAL are coordinates along the AXIS to identify geometry objects
%   to scale. Geometry objects for which P_SCAL falls within the
%   bounding boxes will be scaled with the corresponding S_SCAL entry.
%
%   Examples:
%
%      1) 3x1 rectangle scaled by 1/3 in the x-direction.
%
%      grid = gridgen_scale( gobj_rectangle(0,3), 's_scal', 1/3 );
%      plotgrid( grid )
%
%      2) Polygon scaled by 1/5 in the y-direction.
%
%      grid = gridgen_scale( gobj_polygon([0 0;1 0;0.8 2;1.1 5; 0 5]), 'axis', 2, 's_scal', 1/5 );
%      plotgrid( grid )
%
%      3) Cylinder scaled by 1/6 in the z-direction.
%
%      grid = gridgen_scale( gobj_cylinder([0,0,0],0.5,5,3), 'axis', 3, 's_scal', 1/6 );
%      plotgrid( grid )
%
%      4) Scaling three parallel cylinders in a block (only scaling section with cylinders).
%
%      b1 = gobj_block(-0.75,3.25,-0.75,0.75,-0.5,5.5);
%      c1 = gobj_cylinder([0,0,0],0.5,5,3);
%      c2 = gobj_cylinder([1.25,0,0],0.5,5,3);
%      c3 = gobj_cylinder([2.5,0,0],0.5,5,3);
%      geom.objects = {b1, c1, c2, c3};
%      grid = gridgen_scale( geom, 'axis', 3, 's_scal', 1/5 );
%      plotgrid( grid, 'selcells', 'y>0' )
%      plotsubd( grid, 'linewidth', 1, 'labels', 'off' )
%      plotedg( grid, 'colors', {'k'}, 'labels', 'off', 'linewidth', 1 );
%      view(3)
%
%      4) Scale two blocks in a cylinder (only scaling sections with blocks).
%
%      s1 = gobj_sphere();
%      b1 = gobj_block(-1,1,-1,1,1.5,4.5);
%      b2 = gobj_block(-1,1,-1,1,-3.5,-1.5,'B2');
%      c1 = gobj_cylinder([0,0,-4],1.75,9,3);
%      geom.objects = {s1, b1, b2, c1};
%      grid = gridgen_scale( geom, 'axis', 3, 'p_scal', [-2.5, 3], 's_scal', [1/2, 1/4] );
%      plotgrid( grid, 'selcells', 'y>0' )
%      plotsubd( grid, 'linewidth', 1, 'labels', 'off' )
%      plotedg( grid, 'colors', {'k'}, 'labels', 'off', 'linewidth', 1 );
%      view(3)
%
%   See also GRIDGEN

% J.S. Hysing 210814.
% Copyright 2013-2021 Precise Simulation, Ltd.
if( ~(nargin || nargout) ),help gridgen_scale, return, end


% Assign input to geometry struct.
if( matlabVerLessThan([7,10,0]) )
  warning('off','MATLAB:divideByZero');
  warning('off','MATLAB:log:logOfZero');
end
if( iscell(sin) )
  geom.objects = sin;
elseif( isstruct(sin) )
  if( isfield(sin,'geom') )
    geom = sin.geom;
  elseif( isfield(sin,'objects') )
    geom = sin;
  else
    geom.objects = { sin };
  end
else
  geom = sin;
end
n_sdim = 0;
for i=1:length(geom.objects)
  n_sdim = max([n_sdim,geom.objects{i}.nsdim]);
end


% Parse input property/value pairs.
cOptDef = {'axis'     1;
           'p_scal'   [];
           's_scal',  [];
           'debug',   false;
           'fid',     1 };
if( nargin==0 )
  grid = cOptDef;
  stats = [];
  return;
end
[got,val] = parseopt( cOptDef, varargin{:} );


% Definition of scaling points (must be sorted in increasing order).
axis   = val.axis;     % Scaling axis
p_scal = val.p_scal;   % Scaling points (points to identify objects to scale).
s_scal = val.s_scal;   % Scaling factors.

if( isempty(p_scal) )   % Determine default scaling points.
  for j=1:length(geom.objects)
    o = geom.objects{j};
    p_scal(j) = 0.5*(o.bbox(2,axis) + o.bbox(1,axis));
  end
end
p_scal = unique(p_scal);
l_scal = inf(size(p_scal));              % Scaling box lengths.
b_scal = {};                             % Scaling bounding box.
i_scal = cell(1,length(geom.objects));   % Index to scalings to apply on each object.
m_scal = cell(1,length(geom.objects));   % Index to move object due to other scalings.
i_blk  = [];                             % Index to containing "blocks".


% Loop over scaling points/planes to determine scaling boxes (l_scal/b_scal).
if( isscalar(val.fid) )
  fprintf( val.fid, 'Determining scaling boxes...\n' );
end
for i=1:length(p_scal)

  % Loop over geometry objects.
  for j=1:length(geom.objects)
    o = geom.objects{j};
    pmin = o.bbox(1,axis);
    pmax = o.bbox(2,axis);

    if( pmin < p_scal(i) && p_scal(i) < pmax )
      % Geometry object contains scaling point.

      l = pmax - pmin;
      if( l < l_scal(i) )
        l_scal(i) = l;
        b_scal{i} = [pmin, pmax];
      end
    end
  end

end


% Loop over scaling boxes and determine geometry objects to scale (i_scal/m_scal).
if( isscalar(val.fid) )
  fprintf( val.fid, 'Determining scalings/transformations...\n' );
end
for i=1:length(b_scal)

  % Loop over geometry objects.
  for j=1:length(geom.objects)
    o = geom.objects{j};
    pmin = o.bbox(1,axis);
    pmax = o.bbox(2,axis);

    if( pmin >= b_scal{i}(1) && pmax <= b_scal{i}(2) )
      % Geometry object contained within scaling box
      %   => apply scaling

      i_scal{j} = [i_scal{j}, i];

    elseif( b_scal{i}(2) <= pmin )
      % Geometry object "above" scaling box
      %   => move object "down"

      m_scal{j} = [m_scal{j}, i];

    elseif( b_scal{i}(1) > pmin  && pmax > b_scal{i}(2) )
      % Scaling box enrirely contained within geometry object
      %   => partially scale surrounding object

      i_blk = [i_blk, j];
      i_scal{j} = [i_scal{j}, i];

    end
  end
end


% Loop over geometry objects to apply scalings/transformations.
if( isscalar(val.fid) )
  fprintf( val.fid, 'Applying scalings/transformations...\n' );
end
for j=1:length(geom.objects)
  o = geom.objects{j};
  if( val.debug )
    fprintf('Processing object %d:%s', j, o.tag );
  end
  pmin = o.bbox(1,axis);
  pmax = o.bbox(2,axis);

  offset  = [0,0,0];   % Offset to move in axis direction.
  scaling = [1,1,1];   % Factor to scale object with.

  if( ~isempty(i_scal{j}) && ~any(j==i_blk) )   % Scale (and move back) object.
    is = i_scal{j};

    if( length(is)==1 )
      s = s_scal(is);
      offset(axis) = offset(axis) + (1-s)*o.bbox(1,axis);
      scaling(axis) = s;
      if( val.debug )
        fprintf(' scaling(%d) with s=%f', is, s );
      end
    end

    % else   % Move if object contained within scaling box.

    %   for i=1:length(b_scal)
    %     if( o.bbox(1,axis) >= b_scal{i}(1) && o.bbox(2,axis) <= b_scal{i}(2) )
    %       offset(axis) = offset(axis) + (1-s_scal(i))*(o.bbox(1,axis) - b_scal{i}(1));
    %       m_scal{j} = setdiff(m_scal{j},i);
    %     end
    %   end
  end

  % Translation due to "lower" scaling.
  if( ~isempty(m_scal{j}) )
    for im=m_scal{j}
      s = s_scal(im);
      l = l_scal(im);
      offset(axis) = offset(axis) - (1-s)*l;
    end
    if( val.debug )
      fprintf(' moving %f', offset(axis) );
    end
  end

  % Apply transformation to geometry object.
  if( any(offset) || ~all(scaling==1) )
    geom_j.objects = {o};
    geom_j = geom_apply_transformation( geom_j, geom_j.objects{1}.tag, [0,0,0], scaling(1:o.nsdim) );
    geom_j = geom_apply_transformation( geom_j, geom_j.objects{1}.tag, offset(1:o.nsdim) );
    geom.objects{j} = geom_j.objects{1};
    geom.objects{j}.tag = o.tag;
  end
  if( val.debug )
    fprintf('\n');
  end
end


if( isscalar(val.fid) )
  fprintf( val.fid, 'Scaling containing objects...\n' );
end
i_blk = unique(i_blk);
for j=1:length(i_blk);
  o = geom.objects{i_blk(j)};
  offset  = [0,0,0];   % Offset to move in axis direction.
  scaling = [1,1,1];   % Factor to scale object with.
  l0 = 0; l1 = 0;
  for i=1:length(i_scal{j})
    l0 = l0 + (b_scal{i_scal{j}(i)}(2) - b_scal{i_scal{j}(i)}(1));
    l1 = l1 + s_scal(i_scal{j}(i))*(b_scal{i_scal{j}(i)}(2) - b_scal{i_scal{j}(i)}(1));
  end

  s = ((o.bbox(2,axis) - o.bbox(1,axis)) - (l0-l1)) / ...
      (o.bbox(2,axis) - o.bbox(1,axis));
  offset(axis) = offset(axis) + (1-s)*o.bbox(1,axis);
  scaling(axis) = s;

  geom_j.objects = {o};
  geom_j = geom_apply_transformation( geom_j, geom_j.objects{1}.tag, [0,0,0], scaling(1:o.nsdim) );
  geom_j = geom_apply_transformation( geom_j, geom_j.objects{1}.tag, offset(1:o.nsdim) );

  geom.objects{i_blk(j)} = geom_j.objects{1};
  geom.objects{i_blk(j)}.tag = o.tag;
end


% Generate grid for scaled geometry.
[grid, stats] = gridgen( geom, varargin{:} );


% Stretch grid to match original geometry.
if( isscalar(val.fid) )
  fprintf( val.fid, 'Re-scaling grid...\n' );
end
for i=length(p_scal):-1:1
  for j=1:length(i_scal)
    if( ~any(j==i_blk) && isequal(i_scal{j},i) )
      b2_scal = geom.objects{j}.bbox(:,axis);   % Scaled bounding box.
      break;
    end
  end
  % clf, plot3(grid.p(1,:), grid.p(2,:), grid.p(3,:), 'k.'), hold on, view(3)

  px = grid.p(axis,:);

  % Add offset for partial scaling.
  ix = b2_scal(1) <= px & px <= b2_scal(2);
  % plot3(grid.p(1,ix), grid.p(2,ix), grid.p(3,ix), 'ro'), pause
  grid.p(axis,ix) = px(ix) + (1/s_scal(i))*(px(ix) - b2_scal(1));
  % plot3(grid.p(1,ix), grid.p(2,ix), grid.p(3,ix), 'yo'), pause

  % Add offset for completely scaled boxes.
  ix = px > b2_scal(2);
  % plot3(grid.p(1,ix), grid.p(2,ix), grid.p(3,ix), 'bv'), pause
  grid.p(axis,ix) = grid.p(axis,ix) + (1/s_scal(i))*(b2_scal(2)-b2_scal(1));
  % plot3(grid.p(1,ix), grid.p(2,ix), grid.p(3,ix), 'gs'), pause

end


assert( gridcheck(grid) == 0, 'Grid generation with scaling returned invalid grid.' )


% f = figure;
% plotgrid( grid, 'selcells', 'y>0.03', 'parent', f )
% plotsubd( grid, 'linewidth', 1, 'labels', 'off', 'parent', f )
% plotedg( grid, 'colors', {'k'}, 'labels', 'off', 'linewidth', 1, 'parent', f );
% view(3)
% drawnow
