/sci/ - Science & Math

here's another one, "an implicit meshing algorithm".Please excuse the lack of vectorization. Also, I am aware that all the forces are being computed twice, I have versions that run faster and on any 2D boundary, this is just the compact version. Enjoy the animation!

clear all, close all; clc
figure('units','normalized','position',[.2 .2 .6 .6])

r = .5; % Boundary radius
k = -.5; % Spring constant
h = .5; % 'Time' step
len = .08; % Spring element rest length
N = 200; % Number of nodes
A = rand(N,2); % Generating N random x-y points [0 to 1]
A = (A*(r*2) - r)/2; % Scaling and shifting points
for it = 1:75
fx = zeros(1,N);
fy = zeros(1,N);
Tri = delaunay(A(:,1),(A(:,2)));
for n = 1:N
for i = 1:length(Tri)
for j = 1:3
if Tri(i,j)==n
l = [1,2,3,1,2];
lenx1 = (A(Tri(i,l(j+1)),1) - A(Tri(i,j),1));
lenx2 = (A(Tri(i,l(j+2)),1) - A(Tri(i,j),1));
leny1 = (A(Tri(i,l(j+1)),2) - A(Tri(i,j),2));
leny2 = (A(Tri(i,l(j+2)),2) - A(Tri(i,j),2));
nrm1 = (lenx1^2 + leny1^2)^.5;
nrm2 = (lenx2^2 + leny2^2)^.5;
fx(n) = fx(n)+k*((len-nrm1)*(lenx1/nrm1) + (len-nrm2)*(lenx2/nrm2));
fy(n) = fy(n)+k*((len-nrm1)*(leny1/nrm1) + (len-nrm2)*(leny2/nrm2));
end
end
end
A(n,1)= A(n,1)+h*fx(n);
A(n,2)= A(n,2)+h*fy(n);
end
for n = 1:N
NRM = (A(n,1)^2 + A(n,2)^2)^.5;
if NRM > r
A(n,:) = A(n,:)*r/NRM;
end
end
triplot(Tri,A(:,1),(A(:,2))), hold on
[x,y,z] = cylinder(r,200);
plot(x(1,:),y(1,:),'r')
axis equal, hold off
title('Implicit Meshing of a Circle')
getframe;
end