%-----------------------------------------------------------------------
% Test case and validation of the non-singular displacement gradient
% for DDD-XRD calculations. See Bertin and Cai, CMS, 2018.
%
% This code does the following:
%  - generates a random triangular dislocation loop
%  - computes the displacement gradient field along a line using the
%    non-singular expression provided in Bertin and Cai, CMS, 2018.
%  - compares it with the displacement gradient field obtained by
%    numerically differentiating the displacement field
%  - compares the stress obtained from the displacement gradient field
%    with the non-singular stress expression provided in Cai et al., 2006
%
% This code uses the DDLab data structure to represent the dislocations.
% 
% DDD-XRD
% nbertin@stanford.edu
%-----------------------------------------------------------------------

%---------------------------------------------------------
% INITIALIZATION
%---------------------------------------------------------

% Elasticity parameters
MU = 54.6e9;
NU = 0.324;
LAMBDA = 2*MU*NU/(1-2*NU);
a = 1.0; % Non-singular radius

% Initialize triangular loop
L = 1000;
rn = L*rand(3);
b = [1,1,0];
b = b/norm(b);
n = [1,1,1];
n = n/norm(n);
links = [[1:3]',[[2:3]';1]];
links = [links,repmat(b,3,1),repmat(n,3,1)];

% Plot dislocation loop
clf(figure(1));
hold on;
for i=1:size(links,1)
    n1 = links(i,1);
    n2 = links(i,2);
    plot3(rn([n1,n2],1),rn([n1,n2],2),rn([n1,n2],3),'r.-');
end
hold off;
view(3);


%---------------------------------------------------------
% EVALUATE DISPLACEMENT GRADIENT
%---------------------------------------------------------

% Evaluate strain gradient from non-singular expression 
% along a x-line passing through the center of the loop
N = 50; % number of field points
r = zeros(N,3);
r(:,1) = linspace(0,L,N);
r(:,2) = mean(rn(:,2));
r(:,3) = mean(rn(:,3));
figure(1);
hold on;
plot3(r(:,1),r(:,2),r(:,3),'bo-');
hold off;
dudx = displacement_gradient_structure(rn,links,NU,a,r);

% Evaluate displacement gradient by numerically differentiating
% the displacement field
delta = 2;
rd = r; rd(:,1) = rd(:,1)-delta;
u10 = displacement_structure(rn,links,NU,a,rd);
rd = r; rd(:,1) = rd(:,1)+delta;
u11 = displacement_structure(rn,links,NU,a,rd);
dudx1 = (u11-u10)/2/delta;
rd = r; rd(:,2) = rd(:,2)-delta;
u20 = displacement_structure(rn,links,NU,a,rd);
rd = r; rd(:,2) = rd(:,2)+delta;
u21 = displacement_structure(rn,links,NU,a,rd);
dudx2 = (u21-u20)/2/delta;
rd = r; rd(:,3) = rd(:,3)-delta;
u30 = displacement_structure(rn,links,NU,a,rd);
rd = r; rd(:,3) = rd(:,3)+delta;
u31 = displacement_structure(rn,links,NU,a,rd);
dudx3 = (u31-u30)/2/delta;

% Plot and compare displacement gradients
figure(2);
subplot(3,3,1);
plot(r(:,1),squeeze(dudx(:,1,1)),'-',r(:,1),squeeze(dudx1(:,1)),'+'); title('du_1/dx_1');
subplot(3,3,2);
plot(r(:,1),squeeze(dudx(:,1,2)),'-',r(:,1),squeeze(dudx2(:,1)),'+'); title('du_1/dx_2');
subplot(3,3,3);
plot(r(:,1),squeeze(dudx(:,1,3)),'-',r(:,1),squeeze(dudx3(:,1)),'+'); title('du_1/dx_3');

subplot(3,3,4);
plot(r(:,1),squeeze(dudx(:,2,1)),'-',r(:,1),squeeze(dudx1(:,2)),'+'); title('du_2/dx_1');
subplot(3,3,5);
plot(r(:,1),squeeze(dudx(:,2,2)),'-',r(:,1),squeeze(dudx2(:,2)),'+'); title('du_2/dx_2');
subplot(3,3,6);
plot(r(:,1),squeeze(dudx(:,2,3)),'-',r(:,1),squeeze(dudx3(:,2)),'+'); title('du_2/dx_3');

subplot(3,3,7);
plot(r(:,1),squeeze(dudx(:,3,1)),'-',r(:,1),squeeze(dudx1(:,3)),'+'); title('du_3/dx_1');
subplot(3,3,8);
plot(r(:,1),squeeze(dudx(:,3,2)),'-',r(:,1),squeeze(dudx2(:,3)),'+'); title('du_3/dx_2');
subplot(3,3,9);
plot(r(:,1),squeeze(dudx(:,3,3)),'-',r(:,1),squeeze(dudx3(:,3)),'+'); title('du_3/dx_3');


%---------------------------------------------------------
% EVALUATE STRESS FIELD
%---------------------------------------------------------

% Evaluate stress from displacement gradient
Sdudx = zeros(size(dudx));
for i=1:size(dudx,1)
    Gkk = dudx(i,1,1)+dudx(i,2,2)+dudx(i,3,3);
    Sdudx(i,1,1) = 2*MU*dudx(i,1,1)+LAMBDA*Gkk;
    Sdudx(i,2,2) = 2*MU*dudx(i,2,2)+LAMBDA*Gkk;
    Sdudx(i,3,3) = 2*MU*dudx(i,3,3)+LAMBDA*Gkk;
    Sdudx(i,1,2) = MU*(dudx(i,1,2)+dudx(i,2,1));
    Sdudx(i,1,3) = MU*(dudx(i,1,3)+dudx(i,3,1));
    Sdudx(i,2,3) = MU*(dudx(i,2,3)+dudx(i,3,2));
end

% Evaluate stress field from analytical expression
Sns = stress_structure(rn,links,MU,NU,a,r);

% Plot and compare stress fields
figure(3);
subplot(2,3,1);
plot(r(:,1),squeeze(Sns(:,1,1)),'-',r(:,1),squeeze(Sdudx(:,1,1)),'+'); title('S_{11}');
subplot(2,3,2);
plot(r(:,1),squeeze(Sns(:,2,2)),'-',r(:,1),squeeze(Sdudx(:,2,2)),'+'); title('S_{22}');
subplot(2,3,3);
plot(r(:,1),squeeze(Sns(:,3,3)),'-',r(:,1),squeeze(Sdudx(:,3,3)),'+'); title('S_{33}');

subplot(2,3,4);
plot(r(:,1),squeeze(Sns(:,1,2)),'-',r(:,1),squeeze(Sdudx(:,1,2)),'+'); title('S_{12}');
subplot(2,3,5);
plot(r(:,1),squeeze(Sns(:,1,3)),'-',r(:,1),squeeze(Sdudx(:,1,3)),'+'); title('S_{13}');
subplot(2,3,6);
plot(r(:,1),squeeze(Sns(:,2,3)),'-',r(:,1),squeeze(Sdudx(:,2,3)),'+'); title('S_{23}');
