Error using repmat to apply periodic boundary conditions to voronoi diagram.

Hello everyone

I have been trying to apply periodic boudnary conditions to a voronoi diagram using the built-in command 'repmat'. I want the classic primary domain with its repeated neighbours surrounding it like a 3-by-3 matrix but I do not why it is not performing well. The results show that the points are dividing themselves at the early timesteps before they look like the hexagonal shape that I want but the boundaries do not look like they are "wrapping around" like what a periodic boundary condition is supposed to do. Please, if anyone knows anything, I would appreciate it. I have been stuck on this problem for 2 weeks now. I have pasted the code below.

Thanks.

numRows = 4;
numCols = 4;

x = [];
y = [];

for row = 0:numRows-1
    for col = 0:numCols-1
        x_center = col * (sqrt(3)/2);
        y_center = row;
        
    if mod(col, 2) == 1
        y_center = y_center + 0.5;
    end
    
    x = [x; x_center];
    y = [y; y_center];
    end
end

for i = 1:length(x)
    noise_x = 0.1 - rand(1) * 0.2;
    noise_y = 0.1 - rand(1) * 0.2;
    x(i) = x(i) + noise_x;
    y(i) = y(i) + noise_y;
end

x = double(x(:));
y = double(y(:));

min_x = min(x);
max_x = max(x);
min_y = min(y);
max_y = max(y);

x_repm = repmat(x,3,3);
y_repm = repmat(y,3,3);

x_rep = reshape(x_repm, [],1);
y_rep = reshape(y_repm, [],1);

% Parameters
eta = 1.0;
deltat = 0.001;
T = 0.5;
mu_ij = 50;
s_ij = 1;

figure;

voronoi(x_rep,y_rep);
xlim([0 10]);
ylim([0 10]);

for iter = 0:deltat:T
    F = zeros(length(x_rep), 2);
    dt = delaunay(x_rep, y_rep);
    tri = triangulation(dt, x_rep, y_rep);
    
    for i = 1:length(x_rep)
        neighbours = unique(dt(dt(:,1) == i | dt(:,2) == i, :));
        unique_neighbours = unique(neighbours(:));
        
        for j = unique_neighbours(:)'
            if i ~= j
                rij_x = x_rep(j) - x_rep(i);
                rij_y = y_rep(j) - y_rep(i);
                dist = [rij_x, rij_y];
                r_ij = norm(dist);
                
                F_ij = mu_ij * (dist / r_ij) * (r_ij - s_ij);
                F(i,:) = F(i,:) + F_ij;
                F(j,:) = F(j,:) - F_ij;
            end
        end
    end
    
    x_new = x_rep + (deltat / eta) * F(:, 1);
    y_new = y_rep + (deltat / eta) * F(:, 2);
    
    x_rep = x_new;
    y_rep = y_new;
    
    clf;
    voronoi(x_rep, y_rep);
    xlim([min_x max_x]);
    ylim([min_y max_y]);
    
    hold on;
    title(sprintf('Time: %.2f', iter));
    dt = delaunay(x_rep, y_rep);
    triplot(dt, x_rep, y_rep, '--r');
    drawnow;
end
hold off;

Answer 1

[Vinay]

Hi @Aloe,

The voronoi diagram shows the intial 16 points stored as (x,y) coordinates in the code. The 'repmat' creates an 48*3 array of the 'x' and 'y' coordinates reshaped to 144*1 but the values of the coordinates are copy of the initial 16 coordinates causing overlapping of the points.

To resolve the issue, we have to shift the points in both direction by adding or subtracting the domain width and height to create the periodic effect

%  domain size defined as the x_center and y_center
Lx = (numCols - 1) * (sqrt(3)/2);
Ly = numRows;

% replicated points with offsets
[x_offset, y_offset] = meshgrid(-1:1, -1:1);
x_repm = repmat(x, 9, 1) + Lx * x_offset(:);
y_repm = repmat(y, 9, 1) + Ly * y_offset(:);

I hope this helps!