% Nodes generation. Only uniformly distributed nodes

% number of nodes along two directions
nnx = 5;
nny = 5;
numnode = nnx*nny ;

% domain dimension
L = 4 ;
W = 4 ;

% node generation
node = square_node_array([-L/2 -W/2],[L/2 -W/2],[L/2 W/2],[-L/2 W/2],nnx,nny);

% plot the nodes
colordef white
h = plot(node(:,1),node(:,2),'b.');
set(h,'MarkerSize',15);
hold on
plot([node(1,1) node(nnx,1)],[node(1,2) node(nnx,2)],'b-');
plot([node(nnx,1) node(nnx*nny,1)],[node(nnx,2) node(nnx*nny,2)],'b-');
plot([node(nnx*nny,1) node(nnx*nny-nnx+1,1)],[node(nnx*nny,2) node(nnx*nny-nnx+1,2)],'b-');
plot([node(nnx*nny-nnx+1,1) node(1,1)],[node(nnx*nny-nnx+1,2) node(1,2)],'b-');
axis off
axis equal

% plot circular domain of influence
di = 1.4 ; % radius of support
xc = node(13,1) ;
yc = node(13,2) ;
theta = -pi:0.1:pi ;
x  = xc + di * cos(theta) ;
y  = yc + di * sin(theta) ;
hc = plot(x,y,'r-') ;
set(hc,'LineWidth',3);

% highlight the central node at which weight function and shape function
% are plotted
hn = plot(node(13,1),node(13,2),'r*');
set(hn,'MarkerSize',15);

% compute weight function = spline of order 4, use refined node
% distribution to get a smooth plot
x = linspace(-2,2,8*nnx);
y = linspace(-2,2,8*nny);
[x,y] = meshgrid(x,y); 
w = zeros(8*nnx,8*nny);
for i = 1:8*nnx
    for j = 1:8*nny
        r = sqrt( (x(1,i) - xc)^2 + (y(j,1) - yc)^2 )/di ;
        if (r <= 1.0)
            w(i,j) = 1 - 6*r*r + 8*r*r*r - 3*r*r*r*r;
        else
            w(i,j) = 0.0 ;
        end
    end
end

% plot weight function as 3D plot
figure
surf(x,y,w);

% compute MLS shape function at central nodes
% linear intrinsic basis: p = [1 x y]

% initialization
phi = zeros(8*nnx,8*nny);
dphidx = zeros(8*nnx,8*nny);
dphidy = zeros(8*nnx,8*nny);

form  = 'quartic_spline' ;
nodeI = [xc yc];  % node at which shape function is to be evaluated
for i = 1:8*nnx
    for j = 1:8*nny
       pt = [x(1,i) y(j,1)]; 
       [index] = define_support(node,pt,di);
       [Phi,phidx,phidy] = MLS_ShapeFunction(pt,nodeI,index,node,di,form);
       phi(i,j) = Phi ;
       dphidx(i,j) = phidx ;
       dphidy(i,j) = phidy ;
    end
end

% plot shape function
figure
surf(x,y,phi);
figure
surf(x,y,dphidx);
figure
surf(x,y,dphidy);

nnx= 41;
x = linspace(0,10,nnx);
xa = 5;
rmax = 3;
phi = zeros(nnx);
for i = 1:nnx
    r    = abs(x(i) - xa)/rmax ;
    drdx = sign(x(i) - xa)/rmax ;
    if (r <= 1.0)
        phi(i) = -1/6*r^3 +0.5*r*r -0.5*r + 1/6;
        dphidr = -0.5*r^2 +r -0.5;
    else
        phi(i) = 0.0 ;
        dphidr = 0.0 ;
    end
    dphi(i) = dphidr*drdx;
end

figure
hold on
plot(x,phi,'b-')
plot(x,dphi,'r.-')

figure
hold on
plot(x,phi,'b-')
plot(x,dphi,'r.-')

% Plot branch function
x = linspace(0,10,40);
y = linspace(0,10,40);
[x,y] = meshgrid(x,y); 
phi1 = zeros(40,40);
phi2 = zeros(40,40);
phi3 = zeros(40,40);
phi4 = zeros(40,40);
dphi1dx = zeros(40,40);
for i = 1:40
    for j = 1:40        
        pt = [x(1,i) y(j,1)];
        xp=QT*(pt-xTip)';                  % local coordinates
        [theta,r]=cart2pol(xp(1),xp(2));  % local polar coordinates
        [f,dfdx,dfdy] = branch(r,theta,alpha);
        phi1(i,j) = f(1); dphi1dx(i,j) = dfdx(1);
        phi2(i,j) = f(2);
        phi3(i,j) = f(3);
        phi4(i,j) = f(4);
    end
end

% plot weight function as 3D plot
figure
surf(x,y,phi1);
title('B_1 function')
opts = struct('Color','rgb','Bounds','tight');
exportfig(gcf,'B1.eps',opts)

figure
surf(x,y,phi2);
title('B_2 function')
exportfig(gcf,'B2.eps',opts)

figure
surf(x,y,phi3);
title('B_3 function')
exportfig(gcf,'B3.eps',opts)

figure
surf(x,y,phi4);
title('B_4 function')
exportfig(gcf,'B4.eps',opts)

figure
surf(x,y,dphi1dx);
title('B_1 function')
opts = struct('Color','rgb','Bounds','tight');
exportfig(gcf,'B1.eps',opts)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 4 point quadrature for each "1D element"
[W1,Q1]=quadrature(4, 'GAUSS', 1 );

Wt = [];
Qt = [];
Jt = [];

for i = 1 : (length(topNodes) - 1)
    sctr = [topNodes(i) topNodes(i+1)];
    for q = 1:size(W1,1)
        pt = Q1(q,:);
        wt = W1(q);
        [N,dNdxi] = lagrange_basis('L2',pt);
        J0=dNdxi'*node(sctr,:);
        Qt = [Qt; N' * node(sctr,:)];
        Wt = [Wt; wt];
        Jt = [Jt; norm(J0)];
    end % of quadrature loop
end

for igp = 1 : size(Wt,1)
    pt = Qt(igp,:);                             % quadrature point
    wt = Wt(igp);                               % quadrature weight
    [index] = define_support(node,pt,di);
    en = zeros(1,2*size(index,2));
    force = zeros(1,2*size(index,2));
    [phi,dphidx,dphidy] = MLS_ShapeFunction(pt,index,node,di,form);

    % prescribed traction vector t_bar
    tx = 0 ;
    ty = 1000;
    for j = 1 : size(index,2)
        en(2*j-1) = u_pos(index(j))   ;
        en(2*j  ) = u_pos(index(j))+1 ;
        force(2*j-1) = tx*phi(j);
        force(2*j  ) = ty*phi(j);
    end
    f(en) = f(en) + Jt(igp) * wt * force' ;
end