%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to CSE Your Homework Assignment Project 8    %
%  Elastoplastic Torsion: Twist and Stress               %
%  problem4.m Dianne P. O'Leary   04/04                  %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Determine the stress function u(x,y) on the ellipsoidal 
% cross-section of a bar when G is the shear modulus and 
% theta is the angle of twist per unit length.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If the bar is purely elastic,  
% we model the problem using Poisson's equation
%        -u_xx - u_yy =
%        -div(grad(u))= 2 G theta
% on the ellipse
%    (x/myalpha)^2 + (y/mybeta)^2 <= 1
% When discretized, this problem is equivalent to 
% minimizing a potential function
%     min 1/2 u^T K u - u^T rhs
% where K is a discrete approximation to the Laplace 
% operator and rhs is a vector determined by the 
% right-hand side of the differential equation.  
% Setting the gradient to zero, we find that the
% solution is determined by K u = rhs.
% (See Problem 1.)
%
% Problem 2:
% Now suppose that the bar is elasto-plastic, and the norm of 
% the gradient of the stress function (called the shear 
% stess) must be kept less than the yield stress.  
% We solve this by minimizing the potential function
% (actually, a discretization of it) subject to keeping the
% stress function's value bounded by its distance to the 
% boundary of the ellipse:
%     min 1/2 u^T K u - u^T rhs
%        -d <= u <= d
% where d is the distance between the mesh point and the
% boundary of the ellipse.
% For physical reasons, the lower bounds of -d can be
% replaced by 0 without changing the solution, and these
% lower bounds will not be active except at the boundary.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Global variables:
% myalpha, mybeta = parameters defining the ellipse
%    (x/myalpha)^2 + (y/mybeta)^2 <= 1
% myz = parameter used in find_dist
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global myalpha mybeta myz

clc

ratio = [1 .8 .65 .5 .2];

a = 1 ./ ratio;

figure(1)
jplot = 0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Loop over various aspect ratios for the ellipse
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for ii=1:length(ratio),
   disp(sprintf('Aspect ratio for ellipse = %f',ratio(ii)))
   mybeta = 1;
   myalpha = a(ii);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem definition: g gives the boundary definition
% and bndrycond gives the Dirichlet boundary condition.
% We use Matlab's PDE Toolbox to form the matrix K.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   g='ellipse'; 
   bndrycond='circleb1'; % 0 on the boundary
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The equation specified in assemblepde is 
%              -div(c*grad(u))+a*u=f .
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   pdec=1;
   pdea=0;

   G = 1;
   theta = 1;
   f=2*G*theta;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate a coarse initial mesh.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   [pmesh,emesh,tmesh]=initmesh(g,'hmax',1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Refine the mesh uniformly several times.
% The function assempde gives the mass matrix K and 
% the right-hand side rhs to specify the solution at 
% the internal nodes, and then the full solution is 
% obtained using the node information in the
% matrix B and the boundary information in ud.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   for i=1:3,
     [pmesh,emesh,tmesh]=refinemesh(g,pmesh,emesh,tmesh);
     [K,rhs,B,ud]=assempde(bndrycond,pmesh,emesh,tmesh, ...
                           pdec,pdea,f);
   end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve for the distances d between each interior mesh 
% point and the boundary of the ellipse.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   is = 0;
   [m,n]=size(B);
   for i=1:m
     if (norm(B(i,:))>0) % this is an interior point
       is = is + 1;
       d(is,1) = dist_to_ellipse(pmesh(:,i),10*eps);
     end
   end
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the quadratic programming problem, assemble the 
% solution, and plot.
% Note that rhs'*unc is an estimate of the integral of u
% over the domain.
% The variable atheta is alpha times theta and G=1.
% athetelas records the highest atheta value that leads
%   to an elastic solution (no variables at bounds).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%  sigma_0 = 1
   atheta = 0:.25:5;
   for jj=1:length(atheta)
      unc = quadprog(K,-atheta(jj)/a(ii)*rhs, ...
                     [],[],[],[],zeros(is,1),d);
      if ((jj == 2)|(jj == 4)) &(ii>1)
         jplot = jplot + 1;
         subplot(4,2,jplot)
         uc = B*unc + ud; 
         pdeplot(pmesh,emesh,tmesh,'xydata',uc),axis equal;
         caxis([0,.35])
      end
      torque(ii,jj) = rhs'*unc/(myalpha^3);
      count = sum(unc<d-1.e-8); % = number of variables not at bound
      if (count == n)
         athetelas(ii) = atheta(jj);
         torqelas(ii) = torque(ii,jj);
      end
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot Torque/(sigma_0 alpha^3) vs G alpha theta / sigma_0.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2)

gg(1).color = 'b';
gg(2).color = 'g-.';
gg(3).color = 'm:';
gg(4).color = 'k--';
gg(5).color = 'c';

for ii=1:length(ratio)
   lab(ii,:) = sprintf('b/a = %5.2f',ratio(ii));
   plot(atheta,torque(ii,:),gg(ii).color)
   hold on
end
plot(athetelas,torqelas,'r*')
xlabel('G \alpha \theta / \sigma_0')
ylabel('Torque/ ( \sigma_0 \alpha^3)')
legend(lab) 

figure(1)
print -depsc p4soln.eps
figure(2)
print -depsc p4torq.eps