If you have the Jacobians from joint velocities to link body velocities (you can calculate these by taking each link $i$ as if it were the end effector and calculating its body Jacobian $J^{b}_{i}$), then you can get the mass matrix for each link with respect to the joint angles as the transpose of the link's Jacobian multiplied by the link's mass and inertia matrix $mu_{i}$ multiplied by the link's Jacobian,
$$M_{i} = (J^{b}_{i})^{T} \mu_{i} J^{b}_{i},$$
and the total mass matrix with respect to the joint velocities as the sum of the individual links' mass matrices,
$$M = \sum_{i} (J^{b}_{i})^{T} \mu_{i} J^{b}_{i}.$$
Assuming that we put the link body frames at their centers of mass, the link mass matrices have the form
$$ \mu_{i} = \begin{bmatrix} m_{i} & & \\ & m_{i} & \\ && I_{i} \end{bmatrix}$$
for planar systems (where $m_{i}$ is the mass of the link and $I_{i}$ is its moment of inertia), or
$$ \mu_{i} = \begin{bmatrix} m_{i} & & &&& \\ & m_{i} & &&& \\ &&m_{i} &&& \\ &&& \phantom{-}I_{xx,i} & -I_{xy,i} & -I_{xz,i}\\&&& -I_{yx,i} & \phantom{-}I_{yy,i} & -I_{yz,i} \\&&& -I_{zx,i} & -I_{zy,i} & \phantom{-}I_{zz,i} \end{bmatrix}$$
for systems moving in three-dimensional space. If the link body frames are aligned with the principle axes (longest and shortest mass distributions) of the links, then the off-diagonal $I$ terms in the local expression of the three-dimensional moment of inertia become zero,
$$ \mu_{i} = \begin{bmatrix} m_{i} & & &&& \\ & m_{i} & &&& \\ &&m_{i} &&& \\ &&& I_{xx,i} & & \\&&& & I_{yy,i} & \\&&& & & I_{zz,i} \end{bmatrix}.$$