function [PHI, DPHI, DDPHI] = MLS1DShape(m, nnodes, xi, npoints, x, dm, wtype, para) % SHAPE FUNCTION OF 1D MLS APPROXIMATION % % SYNTAX: [PHI, DPHI, DDPHI] = MLS1DShape(m, nnodes, xi, npoints, x, dm, wtype, para) % % INPUT PARAMETERS % m - Total number of basis functions (1: Constant basis; 2: Linear basis; 3: Quadratic basis) % nnodes - Total number of nodes used to construct MLS approximation % npoints - Total number of points whose MLS shape function to be evaluated % xi(nnodes) - Coordinates of nodes used to construct MLS approximation % x(npoints) - Coordinates of points whose MLS shape function to be evaluated % dm(nnodes) - Radius of support of nodes % wtype - Type of weight function % para - Weight function parameter % % OUTPUT PARAMETERS % PHI - MLS Shpae function % DPHI - First order derivatives of MLS Shpae function % DDPHI - Second order derivatives of MLS Shpae function % % INITIALIZE WEIGHT FUNCTION MATRICES wi = zeros (1, nnodes); % Weight funciton dwi = zeros (1, nnodes); ddwi = zeros (1, nnodes);
% LOOP OVER ALL EVALUATION POINTS TO CALCULATE VALUE OF SHAPE FUNCTION Fi(X) for j = 1 : npoints
% DETERMINE WEIGHT FUNCTIONS AND THEIR DERIVATIVES AT EVERY NODE for i = 1 : nnodes di = x(j) - xi(i); [wi(i), dwi(i), ddwi(i)] = Weight(wtype, para, di, dm(i)); end
% EVALUATE BASIS p, B MATRIX AND THEIR DERIVATIVES if (m == 1) % Shepard function p = [ones(1, nnodes)]; px = [1]; dpx = [0]; ddpx = [0];
B = p .* [wi]; DB = p .* [dwi]; DDB = p .* [ddwi]; elseif (m == 2) p = [ones(1, nnodes); xi]; px = [1; x(j)]; dpx = [0; 1]; ddpx = [0; 0];
B = p .* [wi; wi]; DB = p .* [dwi; dwi]; DDB = p .* [ddwi; ddwi]; elseif (m == 3) p = [ones(1, nnodes); xi; xi.*xi]; px = [1; x(j); x(j)*x(j)]; dpx = [0; 1; 2*x(j)]; ddpx = [0; 0; 2];
B = p .* [wi; wi; wi]; DB = p .* [dwi; dwi; dwi]; DDB = p .* [ddwi; ddwi; ddwi]; else error('Invalid order of basis.'); end
% EVALUATE MATRICES A AND ITS DERIVATIVES A = zeros (m, m); DA = zeros (m, m); DDA = zeros (m, m); for i = 1 : nnodes pp = p(:,i) * p(:,i)';
A = A + wi(i) * pp; DA = DA + dwi(i) * pp; DDA = DDA + ddwi(i) * pp; end