function [X, Y, VOL, U1, U2, W3, U_L2, W_L2, F_L2, M_L2] = ...
    MP_quasistatic_Model_I_True(mvalue, ndiv, ajustRatio, gfunOrder, ntMax, Ftol, Utol)
addpath('../../../SharedFiles')
sideLengthX = 1;
sideLengthY = 1;
ndivx = ndiv;
ndivy = ndiv;
x = linspace(0, sideLengthX, ndivx+1);
y = linspace(0, sideLengthY, ndivy+1);
dx = sideLengthX/ndivx;
dy = sideLengthY/ndivy;
maxdx = max(dx, dy);
[x, y] = meshgrid(x, y);
coord  = [x(:), y(:)];
[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)];
if(gfunOrder==1)
    g1 = firstOrderGfun2D(xi, wfvj, ap);
elseif(gfunOrder==2)
    G = secondOrderPDDOfun(xi, ap, wfvj, 2);
    g1 = G(:,1:2);
end
% g1_ = repmat(wfvj, 1, size(g1,2)) .* g1;
% g2 = secondOrderGhatfun2D(xi, wfvj, ap);
% g2_ = repmat(wfvj, 1, size(g2,2)) .* g2;

% 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   = 1;
nu     = 0.2;
% dens   = 1;
lambda = nu * emod / ((1+nu)*(1-2*nu));
mu     = emod / (2*(1+nu));
couplingNumber = 0.5;
mu_r   = couplingNumber^2/(1-couplingNumber^2) * mu;
characteristicLength = sideLengthY*1;
gamma  = 2 * mu * characteristicLength^2;
% microInertia = 1;

nt = ntMax;
crackLength = 0.0;
fail = abs(coord(ih,1) - 0.5*sideLengthX)<0.5*crackLength & ...
    abs(coord(jh,1) - 0.5*sideLengthX)<0.5*crackLength & ...
    (coord(ih,2) - 0.5*sideLengthY) .*(coord(jh,2) - 0.5*sideLengthY) <= 0;
dmg = accumarray(ih,fail) ./ accumarray(ih,true(size(fail)));
% scatter(coord(:,1), coord(:,2), 20, dmg(:,1), 'filled')
% colormap(jet)
% colorbar
% axis equal
% pause(0.001)
loadThickRatio = 1/3.015;
upNode    = coord(:,2) > sideLengthY-loadThickRatio*delta;
leftNode  = coord(:,1) < loadThickRatio*delta;
rightNode = coord(:,1) > sideLengthX-loadThickRatio*delta;
downNode  = coord(:,2) < loadThickRatio*delta;
innerNode = ~(upNode|downNode);
vel_y = 0.01;
nonfailzone = coord(ih,2)>sideLengthY+5*delta | ...
    coord(ih,2) < -5*delta;
wfvj0 = wfvj;

%
pforce    = zeros(np,2);
pforceold = zeros(np,2);
pcouple    = zeros(np,1);
pcoupleold = zeros(np,1);
velhalfold = zeros(np,2);
rotvelhalfold = zeros(np,1);
rotvelhalf = zeros(np,1);
bforce = zeros(np,2);
bcouple = zeros(np,1);

% massvec1 = 8*mu / delta^2 * dt^2 * (0.25) * 500;
% massvec2 = pi*delta^2 * massvec1 + 4*gamma/delta^2 * dt^2 * 0.25 * 500;

% dt = 0.5*min([dx*sqrt(massvec1/emod), dx*sqrt(massvec2/gamma)]);

dt = 1;
E_max_u = max(lambda+2*mu, mu+mu_r);
E_max_theta = max(gamma, gamma);
c_CFL = 0.5;
massvec1 = dt^2 / (c_CFL * maxdx)^2 * E_max_u;
massvec2 = dt^2 / (c_CFL * maxdx)^2 * E_max_theta;
U_L2 = nan(nt,1);
F_L2 = nan(nt,1);
W_L2 = nan(nt,1);
M_L2 = nan(nt,1);
for tt = 1:1:nt
    % 能用sqrt(sum())就不用vecnorm
    % 能用accumarray就不用full(sparse())
    % 能用jh的地方不另外搞一个j的东西
    currentTime = tt * dt;
    if(vel_y * currentTime<0.01*sideLengthY)
        dis(upNode,2)   = 0 * currentTime;
        dis(upNode,1)   = vel_y * currentTime;
        dis(downNode,2) = 0;
        dis(downNode,1) = 0;
        rot(upNode,1)   = 0;
        rot(downNode,1) = 0;
    else
        dis(upNode,2)   = 0;
        dis(upNode,1)   = 0.01*sideLengthY;
        dis(downNode,2) = 0;
        dis(downNode,1) = 0;
        rot(upNode,1)   = 0;
        rot(downNode,1) = 0;
    end
    %     dis = 0.1 * coord + [0.2*coord(:,2), 0.3*coord(:,1)];
    rot(~innerNode) = 0;
    %r   = coord + dis;
    rotavg = 0.5 * (rot(jh,1)+rot(ih,1));
    drot  = rot(jh,1)-rot(ih,1);
    eta = dis(jh,:)-dis(ih,:) + [xi(:,2).*rotavg, -xi(:,1).*rotavg];
    bl  = sqrt(sum((xi+eta).^2,2));
    %wfvj0(fail) = 0;
    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;
    fail  = fail | s>0.1 & (~nonfailzone);
    
    exx = accumarray(ih, g1_(:,1).*eta(:,1));
    exy = accumarray(ih, g1_(:,1).*eta(:,2));
    eyx = accumarray(ih, g1_(:,2).*eta(:,1));
    eyy = accumarray(ih, g1_(:,2).*eta(:,2));
    kxz = accumarray(ih, g1_(:,1).*drot);
    kyz = accumarray(ih, g1_(:,2).*drot);
    sxx = (lambda+2*mu) * exx + lambda*eyy;
    syy = (lambda+2*mu) * eyy + lambda*exx;
    sxy = (mu+mu_r) * exy + (mu-mu_r) * eyx;
    syx = (mu+mu_r) * eyx + (mu-mu_r) * exy;
    mxz = gamma * kxz;
    myz = gamma * kyz;
    
    zu = eta - [exx(ih).*xi(:,1) + eyx(ih).*xi(:,2), exy(ih).*xi(:,1) + eyy(ih).*xi(:,2)];
    zt = drot - (kxz(ih).*xi(:,1) + kyz(ih).*xi(:,2));
    
    bl = sqrt(sum(xi.^2,2));
    m   = accumarray(ih,wfvj.*bl.^2);
    
    cxx = ajustRatio * 8*(mu-mu_r)*wfvj.*xi(:,1).^2./(bl.^2.*m(ih));
    cxy = ajustRatio * 8*(mu-mu_r)*wfvj.*xi(:,1).*xi(:,2)./(bl.^2.*m(ih));
    cyx = ajustRatio * 8*(mu-mu_r)*wfvj.*xi(:,2).*xi(:,1)./(bl.^2.*m(ih));
    cyy = ajustRatio * 8*(mu-mu_r)*wfvj.*xi(:,2).*xi(:,2)./(bl.^2.*m(ih));
    Tij_s = [cxx.*zu(:,1) + cxy.*zu(:,2),cyx.*zu(:,1) + cyy.*zu(:,2)] + ajustRatio * 4 *mu_r * [wfvj.*zu(:,1)./m(ih), wfvj.*zu(:,2)./m(ih)];
    Mij_s = ajustRatio * 2*gamma*wfvj.*zt(:,1)./(m(ih));
    sxxp = accumarray(ih,Tij_s(:,1).*xi(:,1));
    syxp = accumarray(ih,Tij_s(:,2).*xi(:,1));
    sxyp = accumarray(ih,Tij_s(:,1).*xi(:,2));
    syyp = accumarray(ih,Tij_s(:,2).*xi(:,2));
    mxzp = accumarray(ih,Mij_s(:,1).*xi(:,1));
    myzp = accumarray(ih,Mij_s(:,1).*xi(:,2));
    sxx = sxx - sxxp;
    syy = syy - syyp;
    sxy = sxy - sxyp;
    syx = syx - syxp;
    mxz = mxz - mxzp;
    myz = myz - myzp;
    Tij = [sxx(ih).*g1_(:,1) + syx(ih).*g1_(:,2), sxy(ih).*g1_(:,1) + syy(ih).*g1_(:,2)] + Tij_s; % + 10*mu/(1.19*delta^3)*zu.*[wfvj, wfvj];
    Mij = mxz(ih).*g1_(:,1) + myz(ih).*g1_(:,2) + Mij_s; % + 10*gamma/(1.19*delta^3)*zt.*wfvj;
    Mt  = (xi(:,1).*Tij(:,2) - xi(:,2).*Tij(:,1)) * 0.5;
    % Mij = Mij + Mt;
    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)) ;
    
    % 动态松弛法，参照Madenci的程序
    pforceold = pforce;
    pforce = Force;
    pcoupleold = pcouple;
    pcouple = Couple;
    cn_u = zeros(np, 1);
    for i = 1:np
        %kii = zeros(1,nd);
        cn1_u = 0;
        cn2_u = 0;
        for j = 1:nd
            if(abs(velhalfold(i,j))>1e-15)
                cn1_u = cn1_u - dis(i,j)*dis(i,j)*(pforce(i,j)/massvec1-pforceold(i,j)/massvec1) / (dt * velhalfold(i,j));
            end
            cn2_u = cn2_u + dis(i,j)*dis(i,j);
        end
        if(abs(cn2_u)>1e-12)
            if ((cn1_u / cn2_u) > 0.0)
                cn_u(i) = 2.0 * sqrt(cn1_u / cn2_u);
            else
                cn_u(i) = 0.0;
            end
        else
            cn_u(i) = 0.0;
        end
        if (cn_u(i) > 2.0d0)
            cn_u(i) = 1.9d0;
        end
    end
    
    cn_w = zeros(np, 1);
    for i = 1:np
        %kii = zeros(1,nd);
        cn1_w = 0;
        cn2_w = 0;
        for j = 1:1
            if(abs(rotvelhalfold(i,j))>1e-15)
                cn1_w = cn1_w - rot(i,j)*rot(i,j)*(pcouple(i,j)/massvec2-pcoupleold(i,j)/massvec2) / (dt * rotvelhalfold(i,j));
            end
            cn2_w = cn2_w + rot(i,j)*rot(i,j);
        end
        if(abs(cn2_w)>1e-12)
            if ((cn1_w / cn2_w) > 0.0)
                cn_w(i) = 2.0 * sqrt(cn1_w / cn2_w);
            else
                cn_w(i) = 0.0;
            end
        else
            cn_w(i) = 0.0;
        end
        if (cn_w(i) > 2.0d0)
            cn_w(i) = 1.9d0;
        end
    end
    
    
    if (tt==1)
        velhalf = 1.0d0 * dt / massvec1 * (pforce + bforce) / 2.0;
        rotvelhalf = 1.0d0 * dt / massvec2 * (pcouple) / 2.0;
    else
        for i = 1:nd
            velhalf(:,i) = ((2.0d0 - cn_u * dt) .* velhalfold(:,i) + 2.0 * dt / massvec1 * (pforce(:,i) + bforce(:,i))) ./ (2.0 + cn_u * dt);
        end
        for i = 1:1
            rotvelhalf(:,i) = ((2.0d0 - cn_w * dt) .* rotvelhalfold + 2.0 * dt / massvec2 * (pcouple(:,i))) ./ (2.0 + cn_w * dt);
        end
    end
    % 时间积分
    %vel = 0.5d0 * (velhalfold + velhalf);
    dis = dis + velhalf * dt;
    velhalfold = velhalf;
    
    
    %rotvel = 0.5d0 * (rotvelhalfold + rotvelhalf);
    rot = rot + rotvelhalf * dt;
    rotvelhalfold = rotvelhalf;
    
    
    
    
    if(vel_y * currentTime<0.01*sideLengthY)
        dis(upNode,2)   = 0 * currentTime;
        dis(upNode,1)   = vel_y * currentTime;
        dis(downNode,2) = 0;
        dis(downNode,1) = 0;
        rot(upNode,1)   = 0;
        rot(downNode,1) = 0;
    else
        dis(upNode,2)   = 0;
        dis(upNode,1)   = 0.01*sideLengthY;
        dis(downNode,2) = 0;
        dis(downNode,1) = 0;
        rot(upNode,1)   = 0;
        rot(downNode,1) = 0;
    end
    
    fprintf('tt=%d\n',tt);
    rot(~innerNode) = 0;
    
    ddis = dis-dis_old;
    drot = rot-rot_old;
    dis_old = dis;
    rot_old = rot;
    %     dpforce = pforce - pforceold;
    %     dpcouple = pcouple - pcoupleold;
    U_L2(tt) = sqrt(sum(sum(ddis.*ddis)))/sqrt(sum(sum(dis.*dis)));
    W_L2(tt) = sqrt(sum(sum(drot.*drot)))/sqrt(sum(sum(rot.*rot)));
    F_L2(tt) = sqrt(sum(sum(pforce(innerNode,:).^2)))/sqrt(sum(sum(pforce.^2)));
    M_L2(tt) = sqrt(sum(sum(pcouple(innerNode,:).^2)))/sqrt(sum(sum(pcouple.^2)));
    
    %     if(mod(tt,100)==0)
    %         figure(400); clf
    %         scatter(coord(:,1), coord(:,2), 20, dis(:,2), 'filled')
    %         colormap(jet)
    %         colorbar
    %         axis equal
    %         title('u_2')
    %         pause(0.001)
    %         figure(500); clf
    %         dmg = accumarray(ih,fail) ./ accumarray(ih,true(size(fail)));
    %         scatter(coord(:,1), coord(:,2), 20, dis(:,1), 'filled')
    %         colormap(jet)
    %         colorbar
    %         axis equal
    %         title('u_1')
    %         pause(0.001)
    %         figure(600); clf
    %         scatter(coord(:,1), coord(:,2), 20, rot(:,1), 'filled')
    %         colormap(jet)
    %         colorbar
    %         title('\omega_z')
    %         axis equal
    %         pause(0.001)
    %         figure(700); clf
    %         plot(1:1:nt,log10(U_L2(1:1:nt)), 'r-')
    %         hold on
    %         plot(1:1:nt,log10(W_L2(1:1:nt)), 'b-')
    %         xlabel('Iterations'); ylabel('L2')
    %         figure(800); clf
    %         plot(1:1:nt,log10(F_L2(1:1:nt)), 'r-')
    %         hold on
    %         plot(1:1:nt,log10(M_L2(1:1:nt)), 'b-')
    %         xlabel('Iterations'); ylabel('L2')
    %         pause(0.001)
    %     end
    if(tt>=10 && U_L2(tt)<Utol && W_L2(tt)<Utol && F_L2(tt)<Ftol && M_L2(tt)<Ftol)
        U_L2 = U_L2(1:tt);
        W_L2 = W_L2(1:tt);
        F_L2 = F_L2(1:tt);
        M_L2 = M_L2(1:tt);
        break;
    end
end
X = reshape(coord(:,1), ndiv+1, ndiv+1);
Y = reshape(coord(:,2), ndiv+1, ndiv+1);
VOL = reshape(particleVolume, ndiv+1, ndiv+1);
U1 = reshape(dis(:,1), ndiv+1, ndiv+1);
U2 = reshape(dis(:,2), ndiv+1, ndiv+1);
W3 = reshape(rot(:,1), ndiv+1, ndiv+1);
end