Hello yc,
If you are willing to make an approximation to the area of each grid quadrilateral then I think there is a reasonably viable solution. Consider a square with corner points
a b
c d
and small x and y displacements ax, ay of the top left corner, similarly for b,c,d. The area of the resulting quadrilateral is a complicated function of the 8 deflections but if you are willing to go with the area to lowest order in the deflections and drop higher order terms such as ax*by and so forth, which is reasonable if the deflections are on the order of .1, then the change in area is simply (1/2)* [(bx-ax) + (dx-cx) + (ay -cy) + (by-dy)] and for each quadrilateral to not change in area from its original value, then the sum of all those terms is 0. If you have n small squares per side, then there are (n+1)^2 points, 2*(n+1)^2 deflections and n^2 constraints on the quadrilateral areas. That leaves 2*(n+1)^2 - n^2 random degrees of freedom. I think the following code gets it done. I have not tried to optimise it for speed.
% mesh (n+1)x(n+1) of n squares per side c d.
% 1,1 index is in upper left corner, row index is y, column index is x.
% After constraints are taken into account, the x and y displacements
% of the corners of the squares are derived from uniform random variables
% with -g < x,y < g.
n = 10;
g = .2;
m = zeros(n^2,2*(n+1)^2); % constraint matrix, operates on a column vector of x deflections
% concatenated with a column vector of y deflections
for j = 1:n
for k = 1:n
mx = zeros(n+1,n+1);
mx(j:j+1,k:k+1) = [-1 1;-1 1];
my = zeros(n+1,n+1);
my(j:j+1,k:k+1) = [1 1; -1 -1];
m(j+(k-1)*n,1:(n+1)^2) = mx(:)';
m(j+(k-1)*n,(n+1)^2+1:2*(n+1)^2) = my(:)';
end
end
nullm = null(m);
df = 2*(n+1)^2-n^2; % number of degrees of freedom
r = g*(2*rand(df,1)-1); % random numbers about 0
delxy = nullm*r;
delx = delxy(1:(n+1)^2);
dely = delxy((n+1)^2+1:2*(n+1)^2);
delx = reshape(delx,n+1,n+1);
dely = reshape(dely,n+1,n+1);
% <new calculation>
% check that Alin (the linear approx to A) = 1 for all quadrilaterals.
% find the fractional difference of Alin compared to the true area, A
% D = (Alin-A)/A
% for this calc go to complex variables. visually, +1 is to the right, +i is down
u = delx-i*dely;
da = u(2:end,2:end)-u(1:end-1,1:end-1);
bc = u(1:end-1,2:end)-u(2:end,1:end-1);
A = (1/2)*imag((da+1+i).*conj(bc+1-i));
% multiply this out, discard quadratic terms
Alin = 1 + (1/2)*imag((da+conj(bc))*(1+i));
D = Alin./A -1;
max(D,[],'all')
min(D,[],'all')
