% FEATool Multiphysics axisymmetric stress-strain with custom equation parsing.

% Copyright 2013-2019 Precise Simulation, Ltd.
close all, clear all, clc


test_case = 2;   % 0 = Single triangle test.
                 % 1 = Pressurized hollow cylinder.
                 % 2 = Pressurized hollow sphere.
                 % 3 = Disc with fixed edge and central point load.


sfun = 'sflag1';
p    = 20e4;    % Load force.
nu   = 0.3;
E    = 200e9;


% Geometry and grid.
switch( test_case )

  case 0   % Single triangle test.
    fea.grid.p = [ 3 4 3 ;
                   0 0 1 ];
    fea.grid.c = [1:3]';
    fea.grid.a = zeros(3,1);
    fea.grid.s = 1;
    fea.grid.b = gridbdr( fea.grid.p, fea.grid.c, fea.grid.a );

  case 1   % Pressurized hollow cylinder.
    a = 1;
    b = 2;
    fea.grid = rectgrid( 10, 10, [a b;0 1] );

  case 2   % Pressurized hollow sphere.

    a = 1;
    b = 2;
    fea.grid = ringgrid(12,72,a,b);
    fea.grid = delcells( fea.grid, selcells( fea.grid, '(x<=eps) + (y<=eps)') );

  case 3   % Disc with fixed edge and central point load.

    p  = 10;
    nu = 1/3;
    E  = 1e3;

    a = 0;
    b = 10;
    fea.grid = rectgrid( 20, 20, [a b;-0.5 0.5] );

end


% Equation definitions.
eqn_r = '- c1*( (1-nu)*ur_r + nu*wz_r + nu/r*( u_r - ur_t - wz_t ) - (1-nu)/r^2*u_t ) - c1*c2*( uz_z + wr_z ) = 2*pi*r*Fr';
eqn_z = '- c1*( nu*ur_z + nu/r*u_z + (1-nu)*wz_z ) - c1*c2*( uz_r + wr_r ) = 2*pi*r*Fz';


% fea struct.
fea.sdim = { 'r' 'z' };
fea.dvar = { 'u' 'w' };
fea.sfun = { sfun sfun };
fea.eqn  = parseeqn( { eqn_r eqn_z }, fea.dvar, fea.sdim );


% Constants and expressions.
c1 = E / ((1+nu)*(1-2*nu) );
c2 = (1-2*nu) / 2;
fea.expr = { 'c1' [] [] {[num2str(c1*2*pi),'*r']} ;
             'c2' [] [] {c2} ;
             'nu' [] [] {nu} ;
             'Fr' [] [] {0} ;
             'Fz' [] [] {0} };


% Boundary conditions.
switch( test_case )

  case 0
    fea.bdr.d = { [] [] [] ;
                  0  [] [] };
    fea.bdr.n = { [] [num2str(p*2*pi),'*r'] [] ;
                  [] [] [] };

  case 1
    fea.bdr.d = { [] [] [] [] ;
                  0  0  0  0  };
    fea.bdr.n = { [] [] [] [num2str(p*2*pi),'*r'] ;
                  [] [] [] [] };

  case 2

    fea.bdr.d = { [] [] []  [0] ;
                  [] [] [0] []  };
    fea.bdr.n = { [num2str(p*2*pi),'*r*-nr'] [] [] [] ;
                  [num2str(p*2*pi),'*r*-nz'] [] [] [] };

  case 3
    fea.bdr.d = { [] [] [] []  ;
                  [] [] [] [] };
    fea.bdr.n = { [] [] [] []    ;
                  [] [] [] [-p] };

    % Set point constraints.
    [~,ix] = min( (fea.grid.p(1,:)-b).^2 + (fea.grid.p(2,:)-0).^2 );
    fea.pnt(1).index = ix;
    fea.pnt(1).type  = 'constr';
    fea.pnt(1).dvar  = 2;
    fea.pnt(1).expr  = 0;

end


% Solve.
fea       = parseprob( fea );
fea.sol.u = solvestat( fea, 'icub', 5 );


% Postprocessing.
switch( test_case )

  case 0

    syms r z real positive

    B = [ -1  0 1 0 0 0 ;
          (4-r-z)/r 0 (-3+r)/r 0 z/r 0 ;
           0 -1 0 0 0 1 ;
           -1 -1 0 1 1 0 ];

    C = E/((1+nu)*(1-2*nu))*[ 1-nu   nu   nu          0 ;
                                nu 1-nu   nu          0 ;
                                nu   nu 1-nu          0 ;
                                0    0    0  (1-2*nu)/2 ];
    K = simplify( B'*C'*B*r );
    A = eval( 2*pi*int( int( K, z, 0, 4-r ), r, 3, 4 ) );
    A([2 4],:) = 0; A(:,[2 4]) = 0; A(2,2) = 1; A(4,4) = 1;

    f = zeros(6,1);
    f([1 5]) = 3*pi*p;

    u_ref = A\f;

    err_norm = norm( fea.sol.u - u_ref([1 3 5 2 4 6]) )

  case 1

    n = 20;
    r = linspace(a,b,n);
    z = 0.5*ones(1,n);
    u_ref = p * a^2*(1+nu)*(b^2+r.^2*(1-2*nu)) ./ (E*(b^2-a^2)*r);
    u = evalexpr( 'u', [r;z], fea );

    subplot(1,2,1)
    postplot( fea, 'surfexpr', 'u' )
    title('r-displacement')
    subplot(1,2,2), hold on
    plot(u_ref,r,'r-')
    plot(u,r,'b.')
    legend('exact solution','computed solution')
    xlabel('r')
    grid on

  case 2

    n = 20;
    r = linspace(a,b,n);
    z = zeros(1,n);
    u_ref = 1./(2*E*(b^3-a^3)*r.^2) .* (2*(p*a^3)*(1-2*nu)*r.^3+p*(1+nu)*b^3*a^3);
    u = evalexpr( 'u', [r;z], fea );

    subplot(1,2,1)
    postplot( fea, 'surfexpr', 'sqrt(u^2+w^2)', 'arrowexpr', {'u' 'w'} )
    title('computed displacement')
    subplot(1,2,2), hold on
    plot(u_ref,r,'r-')
    plot(u,r,'b.')
    title( 'radial displacement')
    legend('exact solution','computed solution')
    xlabel('r')
    grid on

  case 3

    subplot(2,2,1)
    postplot( fea, 'surfexpr', 'u' )
    title('computed r-displacement')

    subplot(2,2,2)
    D = E*1^3/(12*(1-nu^2));
    u_ref = [num2str(-p),'/(8*pi*',num2str(D),')*(',num2str((3+nu)/(1+nu)-1),'-2*log(r/',num2str(b),'))*r*z'];

    postplot( fea, 'surfexpr', u_ref )
    title('exact r-displacement')

    subplot(2,2,3)
    postplot( fea, 'surfexpr', 'w' )
    title('computed z-displacement')

    subplot(2,2,4)
    D = E*1^3/(12*(1-nu^2));
    u_ref = [num2str(-p),'/(16*pi*',num2str(D),')*(',num2str((3+nu)/(1+nu)),'*(',num2str(b^2),'-r^2)+2*r^2*log(r/',num2str(b),'))'];
    postplot( fea, 'surfexpr', u_ref )
    title('exact z-displacement')

end
