% 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;