Computing gradient of mean curvature on a mesh using gp toolbox, unexpected error. What am I missing?

I am trying to compute this quantity here H grad(H) where H is the mean curvture

The code breaks at the last line because there is a dimension mismatch between H and grad(H). size(grad(H))= [960 3] and size(V)=[162 3] and size(mean_curvture)=[162 3]. Why is the grad(H) has that dimension? where is the mistake?

% Load a mesh
[V, F] = load_mesh(‘sphere.off’);

% Compute Mean Curvature and Normal Vector

% 2. Compute Cotangent Laplacian (L) and Mass Matrix (M)
L = cotmatrix(V, F);
M = massmatrix(V, F, ‘barycentric’);

% 3. Compute Mean Curvature Vector (H)
H = -inv(M) * (L * V);

% 4. Compute the Magnitude of Mean Curvature (mean curvature at each vertex)
mean_curvature = sqrt(sum(H.^2, 2));

% Compute the Normal Vector
normals = per_vertex_normals(V, F);

% Compute Gaussian Curvature
gaussian_curvature = discrete_gaussian_curvature(V,F);

% Compute Laplacian of Mean Curvature
laplacian_H = L * mean_curvature;

% Compute the Gradient of Mean Curvature
% Use gptoolbox’s gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;

% Compute the new term: newthing
%H_squared = mean_curvature .^ 2;
%newthing = laplacian_H + ((H_squared – 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;
% Compute the new term: newthing
H_squared = mean_curvature .^ 2;
newthing = laplacian_H + ((H_squared – 2 * gaussian_curvature) .* mean_curvature) .* normals + mean_curvature .* grad_H;

% Define small arc equation
h = 0.1; % Small parameter h
vertex = V; % V contains the vertices

% Compute the small arc
f_h = vertex + (mean_curvature .* normals) * h + newthing * (h^2 / 2);

% Visualization of small arc
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), ‘FaceColor’, ‘cyan’, ‘EdgeColor’, ‘none’);
hold on;
plot3(f_h(:,1), f_h(:,2), f_h(:,3), ‘r.’, ‘MarkerSize’, 10);
axis equal;
title(‘Small Arc Visualization’);
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘Z’);
the grad(H) script
% Step 1: Load a Mesh
[V, F] = load_mesh(‘sphere.off’);

% Step 2: Compute Laplace-Beltrami Operator
L = cotmatrix(V, F); % cotangent Laplace-Beltrami operator

% Step 3: Compute the Mean Curvature Vector (H)
H = -L * V; % Mean curvature normal vector

% Step 4: Compute the Mean Curvature Magnitude
mean_curvature = sqrt(sum(H.^2, 2));

% Step 5: Compute the Gradient of Mean Curvature
% Use gptoolbox’s gradient operator for scalar fields
grad_H = grad(V, F) * mean_curvature;

% Step 6: Visualization or Further Processing
figure;
trisurf(F, V(:,1), V(:,2), V(:,3), mean_curvature, ‘EdgeColor’, ‘none’);
axis equal;
lighting gouraud;
camlight;
colorbar;
title(‘Mean Curvature Gradient’);
xlabel(‘X’);
ylabel(‘Y’);
zlabel(‘Z’);
grad(V, F) from gp toolbox
function [G] = grad(V,F)
% GRAD Compute the numerical gradient operator for triangle or tet meshes
%
% G = grad(V,F)
%
% Inputs:
% V #vertices by dim list of mesh vertex positions
% F #faces by simplex-size list of mesh face indices
% Outputs:
% G #faces*dim by #V Gradient operator
%
% Example:
% L = cotmatrix(V,F);
% G = grad(V,F);
% dblA = doublearea(V,F);
% GMG = -G’*repdiag(diag(sparse(dblA)/2),size(V,2))*G;
%
% % Columns of W are scalar fields
% G = grad(V,F);
% % Compute gradient magnitude for each column in W
% GM = squeeze(sqrt(sum(reshape(G*W,size(F,1),size(V,2),size(W,2)).^2,2)));
%

dim = size(V,2);
ss = size(F,2);
switch ss
case 2
% Edge lengths
len = normrow(V(F(:,2),:)-V(F(:,1),:));
% Gradient is just staggered grid finite difference
G = sparse(repmat(1:size(F,1),2,1)’,F,[1 -1]./len, size(F,1),size(V,1));
case 3
% append with 0s for convenience
if size(V,2) == 2
V = [V zeros(size(V,1),1)];
end

% Gradient of a scalar function defined on piecewise linear elements (mesh)
% is constant on each triangle i,j,k:
% grad(Xijk) = (Xj-Xi) * (Vi – Vk)^R90 / 2A + (Xk-Xi) * (Vj – Vi)^R90 / 2A
% grad(Xijk) = Xj * (Vi – Vk)^R90 / 2A + Xk * (Vj – Vi)^R90 / 2A +
% -Xi * (Vi – Vk)^R90 / 2A – Xi * (Vj – Vi)^R90 / 2A
% where Xi is the scalar value at vertex i, Vi is the 3D position of vertex
% i, and A is the area of triangle (i,j,k). ^R90 represent a rotation of
% 90 degrees
%
% renaming indices of vertices of triangles for convenience
i1 = F(:,1); i2 = F(:,2); i3 = F(:,3);
% #F x 3 matrices of triangle edge vectors, named after opposite vertices
v32 = V(i3,:) – V(i2,:); v13 = V(i1,:) – V(i3,:); v21 = V(i2,:) – V(i1,:);

% area of parallelogram is twice area of triangle
% area of parallelogram is || v1 x v2 ||
n = cross(v32,v13,2);

% This does correct l2 norm of rows, so that it contains #F list of twice
% triangle areas
dblA = normrow(n);

% now normalize normals to get unit normals
u = normalizerow(n);

% rotate each vector 90 degrees around normal
%eperp21 = bsxfun(@times,normalizerow(cross(u,v21)),normrow(v21)./dblA);
%eperp13 = bsxfun(@times,normalizerow(cross(u,v13)),normrow(v13)./dblA);
eperp21 = bsxfun(@times,cross(u,v21),1./dblA);
eperp13 = bsxfun(@times,cross(u,v13),1./dblA);

%g = …
% ( …
% repmat(X(F(:,2)) – X(F(:,1)),1,3).*eperp13 + …
% repmat(X(F(:,3)) – X(F(:,1)),1,3).*eperp21 …
% );

GI = …
[0*size(F,1)+repmat(1:size(F,1),1,4) …
1*size(F,1)+repmat(1:size(F,1),1,4) …
2*size(F,1)+repmat(1:size(F,1),1,4)]’;
GJ = repmat([F(:,2);F(:,1);F(:,3);F(:,1)],3,1);
GV = [eperp13(:,1);-eperp13(:,1);eperp21(:,1);-eperp21(:,1); …
eperp13(:,2);-eperp13(:,2);eperp21(:,2);-eperp21(:,2); …
eperp13(:,3);-eperp13(:,3);eperp21(:,3);-eperp21(:,3)];
G = sparse(GI,GJ,GV, 3*size(F,1), size(V,1));

%% Alternatively
%%
%% f(x) is piecewise-linear function:
%%
%% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk
%% ∇f(x) = … = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk
%% = ∇φi fi + ∇φj fj + ∇φk) fk
%%
%% ∇φi = 1/hjk ((Vj-Vk)/||Vj-Vk||)^perp =
%% = ||Vj-Vk|| /(2 Aijk) * ((Vj-Vk)/||Vj-Vk||)^perp
%% = 1/(2 Aijk) * (Vj-Vk)^perp
%%
%m = size(F,1);
%eperp32 = bsxfun(@times,cross(u,v32),1./dblA);
%G = sparse( …
% [0*m + repmat(1:m,1,3) …
% 1*m + repmat(1:m,1,3) …
% 2*m + repmat(1:m,1,3)]’, …
% repmat([F(:,1);F(:,2);F(:,3)],3,1), …
% [eperp32(:,1);eperp13(:,1);eperp21(:,1); …
% eperp32(:,2);eperp13(:,2);eperp21(:,2); …
% eperp32(:,3);eperp13(:,3);eperp21(:,3)], …
% 3*m,size(V,1));

if dim == 2
G = G(1:(size(F,1)*dim),:);
end

% Should be the same as:
% g = …
% bsxfun(@times,X(F(:,1)),cross(u,v32)) + …
% bsxfun(@times,X(F(:,2)),cross(u,v13)) + …
% bsxfun(@times,X(F(:,3)),cross(u,v21));
% g = bsxfun(@rdivide,g,dblA);

case 4
% really dealing with tets
T = F;
% number of dimensions
assert(dim == 3);
% number of vertices
n = size(V,1);
% number of elements
m = size(T,1);
% simplex size
assert(size(T,2) == 4);

if m == 1 && ~isnumeric(V)
simple_volume = @(ad,r) -sum(ad.*r,2)./6;
simple_volume = @(ad,bd,cd) simple_volume(ad, …
[bd(:,2).*cd(:,3)-bd(:,3).*cd(:,2), …
bd(:,3).*cd(:,1)-bd(:,1).*cd(:,3), …
bd(:,1).*cd(:,2)-bd(:,2).*cd(:,1)]);
P = sym(‘P’,[1 3]);
V1 = V(T(:,1),:);
V2 = V(T(:,2),:);
V3 = V(T(:,3),:);
V4 = V(T(:,4),:);
V1P = V1-P;
V2P = V2-P;
V3P = V3-P;
V4P = V4-P;
A1 = simple_volume(V2P,V4P,V3P);
A2 = simple_volume(V1P,V3P,V4P);
A3 = simple_volume(V1P,V4P,V2P);
A4 = simple_volume(V1P,V2P,V3P);
B = [A1 A2 A3 A4]./(A1+A2+A3+A4);
G = simplify(jacobian(B,P).’);
return
end

% f(x) is piecewise-linear function:
%
% f(x) = ∑ φi(x) fi, f(x ∈ T) = φi(x) fi + φj(x) fj + φk(x) fk + φl(x) fl
% ∇f(x) = … = ∇φi(x) fi + ∇φj(x) fj + ∇φk(x) fk + ∇φl(x) fl
% = ∇φi fi + ∇φj fj + ∇φk fk + ∇φl fl
%
% ∇φi = 1/hjk = Ajkl / 3V * (Facejkl)^perp
% = Ajkl / 3V * (Vj-Vk)x(Vl-Vk)
% = Ajkl / 3V * Njkl / ||Njkl||
%

% get all faces
F = [ …
T(:,1) T(:,2) T(:,3); …
T(:,1) T(:,3) T(:,4); …
T(:,1) T(:,4) T(:,2); …
T(:,2) T(:,4) T(:,3)];
% compute areas of each face
A = doublearea(V,F)/2;
N = normalizerow(normals(V,F));

% compute volume of each tet
vol = volume(V,T);

GI = …
[0*m + repmat(1:m,1,4) …
1*m + repmat(1:m,1,4) …
2*m + repmat(1:m,1,4)];
GJ = repmat([T(:,4);T(:,2);T(:,3);T(:,1)],3,1);
GV = repmat(A./(3*repmat(vol,4,1)),3,1).*N(:);
G = sparse(GI,GJ,GV, 3*m,n);

end
end
the sphere.off file (mesh data)I am trying to compute this quantity here H grad(H) where H is the mean curvture