% SPANNER_PARAMETRIC Parametric stress and displacement calculation of a spanner.
%
%   SPANNER_PARAMETRIC Example to import a Gmsh grid and calculate displacements
%   on a fixed spanner. The load force is distributed in the tangential load
%   direction through the force fraction parameter FRAC.

% Copyright 2013-2017 Precise Simulation, Ltd.


frac = linspace(0,1,200);   % Fraction of lateral to inline force distribution 0..1

fs = 8;   % Font size.
lw = 2;   % Line width.

% Define scalings factors.
DSCALE = 2;   % Deformation plot scaling.
USE_METERS = true;
s = double( ~USE_METERS + USE_METERS*1e-3 );   % m to mm conversion factor.


% Import (and scale) grid.
fea.sdim = {'x','y','z'};
fea.grid = impexp_gmsh( 'spanner.msh', 'import' );
fea.grid.p = fea.grid.p*s;


% Add 3D solid linear elasticity physics mode and define material parameters.
fea  = addphys( fea, @linearelasticity );
sfun = 'sflag1';
fea.phys.el.sfun            = { sfun sfun sfun };
fea.phys.el.eqn.coef{1,end} = { 0.29 };
fea.phys.el.eqn.coef{2,end} = { 190e3/s^2 };


% Set all boundaries to no load per default.
n_bdr  = max(fea.grid.b(3,:));
bc_sel = cell(3,n_bdr);
[bc_sel{:}] = deal(0);


% Fix all displacements on boundaries 11 and 14.
i_fix = [11 14];
[bc_sel{:,i_fix}] = deal(1);
fea.phys.el.bdr.coef{5} = bc_sel;


% Apply force for x > 140 mm.
i_force = [ 1 4 ];
force = 1000/(6*s*80*s);
fea.phys.el.bdr.coef{7}{1,i_force(1)} = ['-',num2str((1-frac(1))*force),'*(y>140*',num2str(s),')'];
fea.phys.el.bdr.coef{7}{3,i_force(2)} = [num2str(frac(1)*force),'*(y>140*',num2str(s),')'];


% Parse  problem.
fea = parsephys(fea);
fea = parseprob(fea);


% Main loop.
figure
umax = [];
wmax = [];
n    = 10;   % Number of force lines.
p1   = [ 10*s*ones(n,1) linspace(140*s,220*s,n)' zeros(n,1) ];
p2   = [ zeros(n,1) linspace(140*s,220*s,n)' -3*s*ones(n,1) ];
n1   = [ ones(n,1) zeros(n,1) zeros(n,1) ];
for i=1:length(frac)

  % Split forces in x/z-directions according to frac.
  fea.bdr.n{1}{i_force(1)} = ['-',num2str((1-frac(i))*force),'*(y>140*',num2str(s),')'];
  fea.bdr.n{3}{i_force(2)} = [    num2str(   frac(i) *force),'*(y>140*',num2str(s),')'];

  % Solve problem.
  fea.sol.u = solvestat( fea, 'icub', 1+str2num(strrep(sfun,'sflag','')) );

  % Determine maximum x/z-displacements.
  umax = [ umax max(abs(fea.sol.u(unique(fea.eqn.dofm{1}(:)))/s)) ];
  wmax = [ wmax max(abs(fea.sol.u(unique(sum(fea.eqn.ndof(1:2))+fea.eqn.dofm{3}(:)))/s)) ];

  % Create fea struct with displaced grid.
  dp = zeros(size(fea.grid.p));
  for j=1:3
    dp(j,:) = DSCALE*evalexpr( fea.dvar{j}, fea.grid.p, fea );
  end
  fea_disp = fea;
  fea_disp.grid.p = fea_disp.grid.p + dp;

  % Plot deformed object colored by the total displacement.
  clf
  subplot( 1, 5, [1 2 3] )
  postplot( fea_disp, 'surfexpr', ['sqrt(u^2+v^2+w^2)/',num2str(s)], 'colorbar', false, 'linewidth', 0.1 )
  view( [-60 20] )
  axis( [-0.0272 0.0273 -0.0237 0.2577 -0.0050 0.0767] )
  axis off

  % Construct force line visualizations.
  p12 = p1 + DSCALE*[ evalexpr( 'u', p1', fea ) evalexpr( 'v', p1', fea ) evalexpr( 'w', p1', fea ) ];
  hq1 = quiver3( p12(:,1), p12(:,2), p12(:,3), n1(:,1), n1(:,2), n1(:,3), max(eps,1-frac(i)), 'MaxHeadSize', 0, 'LineWidth', lw );
  p22 = p2 + DSCALE*[ evalexpr( 'u', p2', fea ) evalexpr( 'v', p2', fea ) evalexpr( 'w', p2', fea ) ];
  n2  = [ zeros(n,1) 0.4*frac(i)*ones(n,1) -(1-0.4*frac(i))*ones(n,1) ];
  hq2 = quiver3( p22(:,1), p22(:,2), p22(:,3), n2(:,1), n2(:,2), n2(:,3), max(eps,frac(i)), 'MaxHeadSize', 0, 'LineWidth', lw );

  % Plot displacement curves.
  subplot( 1, 5, [4 5] )
  [ax,h1,h2] = plotyy( frac(1:i), umax, frac(1:i), wmax );
  xlabel( 'Lateral loading fraction', 'FontSize', fs )
  ylabel( ax(1), 'Maximum x-displacement [mm]', 'FontSize', fs )
  axis(   ax(1), [0 1 0 3] )
  ylabel( ax(2), 'Maximum z-displacement [mm]', 'FontSize', fs )
  axis(   ax(2), [0 1 0 40] )
  set( ax, 'FontSize', fs )
  set( [h1 h2], 'LineWidth', lw )
  grid on

  set( hq1, 'Color', get( h1, 'Color' ) )
  set( hq2, 'Color', get( h2, 'Color' ) )

  if( ~(exist([pwd,filesep,'tmp'])==7) )
    mkdir( 'tmp' )
  end
  set( gcf, 'Position', [200,200,853*5/3,480], 'PaperUnits', 'centimeters', 'PaperPosition', [0 0 16 9], 'PaperPositionMode', 'auto' );
  print( '-r200', '-dpng', sprintf( 'tmp/img%03i.png', i ) )
end

cd tmp
% system( 'convert -loop 0 img*.png -resize 640x480 -dither none -colors 8 -fuzz "40%" -layers OptimizeFrame output.gif' )
system( 'ffmpeg -y -i img%03d.png -c:v libx264 -vf "fps=25, format=yuv420p, crop=2000:714:163:0" output.mp4' )
cd ..
