% Solution for final exam problem F.3
% crack nucleation model
% requires fact.m

%initialize parameters
mu = 1; nu = 0.3; 
K = mu/4/pi/(1-nu);
A = 0.1; alpha=1; 
B = 0.05; beta = 0.01; C = 0.04; %0.04
d = 1;
N = 500; 
dim = 2*N + 1;
delta = 0.1;

maxiter =  5000;
plotfreq = 200;
clear u 
x=[-N:N]*d;
%applied stress
sig = B * exp(-beta*x.^2) + C;

ind=[0:2*N];
F=fact(ind);
FM=zeros(dim,dim);
for i=1:dim,
    FM(i,:)=F([i:-1:2,1:dim-(i-1)]) + log(d) - 3/2;
end
figure(2); mesh(FM);
disp('FM constructed');


% initialize solution
if(~exist('u'))
  u=zeros(size(x));
end

%relaxation loop
for iter=0:maxiter,
   
rho = ( [u(2:end),u(end)] - u ) / d ;
gamma = A*(1-(1+alpha*u).*exp(-alpha*u))/alpha^2;
dgammadu = A*u.*exp(-alpha*u);
G = rho*FM*rho' * (-K*d^2) + sum(gamma) * d - sum(sig.*u) * d;

dGdrho = 2*rho*FM * (-K*d^2);
dGdu = ([0,dGdrho(1:end-1)]-[dGdrho(1:end-1),0])/d + dgammadu * d - sig * d;

u = u - dGdu * delta;
u(1)=u(2); u(end)=u(end-1);

if mod(iter,plotfreq)==0
   plot(x,u); drawnow
   disp(sprintf('iter=%d C=%f G = %20.12e',iter,C,G));
end

end

