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!

[Aloe]

Hi @Vinay, it says error when I run. Arrays have incompatible sites for this operation.

[Walter Roberson]

@Aloe

You could get that error if numCols or numRows are vectors or arrays (non-scalar). In such a case, Lx or Ly would be vectors or arrays, and the Lx * x_offset(:) or Ly * y_offset(:) could easily fail.

[Aloe]

@Walter Roberson and plus they are added together so it wouldn't work. Will changing the sizes of Lx work to make them compatible?

[Walter Roberson]

The code requires that numCols and numRows are scalars. Or that they are 1 x 9 vectors. Or that they are (9 times size(x,1)) by 9 arrays.

The code was probably designed for the case where x and y are scalars.

[Aloe]

so I can't do anything?

[Walter Roberson]

x_repm = repmat(x, 9, 1) + repmat(Lx * x_offset(:), size(x,1), 1);
y_repm = repmat(y, 9, 1) + repmat(Ly * y_offset(:), size(y,1), 1);

[Vinay]

Hi @Aloe,

I've included the offset case in the original code. The MATLAB file is attached.

[Aloe]

Thanks @Vinay, that works. Malo