clear;

R = 1;
phi = pi/3;
nR = 100;

mesh = create_spherical_cap(R, phi, nR);

[~, ~, w, elem] = geo2gauss(mesh, 5);

W = accumarray(elem, w);

V_fe = sum(W);
V_an = 2*pi*R^2*(1 - cos(phi));

C2N = cell2node_interp(mesh);

elem = drop_IDs(mesh);

v = zeros(size(mesh.Nodes, 1),1);
for iElem = 1:size(elem,1)
    type = elem(iElem,2);
    nodenum = mod(type,10);
    nodinds = elem(iElem,4+(1:nodenum));
    v(nodinds) = v(nodinds) + 1/nodenum*W(iElem);
end