Problem description: I have a 2D smooth manifold with boundary $M$ embedded in $\mathbb R^3$, discretized by a triangular mesh $T$, and need to find the eigenfunctions of the Laplace-Beltrami operator $\Delta$ on $M$ (the Riemannian metric of $M$ is induced from its $\mathbb R^3$ embedding). Since $M$ is a manifold with boundary, I need the solution to be: $$ \Delta f=\lambda f\ \ \textrm{in}\ \ M,\qquad f=0\ \ \textrm{on}\ \ \partial M. $$
Numerical formulation: For a numerical solution, I apply the finite element method to the triangular mesh $T$. Let $V$ be the vector denoting the values of $f$ at the vertices of $T$. Then $\Delta$ can be discretized by the cotan-Laplacian matrix $L$ (see e.g. here), and a mass matrix $M$ (a diagonal matrix denoting the area of each vertex), i.e. $\Delta \approx M^{-1}L$. The eigenvalue problem (not considering boundary conditions) is then: $$ M^{-1}LV=\lambda V\iff LV=\lambda MV.$$ Since $L$ and $M$ are both symmetric and positive semi-definite, this amounts to a generalized eigenvalue problem, which can be solved with existing numerical linear algebra algorithm/packages.
My question: Is there a (standard?) numerical formulation of the eigenvalue problem that incorporates the desired boundary conditions, which can be solved with existing linear algebra packages? Naively the boundary condition leads to an additional constraint: $$ BV=0 $$ where $B$ is a diagonal matrix, with the $i$-th diag element is $1$ if $V_i$ is a boundary vertex and $0$ otherwise. However, there doesn't seem to be an algorithm that allows incorporating such additional constraints.
More generally, is it possible to use other boundary conditions, e.g., Dirichlet, Neumann, or these two mixed?