addpath('FEMDL/src/2d'); init

t0 = 0; tf = 1;                                     % initial and final time
v = @double_gyre;                                   % vector field
T = @(x) at_tf(flowmap(v, x, [t0 tf]));             % flow map
DL = @(DT) 0.5*(eye(2) + inv(DT)*inv(DT)');         % dynamic Laplacian
DLx = @(x) fapply1(DL, D(T,x));                     % evaluate DL at each row of x
p = grid2(50,50);                                   % nodes on a grid
m = trimesh(p);                                     % triangulation
A = triquad(m,DLx,2);                               % integrate DL on triangles
[K,M] = assemble2(m,A);                             % assemble stiffness and mass matrix
[V,L] = eigs(K,M,6,'SM');                           % solve eigenproblem
[lam,ord] = sort(diag(L),'descend'); V = V(:,ord);  % sort eigenvalues
figure; plotf2(m,normed(V(:,3))); colorbar;         % plot eigenvector
W = kmeans(V(:,1:3),3,'Replicates',20);             % kmeans clustering
figure; scatter(p(:,1),p(:,2),30,W,'filled');       % plot clustering
Tp = T(p);                                          % advect grid
figure; scatter(Tp(:,1),Tp(:,2),30,W,'filled');     % plot advected clustering
colormap(jet(nc));