%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Solution to CSE Your Homework Assignment Project 8    %
%  Elastoplastic Torsion: Twist and Stress               %
%  problem1.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.
%
% In this first problem, the bar is purely elastic.  
% We model this case using Poisson's equation
%        -u_xx - u_yy =
%        -div(grad(u))= 2 G theta
% on the ellipse
%    (x/myalpha)^2 + (y/mybeta)^2 <= 1 
% with zero boundary conditions.
% 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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

myalpha = 1;
mybeta = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Problem definition: g gives the boundary definition
% and bndrycond gives the Dirichlet boundary condition.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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 = 5;
theta = 1;
f=2*G*theta;

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Refine the mesh uniformly several times,
% computing the approximate solution each time.
% 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, both set by assempde.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:4,
  [pmesh,emesh,tmesh]=refinemesh(g,pmesh,emesh,tmesh);
  [K,rhs,B,ud]=assempde(bndrycond,pmesh,emesh,tmesh,pdec,pdea,f);
  un = K \ rhs;
  u = B*un + ud;  % ud=0 for this problem
  energy(i) = .5*un'*K*un - un'*rhs;
  if (i==3)
     figure(1)
     pdeplot(pmesh,emesh,tmesh),axis equal
     xlabel('x')
     ylabel('y')
     figure(2)
     pdeplot(pmesh,emesh,tmesh,'xydata',u),axis equal;
     xlabel('x')
     ylabel('y')
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% See how energy estimates change with the grid,
% approximating the solution by the solution for i=4.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:3,
  disp(sprintf('Grid %d, Est. Energy = %f, Est Error = %e',...
        i, energy(i), energy(i)-energy(4)))
end

figure(1)
print -depsc p1mesh.eps
figure(2)
print -depsc p1soln.eps