The rise of metasurfaces to manipulate the polarization states of light motivates the development of versatile numerical methods able to model and analyze their polarimetric properties. Here we make use of a scattered-field formulation well suited to the Finite Element Method (FEM) to compute the Stokes-Mueller matrix of metasurfaces. The major advantage of the FEM lies in its versatility and its ability to compute the optical properties of structures with arbitrary and realistic shapes, and rounded edges and corners. We benefit from this method to design achiral, pseudo-chiral, and chiral metasurfaces with specific polarimetric properties. We compute and analyze their Mueller matrices. The accuracy of this method is assessed for both dielectric and metallic scatterers hosting Mie and plasmonic resonances.