function v = mesh_volume(model)
% Returns vector v, the volume for each element

elem = drop_IDs(model);
v = zeros(size(elem,1),1);
utypes = unique(elem(:,2));

all_types = [12, 23, 24, 34, 36, 38];
all_orders = [1, 1, 3, 1, 3, 3];

for i_type = 1 : length(all_types)
    type = all_types(i_type);
    order = all_orders(i_type);
    if ismember(type, utypes)
        sel = elem(:, 2) == type;
        [w, ~, ~, ~, e] = fe_quad(model.Nodes, elem(sel, :), type, order);
        v(sel) = accumarray(e, w);
    end
end

end