正方形上的周期二维拉普拉斯算子。

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) ];

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 );

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);

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 );