clc; clear
addpath('..\..\..\SharedFiles');
sideLengthX = 0.10;
sideLengthY = 0.04;
ndivx = 401;
ndivy = 161;
mvalue = 4;
dx = sideLengthX/ndivx;
dy = sideLengthY/ndivy;
maxdx = max(dx, dy);
xl = linspace(0, sideLengthX, ndivx+1);
yl = linspace(-1/2*sideLengthY, 1/2*sideLengthY, ndivy+1);
[x, y] = meshgrid(xl, yl);
[x1, y1] = meshgrid(xl, -1/2*sideLengthY - (1:mvalue)*dy);
[x2, y2] = meshgrid(xl, +1/2*sideLengthY + (1:mvalue)*dy);
coord  = [x(:), y(:); x1(:), y1(:); x2(:), y2(:)];
[np, nd] = size(coord);
thickness = 1;
particleVolume = (dx*dy*thickness) * ones(np,1);
delta = maxdx * (mvalue+0.015);
[ih, jh] = gridBasedNonUniformHorizonSearch(coord, delta, 1.5, 'no');
xi  = coord(jh,:)-coord(ih,:);
bl0 = sqrt(sum(xi.^2,2));
vc = (bl0<delta-maxdx/2)*1 + ...
    (bl0>delta-maxdx/2).*(delta-bl0+maxdx/2)/dx;
vj = particleVolume(jh).*vc;
wf = exp(-bl0.^2/delta^2);
wfvj = wf.*vj;

ap = [0;find(diff(ih)==1); numel(ih)];
g1 = firstOrderGfun2D(xi, wfvj, ap);
g2 = secondOrderGhatfun2D(xi, wfvj, ap);

% field quantities
dis = zeros(np, nd);
vel = zeros(np, nd);
acc = zeros(np, nd);
rot = zeros(np, 1);
rotvel = zeros(np, 1);
rotacc = zeros(np, 1);
bodyForce = zeros(np, nd);
bodyCouple = zeros(np, 1);
rot_old = rot;
dis_old = dis;

% material properties
emod   = 72e9;
nu     = 1/3;
dens   = 2440;
lambda = nu * emod / ((1+nu)*(1-2*nu));
mu     = emod / (2*(1+nu));
lambda = lambda * (1-2*nu) / (1-nu);

couplingNumber = 0.00;
mu_r   = couplingNumber^2/(1-couplingNumber^2) * mu;
characteristicLength = 0.1*sideLengthY;
gamma  = 4 * mu * characteristicLength^2;
microInertia = 1e-2;
pres = 2e6;

nt=1000000;
fail = ~preCrack2D(coord, ih, jh, [-sideLengthX/2,0], [sideLengthX/2-dx/100,0]);
dmg0 = 1 - accumarray(ih,fail) ./ accumarray(ih,true(size(fail)));

upNode    = coord(:,2) > +sideLengthY/2+dx/100;
downNode  = coord(:,2) < -sideLengthY/2-dx/100;
rightNode = coord(:,1) > sideLengthX - dx/100;

wfvj0 = wfvj;
nonfaillayers = 2;
nonfailzone = ...
    coord(ih,2) > +sideLengthY/2-nonfaillayers*mvalue*dx + dx/100 | ...
    coord(ih,2) < -sideLengthY/2+nonfaillayers*mvalue*dx - dx/100 | ...
    coord(jh,2) > +sideLengthY/2-nonfaillayers*mvalue*dx + dx/100 | ...
    coord(jh,2) < -sideLengthY/2+nonfaillayers*mvalue*dx - dx/100;


t0   = 0;
recordTime = 0.05e-6;
nrecord = 1;
resFolder = sprintf('Result_pr_1_3\\CN_%g_LB_%g_J0_%g_S0_%g\\', couplingNumber, ...
    characteristicLength*1e3, microInertia, pres/1e6);
if(~exist(resFolder, 'dir'))
    mkdir(resFolder);
elseif(exist([resFolder,'step1'], 'file'))
    fprintf('Already Setup Folder and numerical implementation\n');
    return;
end
init.x    = x;
init.y    = y;
init.wfvj = wfvj;
init.g    = g1;
init.ghat = g2;
init.dmg0 = dmg0;
init.emod = emod;
init.mu   = mu;
init.mu_r = mu_r;
init.lb   = characteristicLength;
init.gamma = gamma;
init.nu = nu;
init.microInertia = microInertia;
init.pres = pres;
save([resFolder,'init'], 'init');
recordFlag = false;
currentTime = 0;
for tt = 1:1:nt
    % 能用sqrt(sum())就不用vecnorm
    % 能用accumarray就不用full(sparse())
    % 能用jh的地方不另外搞一个j的东西
    dt = 0.25*dx*min(1/sqrt(max(lambda+2*mu, mu+mu_r)/dens), ...
        1/sqrt(gamma/dens/microInertia));
    currentTime = currentTime + dt;
    if(nrecord*recordTime<currentTime)
        dt = dt - (currentTime - nrecord*recordTime);
        currentTime = nrecord*recordTime;
        recordFlag = true;
        nrecord = nrecord + 1;
    end
    
    if(currentTime<t0)
        bodyForce(upNode,2) = +pres*currentTime/t0*sideLengthX*thickness/sum(particleVolume(upNode));
        bodyForce(downNode,2) = -pres*currentTime/t0*sideLengthX*thickness/sum(particleVolume(downNode));
    else
        bodyForce(upNode,2) = +pres*1*sideLengthX*thickness/sum(particleVolume(upNode));
        bodyForce(downNode,2) = -pres*1*sideLengthX*thickness/sum(particleVolume(downNode));
    end
    
    r   = coord + dis;
    rotavg = 0.5 * (rot(jh,1)+rot(ih,1));
    eta = dis(jh,:)-dis(ih,:) + [xi(:,2).*rotavg, -xi(:,1).*rotavg];
    bl  = sqrt(sum((xi+eta).^2,2));
    
    g1_ = repmat(wfvj0, 1, size(g1,2)) .* g1;
    g2_ = repmat(wfvj0, 1, size(g2,2)) .* g2;
    
    theta = accumarray(ih,g1_(:,1).*eta(:,1) + g1_(:,2).*eta(:,2));
    s     = (bl-bl0)./bl0;
    scr = sqrt(4*pi*3.8/(9*emod*delta));
    fail  = (fail & s<scr) | nonfailzone;
    wfvj0 = wfvj0.*fail;
    dmg = 1 - accumarray(ih,wfvj0) ./ accumarray(ih,wfvj);
    
    Tij   = repmat((lambda-mu+mu_r)*theta(ih),1,nd).*g1_ + ...
        (mu+mu_r)*[(g2_(:,1)+g2_(:,2)).*eta(:,1), (g2_(:,1)+g2_(:,2)).*eta(:,2)] + ...
        (mu-mu_r)*[2*g2_(:,1).*eta(:,1)+g2_(:,3).*eta(:,2), 2*g2_(:,2).*eta(:,2)+g2_(:,3).*eta(:,1)];
    Mij   = gamma * (g2_(:,1)+g2_(:,2)) .* (rot(jh,1)-rot(ih,1));
    Mt  = (xi(:,1).*Tij(:,2) - xi(:,2).*Tij(:,1)) * 0.5;
    np = size(coord,1);
    Force = zeros(np,nd);
    for i = 1:nd
        Force(:,i) = accumarray(ih,Tij(:,i))-accumarray(jh,Tij(:,i));
    end
    Couple  = accumarray(ih,Mij(:,1))-accumarray(jh,Mij(:,1)) + accumarray(ih,Mt(:,1))+accumarray(jh,Mt(:,1)) ;
    
    dis_new = dt^2/dens * (Force+bodyForce) + 2*dis - dis_old;
    dis_old = dis;
    dis     = dis_new;
    vel     = (dis_new - dis_old)/(2*dt);
    
    rot_new = dt^2/(dens * microInertia)* (Couple+bodyCouple) + 2*rot - rot_old;
    rot_old = rot;
    rot     = rot_new;
    rotvel  = (rot_new - rot_old)/(2*dt);
    
    fprintf('tt=%d\n',tt);
    if(recordFlag)
        step.currentTime = currentTime;
        np = numel(x);
        step.U1 = reshape(dis(1:np,1), size(x));
        step.U2 = reshape(dis(1:np,2), size(x));
        step.W3 = reshape(rot(1:np,1), size(x));
        step.dmg = reshape(dmg(1:np,1), size(x));
        recordFlag = false;
        save([resFolder,sprintf('step%d', nrecord-1)], 'step');
    end
    
    if(max(dmg(rightNode,1))>0.3)
        fprintf('Crack propagation finished!\n')
        init.nrecord = nrecord-1;
        save([resFolder,'init'], 'init');
        break;
    end
end