% MWLS1D - MESHLESS WEIGHTED LEAST-SQUARE PROGRAM TO SOLVE 1D ELASTIC PROBLEM - A 1D BAR OF UNIT 
%          LENGTH SUBJECTED TO BODY FORCE LINEARLY VARIATING WITH THE COORDINATE. THE LEFT END IS  
%          FIXED AND THE RIGHT IS FREE.
% THE ANALYTICAL SOLUTION IS GIVEN BY
%           u = (x/2 - x^3/6)/E
% WHERE X REPRESENTS THE COORDINATE ALONG THE BAR AND E IS THE ELASTIC MODULUS
% THE FUNCTIONAL IS CALCULATED BY NODAL SUNMMATION, WHICH AVOIDS NUMERICAL INTEGRATION.

clear all

% SET MATERIAL PROPERTIES
E = 1.0;          % ELASTIC MODULUS
area = 1.0;       % THE AREA OF THE CROSS SECTION
L = 1.0           % THE LENGTH OF THE BAR

% SET UP NODE DISTRIBUTION
dx = 0.1;                   % DISTANCE BETWEEN ADJACENT NODES
xi = [0.0 : dx : L];        % NODAL COORDINATES
nnodes = length(xi);        % NUMBER OF NODES

% SET COMPUTATIONAL PARAMETERS, INCLUDING THE DIMENSION OF SUPPORT AND THE VALUE OF PENALTY.
scale = 2.5;                       % SCALE FACTOR
dm = scale*dx*ones(1,nnodes);      % SUPPORT RADIUS
tracpen = 1e5;                     % PENALTY OF TRACTION BOUNDARY
disppen = tracpen*(E/L)^2;         % PENALTY OF DISPLACEMENT BOUNDARY

% MATRICE INITIALIZATION
K = zeros(nnodes,nnodes);           % STIFFNESS MATRIX
P = zeros(nnodes,1);                % LOAD VECTOR

% EVALUATE SHAPE FUNCTION AND ITS DERIVATIVES AT ALL NODES
[phi,dphi,ddphi] = MLS1DShape(3, nnodes, xi, nnodes, xi, dm, 'SPLIN', 0.0);     % SHAPE FUNCTIONS AT NODES

% LOOP OVER ALL NODES IN DOMAIN
for m = 1:nnodes
    K = K + E^2 * ddphi(m,:)' * ddphi(m,:);       % ASSEMBLED TO STIFFNESS MATRIX
    fbody = xi(m) * area;
    P = P - E * ddphi(m,:)' * fbody;              % ASSEMBLED TO LOAD VECTOR    
end

% ENFORCEMENT OF DISPLACEMENT BOUDNARY CONDITIONS
K = K + disppen * phi(1,:)'*phi(1,:);

% ENFORCEMENT OF TRACTION BOUNDARY CONDITIONS. Note: THOUGH THE TRACTION BOUNDARY HERE IS FREE, ITS
% CONTRIBUTION TO STIFFNESS MATRIX SHOULD ALSO BE CONSIDERED
K = K + tracpen*E^2 * dphi(nnodes,:)'*dphi(nnodes,:);  

% SOLVE FOR NODAL PARAMETERS
nodp = K\P;                               % NODAL FICTIOUS VALUE

% CALCULATE NODAL DISPLACEMENTS AND STRESSES
uh = phi*nodp;                            % APPROXIMATE VALUE OF DISPLACEMENTS 
sh = E*dphi*nodp;                         % APPROXIMATE VALUE OF STRESSES

% EVALUATE RELATIVE ERRORS
ue = (xi/2.0 - xi.*xi.*xi/6.0)/E;         % EXACT DISPLACEMENTS
se = (1.0 - xi.*xi)/2.0;                  % EXACT STRESSES

erru = norm(ue'-uh)/norm(ue)*100    % RELATIVE ERROR OF DISPLACEMENTS
errs = norm(se'-sh)/norm(se)*100    % RELATIVE ERROR OF STRESSES

% DRAW RESULT CURVE
figure(1);
subplot(1,2,1);
hu = plot(xi,ue,xi,uh);
xlabel('Coordinate','Fontsize',12);
ylabel('Displacement','Fontsize',12);
legend(hu,'Exact Solution','WLSM Solution');
grid on;

subplot(1,2,2);
hs = plot(xi,se,xi,sh);
xlabel('Coordinate','Fontsize',12);
ylabel('Stress','Fontsize',12);
legend(hs,'Exact Solution','WLSM Solution');
grid on;

% Output nodal displacements and stresses
fid1 = fopen('1DBarDis.dat','w');
fid2 = fopen('1DBarStr.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);