function [faces, eind, cent, norm, r2] = get_face_properties(model)

elements = drop_IDs(model);
faces = get_faces(elements);

% Create a mesh of the faces
face_mesh = struct( ...
    'Nodes', model.Nodes, ...
    'Elements', [ ...
        (1 : size(faces, 1)).', faces(:, 2), ones(size(faces, 1), 2), ...
        faces(:, 3:end), zeros(size(faces, 1), 10-size(faces, 2)) ...
    ], ...
    'Materials', [1 1 0 0 0 0], ...
    'Properties', [1 1 0 0 0 0] ...
);


[cent, norm] = centnorm(face_mesh);

if nargout > 4
% Compute face radii
r2 = -inf(size(faces, 1), 1);

u_types = unique(faces(:, 2));
for i_type = 1 : length(u_types)
    type = u_types(i_type);
    n_nod = mod(type, 10);
    sel = faces(:, 2) == type;
    for i_nod = 1 : n_nod
        d = model.Nodes(faces(sel, 2 + i_nod), 2 : 4) - cent(sel, :);
        r2(sel) = max(r2(sel), dot(d, d, 2));
    end
end

end

% Sort to element indices and faces
eind = faces(:, 1);
faces = faces(:, 2:end);


end