%
% GBMM1D - ONE DIMENSIONAL Galerkin-based MESHLESS PROGRAM FOR SOLVING A 1D BAR OF UNIT LENGTH
%          SUBJECTED TO A LINEAR BODY FORCE OF MAGNITUDE X WHOSE EXACT SOLUTION IS GIVEN BY
%          u = (x/2 - x^3/6)/E
%
%          BACKGROUND CELL QUADRATURE IS EMPLOYED TO EVALUATE INTEGRALS
%              - CELLS ARE COINCIDE WITH THE INTERVALS BETWEEN THE NODES
%              - ONE POINT GAUSS QUADRATURE
%
%          LAGRANGIAN MULTIPLIER METHOD IS EMPLOYER TO IMPOSE THE ESSENTIAL BOUNDARY CONDITIONS
%

clear all

% SET UP NODAL COORDINATES ALONG BAR, DETERMINE NUMBER OF CELLS
dx = 0.1;                % Distance between adjacent nodes
xi = [0.0 : dx : 1.0];   % Nodal coordinates
nnodes = length(xi);
ncells = nnodes-1;

% SET MATERIAL PROPERITES
E = 1.0;     % Elastic modulus
area = 1.0;  % Area of cross section

% DETERMINE RADIUS OF SUPPORTS FOR EACH NODE
scale =30.0;
dm = scale*dx*ones(1,nnodes);

%SET UP GAUSS POINTS, WEIGHTS, AND JACOBIAN FOR EACH CELL
gg = zeros(1,ncells);  % Coordinates of Gauss points
jac = dx/2;            % Jacobian for each cell
weight = 2;            % weight for each Gauss points
gg = dx/2 : dx : 1.0-dx/2;

% INITIALIZE MATRICES
K = zeros(nnodes);
P = zeros(nnodes,1);
G = zeros(nnodes,1);

% LOOP OVER GAUSS POINTS
for j = 1:length(gg)
    xg = gg(j);

    % EVALUATE SHAPE FUNCTIONS AND THEIR DERIVATIVES AT GAUSS POINT xg
    [PHI, DPHI, DDPHI] = MLS1DShape(2, nnodes, xi, 1, xg, dm, 'SPLIN',0.0);

    % ASSEMBLE DISCRETE EQUATIONS
    K = K + (weight*E*area*jac)*(DPHI'*DPHI);

    fbody = area*xg;
    P = P + (weight*fbody*jac)*PHI';
end

% ENFORCE BOUNDARY CONDITION USING LAGRANGE MULTIPLIERS

xg = 0.0;  % Prescribed displacement boundary
[PHI, DPHI, DDPHI] = MLS1DShape(2, nnodes, xi, 1, xg, dm, 'SPLIN',0.0);

G(1:3,1) = -PHI(1:3)';
Q = [0];
M = [K G; G' zeros(1)];

% SOLVE FOR NODAL PARAMETERS
d  = M\[P' Q]';

uh = zeros(nnodes,1);  % Nodal displacements
sh = zeros(nnodes,1);  % Nodal stress

for j=1:nnodes
   [PHI, DPHI, DDPHI] = MLS1DShape(2, nnodes, xi, 1, xi(j), dm, 'SPLIN',0.0);
   uh(j) = PHI * d(1:nnodes);
   sh(j) = E * DPHI * d(1:nnodes);
end

% EVALUATE RELATIVE ERROR NORMS
ue = (xi/2.0 - xi.*xi.*xi/6.0)/E;  % Exact solution
se = (1 - xi.*xi)/2.0;

erru = norm(ue'-uh)/norm(ue)*100
errs = norm(se'-sh)/norm(se)*100

% PLOT RESULTS
figure
subplot(1,2,1);  plot(xi, ue,'-*', xi, uh,'-o');
ylabel('位移');
legend('精確解','GBMM解',2);
subplot(1,2,2);  plot(xi, se,'-*', xi, sh,'-o');
ylabel('應力');
legend('精確解','GBMM解');

% Output nodal displacements and stresses
fid1 = fopen('G1DBarDis.dat','w');
fid2 = fopen('G1DBarStr.dat','w');

fprintf(fid1,'%10s%10s%10s\n', 'x', 'ue','uh');
fprintf(fid2,'%10s%10s%10s\n', 'x', 'se','sh');

for j = 1 : nnodes
   fprintf(fid1,'%10.4f%10.4f%10.4f\n', xi(j), ue(j), uh(j));
   fprintf(fid2,'%10.4f%10.4f%10.4f\n', xi(j), se(j), sh(j));
end
   
fclose(fid1);
fclose(fid2);