advanceBoundary.m
function [Xl] = advanceBoundary(Xl0, u, X, Y, dt, N, h);
[idxs delta] = evalDelta(Xl0, X, Y, N, h);
u1 = u(:,1);
u2 = u(:,2);
Xl = Xl0 + (dt*h*h) * [ dot(u1(idxs), delta, 2) dot(u2(idxs), delta, 2) ];
calcLagrangianForce.m
function [F] = calcLagrangianForce(Xl, L, dtheta, XT, K, Kp);
persistent Lap isInit L0Up L0Dn idxUp idxDn tmp LUp LDn TUp TDn numL;
if( size( isInit ) == 0 )
isInit = 1;
e = ones(L,1);
Lap = spdiags([e -2*e e], [-1:1], L, L);
% make periodic:
Lap(1,L) = 1;
Lap(L,1) = 1;
Lap = (K / (dtheta * dtheta)) * Lap;
numL = length(XT(:,1));
idxUp = [2:numL 1]';
idxDn = [numL 1:(numL-1)]';
tmp = XT(idxUp,:) - XT;
L0Up = sqrt(dot(tmp,tmp,2));
tmp = XT - XT(idxDn,:);
L0Dn = sqrt(dot(tmp,tmp,2));
LUp = zeros(numL,1);
LDn = zeros(numL,1);
TUp = zeros(numL,2);
TDn = zeros(numL,2);
end
% no rest length springs
F = [( Lap * Xl(:,1) ) ( Lap * Xl(:,2) )];
% springs with rest length:
%tmp = Xl(idxUp,:) - Xl;
%LUp = sqrt(dot(tmp,tmp,2));
%TUp = [(tmp(:,1) ./ LUp) (tmp(:,2) ./ LUp)];
%tmp = Xl - Xl(idxDn,:);
%LDn = sqrt(dot(tmp,tmp,2));
%TDn = [(tmp(:,1) ./ LDn) (tmp(:,2) ./ LDn)];
%F = zeros(numL, 2);
%F(:,1) = (K / (dtheta*dtheta)) * ( (LUp - L0Up) .* TUp(:,1) - (LDn - L0Dn) .* TDn(:,1) );
%F(:,2) = (K / (dtheta*dtheta)) * ( (LUp - L0Up) .* TUp(:,2) - (LDn - L0Dn) .* TDn(:,2) );
% target points
%F = F + Kp * ( XT - Xl );
D02DPeriodic.m
function [D0x, D0y] = D02DPeriodic(N, h);
M = N * N;
e = ones(M,1);
eu1 = ones(M,1);
eu1((N+1):N:M) = 0;
ed1 = ones(M,1);
ed1(N:N:M) = 0;
% Periodic 2D Centered Difference Operators on Square:
D0x = spdiags([-ed1 eu1], [-1 1], M, M);
D0y = spdiags([-e e], [-N N], M, M);
for( i = 1:N )
D0x( (i-1)*N + 1, i*N ) = -1;
D0x( i*N, (i-1)*N + 1 ) = 1;
D0y( i, N*N-N+i ) = -1;
D0y( N*N-N+i, i ) = 1;
end
D0x = D0x / (2*h);
D0y = D0y / (2*h);
D02DPeriodicFT.m
function [D0xT, D0yT] = D02DPeriodicFT(N, h);
M = N * N;
% 2D Fourier Space Indexes:
l = repmat( (0:(N-1))', N, 1 );
m = reshape( repmat( (0:(N-1)), N, 1 ), M, 1);
% Fourier Transformed 2D Centered Difference Operators on Square:
D0xT = (sqrt(-1) / h) * sin( (2 * pi / N) * l );
D0yT = (sqrt(-1) / h) * sin( (2 * pi / N) * m );
本文分享自 图像处理与模式识别研究所 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!