cbmm1d.m

自己整理出来的几个常用的一维无网格方法的matlab程序。

%
% CBMM1D - ONE DIMENSIONAL Collocation-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
%

clear all

% SET UP NODAL COORDINATES ALONG BAR
dx = 0.1;                % Distance between adjacent nodes
xi = [0.0 : dx : 1.0];   % Nodal coordinates
nnodes = length(xi); 

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

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

K = zeros(nnodes);    % Coefficient matrix
P = zeros(nnodes,1);  % Right-hand side vector

[PHI, DPHI, DDPHI] = MLS1DShape(2, nnodes, xi, nnodes, xi, dm, 'WNEW',0.0);

K(1,:) = PHI(1,:);    % The 1st equation is the prescribed displacement condition at
P(1)   = 0.0;         % the left end of the bar

for i = 2:nnodes-1              % The 2nd to (nnodes-1)th equations are the equilibrium
   K(i,:) = E * DDPHI(i,:);     % equations at all inner nodes
   P(i)   = -xi(i);
end

K(nnodes,:) = DPHI(nnodes,:);   % The last equation is the prescribed traction condition
P(nnodes)   = 0.0;              % at the right end of the bar

d = K \ P;

uh = PHI * d;        % Nodal displacements
sh = E * DPHI * d;   % Nodal stress
   
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('精确解','CBMM解',2);
subplot(1,2,2);  plot(xi, se,'-*', xi, sh,'-o');
ylabel('应力');
legend('精确解','CBMM解');
% Output nodal displacements and stresses
fid1 = fopen('C1DBarDis.dat','w');
fid2 = fopen('C1DBarStr.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);