% Generate the penalty-function matrix G and place it in a struct.

fem.xmesh = meshextend(fem);

% Extract Gauss-point coordinates and mesh-element numbers from fem:
me_td = posteval(fem,'meshelement','geomnum',2,'spoint',postgp('tri',1));
me_qd = posteval(fem,'meshelement','geomnum',2,'bpoint',postgp('quad',1));
  
p = [me_td.p, me_qd.p]; % first-order Gauss-point coordinates
me_numbers = [me_td.d, me_qd.d]; % mesh-element numbers 

n = length(me_numbers); % number of parameters

for i = 1:n
  pos(i) = find(me_numbers == i);
end

[x1,x2] = meshgrid(p(1,pos(:)),p(1,pos(:)));
xdist2 = (x2-x1).^2;
[y1,y2] = meshgrid(p(2,pos(:)),p(2,pos(:)));
ydist2 = (y2-y1).^2;
dist = (xdist2 + ydist2).^(1/2);

% Calculate the covariance matrix and its inverse:
sigma2 = 1; l = 50; % variogram sill and range parameters
Q = 2*sigma2*(1-exp(-dist/l));
Q_inv = inv(Q); Q_inv = 0.5*(Q_inv+Q_inv'); % restore symmetry

% Calculate G:
X = ones(n,1);
G = Q_inv - Q_inv*X*X'*Q_inv/(X'*Q_inv*X); G = 0.5*(G+G');

% Define structs for Q_inv and G:
rng = 1:n;
Q_inv_struct.x = rng; Q_inv_struct.y = rng; Q_inv_struct.data = Q_inv;
G_struct.x = rng; G_struct.y = rng; G_struct.data = G;