function M = AddMass_mult_rot(model, M, c_wall, c_fluid, c_air)
mesh = model.Mesh;
nodes = mesh.Nodes;
elements = mesh.Elements;

Nelems = size(elements,2);

regions = zeros(Nelems,1);
regions(findElements(mesh,'region','Face',1)) = 1; 
regions(findElements(mesh,'region','Face',2)) = 2;
regions(findElements(mesh,'region','Face',3)) = 3;

for e = 1:Nelems
    vert = elements(:,e);

    x = nodes(1,vert);   % r-coordinates
    y = nodes(2,vert);   % z-coordinates

    Ae = polyarea(x,y);  % element area
    rbar = mean(x);      % mean radius (axisymmetric weight)

    switch regions(e) % Select heat capacity
        case 1
            c = c_wall;
        case 2
            c = c_fluid;
        case 3
            c = c_air;
    end
    % Axisymmetric element mass matrix
    Me = rbar * c * Ae / 12 * [2 1 1; 1 2 1; 1 1 2];
    % Assemble
    M(vert,vert) = M(vert,vert) + Me;
end
end