|
clear, clc, clf
Lx = 5;
Nx = 25;
Ly = 0.5;
Ny = 5;
t = 0.1;
Slimit = 1000;
%CREATE NODES
[X,Y]=meshgrid(0:Lx/Nx:Lx,0:Ly/Ny:Ly);
X=X’;
Y=Y’;
nodes = [X(1:end)’ Y(1:end)’];
%BLOCK 1
N = 1:Nx; n1a = N; n1b = N+1; %horizontal –
n2a = [N Nx+1]; n2b = n2a+Nx+1; %vertical ||
n3a = N; n3a(2:2:end) = n3a(2:2:end)+1; %diagonal //
n3b = N+Nx+1; n3b(1:2:end) = n3b(1:2:end)+1; %diagonal \\
block1 = [n1a n2a n3a; n1b n2b n3b]’;
%BLOCK 2
N = 1:Nx; n1a = N; n1b = N+1; %horizontal —
n2a = [N Nx+1]; n2b = n2a+Nx+1; %vertical ||
n3a = N; n3a(1:2:end) = n3a(1:2:end)+1; %diagonal \\
n3b = N+Nx+1; n3b(2:2:end) = n3b(2:2:end)+1; %diagonal //
block2 = [n1a n2a n3a; n1b n2b n3b]’;
%COMBINE BLOCKS
elem = [];
for i = 1:Ny
if rem(i,2) == 1
elem = cat(1,elem,block1+(Nx+1)*(i-1));
else
elem = cat(1,elem,block2+(Nx+1)*(i-1));
end
end
elem = [elem ; [N’ N’+1]+Ny*(Nx+1)] %top_elements
%PLOTTING
clf
hold on
for i=1:length(elem)
x=nodes(elem(i,:),1);
y=nodes(elem(i,:),2);
plot(x,y,’b.-‘)
%pause(0.1)
end
plot(nodes(:,1),nodes(:,2),’ro’)
hold off
axis equal
|
