function [X, Y, VOL, U1, U2, W3] = MP_Static_Model_II(mvalue, ndiv)

addpath('..\..\..\SharedFiles');

sideLengthX = 1;
sideLengthY = 1;
ndivx = ndiv;
ndivy = ndiv;
dx = sideLengthX/ndivx;
dy = sideLengthY/ndivy;
maxdx = max(dx, dy);
x = linspace(0, sideLengthX, ndivx+1);
y = linspace(0, sideLengthY, ndivy+1);

[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 = vecnorm(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);
g1_ = repmat(wfvj, 1, size(g1,2)) .* g1;
g2 = secondOrderGhatfun2D(xi, wfvj, ap);
g2_ = repmat(wfvj, 1, size(g2,2)) .* g2;

emod = 1;
nu   = 0.2;
% nu = nu/(1+nu);
% emod = emod*(1+nu^2);
lambda = nu*emod/(1-2*nu)/(1+nu);
mu     = emod/(2*(1+nu));
% couplingNumber = sqrt(0.5);
couplingNumber = 0.5;
mu_r   = couplingNumber^2/(1-couplingNumber^2) * mu;
characteristicLength = sideLengthY*1;
gamma  = 2 * mu * characteristicLength^2;
%microInertia = 1;

g11 = (mu+mu_r)*(g2_(:,1)+g2_(:,2)) + (mu-mu_r)*(2*g2_(:,1));
g12 = + (mu-mu_r)*g2_(:,3);
g21 = + (mu-mu_r)*g2_(:,3);
g22 = (mu+mu_r)*(g2_(:,1)+g2_(:,2)) + (mu-mu_r)*(2*g2_(:,2));

gq1 = (g11.*xi(:,2)-g12.*xi(:,1))*0.5;
gq2 = (g21.*xi(:,2)-g22.*xi(:,1))*0.5;
qtgq = (xi(:,2).*gq1 - xi(:,1).*gq2)*0.5;
% thetai = 8*ones(size(ih));
%fff = sparse(ih,1,thetai.*qtgq)*8;
g10 = g1_(:,1);
g01 = g1_(:,2);

gtheta = (g10.*xi(:,2)-g01.*xi(:,1))*0.5;
ggamma = gamma*(g2_(:,1)+g2_(:,2));


B = + sparse(ih, 3*jh-2, g10, np, 3*np) + sparse(ih, 3*jh-1, g01, np, 3*np) ...
    - sparse(ih, 3*ih-2, g10, np, 3*np) - sparse(ih, 3*ih-1, g01, np, 3*np)...
    + sparse(ih, 3*jh-0, gtheta, np, 3*np)...
    + sparse(ih, 3*ih-0, gtheta, np, 3*np); % 仍是加号原因在于平均转角

% 先中间，后两边，先上后下
D = + sparse(3*jh-2, 3*jh-2, g11, 3*np, 3*np) + sparse(3*jh-2, 3*jh-1, g12, 3*np, 3*np) ...
    + sparse(3*jh-1, 3*jh-2, g21, 3*np, 3*np) + sparse(3*jh-1, 3*jh-1, g22, 3*np, 3*np) ...
    - sparse(3*ih-2, 3*jh-2, g11, 3*np, 3*np) - sparse(3*ih-2, 3*jh-1, g12, 3*np, 3*np) ...
    - sparse(3*ih-1, 3*jh-2, g21, 3*np, 3*np) - sparse(3*ih-1, 3*jh-1, g22, 3*np, 3*np) ...
    - sparse(3*jh-2, 3*ih-2, g11, 3*np, 3*np) - sparse(3*jh-2, 3*ih-1, g12, 3*np, 3*np) ...
    - sparse(3*jh-1, 3*ih-2, g21, 3*np, 3*np) - sparse(3*jh-1, 3*ih-1, g22, 3*np, 3*np) ...
    + sparse(3*ih-2, 3*ih-2, g11, 3*np, 3*np) + sparse(3*ih-2, 3*ih-1, g12, 3*np, 3*np) ...
    + sparse(3*ih-1, 3*ih-2, g21, 3*np, 3*np) + sparse(3*ih-1, 3*ih-1, g22, 3*np, 3*np) ...
    ...
    + sparse(3*jh-2, 3*jh, gq1, 3*np, 3*np) + sparse(3*jh-1, 3*jh, gq2, 3*np, 3*np) ...
    + sparse(3*jh-2, 3*ih, gq1, 3*np, 3*np) + sparse(3*jh-1, 3*ih, gq2, 3*np, 3*np) ...
    - sparse(3*ih-2, 3*jh, gq1, 3*np, 3*np) - sparse(3*ih-1, 3*jh, gq2, 3*np, 3*np) ...
    - sparse(3*ih-2, 3*ih, gq1, 3*np, 3*np) - sparse(3*ih-1, 3*ih, gq2, 3*np, 3*np) ...
    ...
    + sparse(3*jh, 3*jh-2, gq1, 3*np, 3*np) + sparse(3*jh, 3*jh-1, gq2, 3*np, 3*np) ...
    - sparse(3*jh, 3*ih-2, gq1, 3*np, 3*np) - sparse(3*jh, 3*ih-1, gq2, 3*np, 3*np) ...
    + sparse(3*ih, 3*jh-2, gq1, 3*np, 3*np) + sparse(3*ih, 3*jh-1, gq2, 3*np, 3*np) ...
    - sparse(3*ih, 3*ih-2, gq1, 3*np, 3*np) - sparse(3*ih, 3*ih-1, gq2, 3*np, 3*np) ...
    ...
    + sparse(3*jh, 3*jh, qtgq, 3*np, 3*np) + sparse(3*jh, 3*ih, qtgq, 3*np, 3*np) ...
    + sparse(3*ih, 3*jh, qtgq, 3*np, 3*np) + sparse(3*ih, 3*ih, qtgq, 3*np, 3*np) ...
    ...
    + sparse(3*jh, 3*jh, ggamma, 3*np, 3*np) - sparse(3*jh, 3*ih, ggamma, 3*np, 3*np) ...
    - sparse(3*ih, 3*jh, ggamma, 3*np, 3*np) + sparse(3*ih, 3*ih, ggamma, 3*np, 3*np);


K = (lambda-mu+mu_r)*(B'*B) + D;
meps      = dx / 100;
upNode    = find(coord(:,2)-1 > -meps);
downNode  = find(coord(:,2) < meps);
leftNode  = find(coord(:,1)<0+meps & coord(:,2)<1-meps & coord(:,2)>meps);
rightNode = find(coord(:,1)>1-meps & coord(:,2)<1-meps & coord(:,2)>meps);
G = [...
    sparse((1:numel(upNode)),   3*upNode-2,   ones(size(upNode)),    numel(upNode),   3*np); ...
    sparse((1:numel(upNode)),   3*upNode-1,   ones(size(upNode)),    numel(upNode),   3*np); ...
    sparse((1:numel(upNode)),   3*upNode ,    ones(size(upNode)),    numel(upNode),   3*np); ...
    sparse((1:numel(downNode)), 3*downNode-2, ones(size(downNode)),  numel(downNode), 3*np);...
    sparse((1:numel(downNode)), 3*downNode-1, ones(size(downNode)),  numel(downNode), 3*np);...
    sparse((1:numel(downNode)), 3*downNode,   ones(size(downNode)),  numel(downNode), 3*np)];
F = zeros(3*np+size(G,1),1);
F(3*np+1:1:3*np+numel(upNode),1) = 0.01;
U = [K,G';G,sparse(size(G,1),size(G,1))]\F;
U1 = U(1:3:3*np);
U2 = U(2:3:3*np);
Omega3 = U(3:3:3*np);
X = x;
Y = y;
U1 = reshape(U1,size(x));
U2 = reshape(U2,size(x));
W3 = reshape(Omega3,size(x));
VOL = reshape(particleVolume,size(x));
rmpath('..\..\..\SharedFiles');
end

