% MLPG1D - A PROGRAM FOR SOLVING 1D ELASTOSTATIC PROBLEM: A BAR OF UNIT LENGTH SUBJECTED TO
%          LINEAR BODY FORCE. THE LEFT END OF THE BAR IS FIXED AND THE RIGHT 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;
%          IN THE COMPUTATION, THE INTEGRAL IS CALCULATED OVER THE SUBDOMAIN; THUS THE BACKGROUND
%          CELLS FOR INTEGRATION ARE NO LONGER REQUIRED.

%          THE PENALTY METHOD IS EMPLOYED TO ENFORCE THE ESSENTIAL BOUNDARY CONDITIONS

clear all

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

% SET UP NODAL COORDINATES
deltax = 0.1;           % DISTANCE BETWEEN ADJACENT NODES
xi = 0.0:deltax:1.0;    % NODAL COORDINATES
nnode = length(xi);     % NUMBER OF NODES

% SET UP SUPPORT RADIUS
spscale = 2.5;    %\ SCALE FACTOR OF SUPPORT
spradius = spscale*deltax*ones(1,nnode);    % SUPPORT RADIUS

% SET UP INTERGRATION SUBDOMAIN RADIUS
sdscale = 0.6;   % SCALE FACTOR OF SUBDOMAIN
sdradius = sdscale*deltax*ones(1,nnode);    % SUBDOMAIN RADIUS

% MATRICE INITIALIZING
K = zeros(nnode);     % STIFFNESS MATRIX
P = zeros(nnode,1);   % FORCE VECTOR
numcellgp = 1;                  % NUMBER OF GAUSS POINTS IN ONE CELL
gsspst = zeros(numcellgp,1);    % POSITION OF GAUSS POINTS
gsswgh = zeros(numcellgp,1);    % WEIGHT OF GAUSS POINTS

switch numcellgp                % INITIALIZE GAUSS POINTS
case 1
    gsswgh(1) = 2.0;
    gsspst(1) = 0.0;
case 2
    gsswgh(1) = 1.0;
    gsswgh(2) = 1.0;
    gsspst(1) = -0.577350269189626;
    gsspst(2) =  0.577350269189626;
case 3
    gsswgh(1) =  0.555555555555556;
    gsswgh(2) =  0.888888888888889;
    gsswgh(3) =  0.555555555555556;
    gsspst(1) = -0.774596669241483;
    gsspst(2) =  0.0;
    gsspst(3) =  0.774596669241483;
end

alpha = 1e5;     % PENALTY

% LOOP OVER EACH SUBDOMAIN
for i = 1:nnode
    sdleft = xi(i)-sdradius(i);        % LEFT BOUNDARY OF SUBDOMAIN
    sdright = xi(i)+sdradius(i);       % RIGHT BOUNDARY OF SUBDOMAIN
    ncell = 8;                         % NUMBER OF INTEGRATION CELL IN A SUBDOMAIN
    celldiam = 2*sdradius(i)/ncell;    % DIAMETER OF EACH INTEGRATION CELL
    jacobi = celldiam/2;               % JACOBI DETERMINANT OF A CELL

    % LOOP OVER GAUSS POINTS
    for j = 1:ncell
        cellleft = sdleft + (j-1)*celldiam;  % LEFT BOUNDARY OF INTEGRATION CELL
        cellright = cellleft + celldiam;     % RIGHT BOUNDARY OF INTEGRATION CELL

        if cellleft>0.0 & cellright<1.0     % INTEGRATION CELL IN THE DOMAIN
            coef1 = (cellright + cellleft)/2;
            coef2 = (cellright - cellleft)/2;

            for m = 1:numcellgp
                xg = coef1 + coef2*gsspst(m);  % GET COORDINATE OF GAUSS POINT

                if xg>0.0 & xg<1.0            % POINT XG IN THE DOMAIN
                    % VALUE OF TEST FUNCTION AT POINT XG
                    [wi,dwi,ddwi] = Weight('SPLIN', 0.0, xg-xi(i), sdradius(i));

                    % VALUE OF SHAPE FUNCTION AT POINT XG
                    [phi,dphi,ddphi] = MLS1DShape(2, nnode, xi, 1, xg, spradius, 'SPLIN', 0.0);

                    K(i,:) = K(i,:) + dwi*E*area*dphi*gsswgh(m)*jacobi;  % ADDED TO STIFFNESS MATRIX

                    fbody = xg*area;                           % BODY FORCE
                    P(i) = P(i) + wi*fbody*gsswgh(m)*jacobi;   % ADDED TO LOAD VECTOR
                end
            end
        end
    end

    % ENFORCEMENT OF ESSENTIAL BOUNDARY CONDITION

    if sdleft <= 0.0   % THE LEFT BOUNDARY OF A SUBDOMAIN IS THE COUNTERPART OF THE WHOLE DOMAIN
        [wi,dwi,ddwi] = Weight('SPLIN', 0.0, -xi(i), sdradius(i));
        [phi,dphi,ddphi] = MLS1DShape(2, nnode, xi, 1, 0.0, spradius, 'SPLIN', 0.0);

        % PENALTY METHOD USED HERE TO ENFORCE ESSENTIAL BOUNDARY CONDITIONS
        K(i,:) = K(i,:) + wi*E*area*dphi;
        K(i,:) = K(i,:) + alpha*wi*area*phi;
    end
end

% SOLVE LINEAR ALGEBRAIC EQUATIONS FOR FICTITIOUS NODAL DISPLACEMENTS

nodp = K\P;   % FICTITIOUS NODAL DISPLACEMENTS

% CALCULATE NODAL DISPLACEMENTS AND STRESSES

% SHAPE FUNCTIONS AT NODES
clear phi dphi ddphi;
[phi,dphi,ddphi] = MLS1DShape(2, nnode, xi, nnode, xi, spradius, 'SPLIN', 0.0);

dsph = phi*nodp;       % APPROXIMATE VALUE OF DISPLACEMENTS
strh = E*dphi*nodp;    % APPROXIMATE VALUE OF STRESSES

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

errdsp = norm(dspe'-dsph)/norm(dspe)*100    % RELATIVE ERROR OF DISPLACEMENTS
errstr = norm(stre'-strh)/norm(stre)*100    % RELATIVE ERROR OF STRESSES

% DRAW RESULT CURVE
figure(1); subplot(1,2,1); hu = plot(xi,dspe,'k-',xi,dsph,'ko');
xlabel('Coordinate','Fontsize',12);
ylabel('Displacement','Fontsize',12); legend(hu,'Exact Solution','MLPG Solution'); grid on;

subplot(1,2,2); hs = plot(xi,stre,'k-',xi,strh,'ko');
xlabel('Coordinate','Fontsize',12);
ylabel('Stress','Fontsize',12); legend(hs,'Exact Solution','MLPG Solution');