We consider minimization problems involving the Dirichlet integral under an arbitrary number of volume constraints on the level sets and a generalized boundary condition. More precisely, given a bounded open domain Ω ⊂ ℝⁿ with smooth boundary, we study the problem of minimizing f_Ω |∇u|² among all those functions u ∈ H¹ that simultaneously satisfy n-dimensional measure constraints on the level sets of the kind |{u = l_i}| = α₁, i = 1,..., k, and a generalized boundary condition u ∈. Here, Κ is a closed convex subset of H¹ such that Κ + H₀¹ = Κ; the invariance of Κ under H₀¹ provides that the condition u ∈ Κ actually depends only on the trace of u along δΩ. By a penalization approach, we prove the existence of minimizers and their Hölder continuity, generalizing previous results that are not applicable when a boundary condition is prescribed. Finally, in the case of just two volume constraints, we investigate the Γ-convergence of the above (rescaled) functionals when the total measure of the two prescribed level sets tends to saturate the whole domain Ω. It turns out that the resulting Γ-limit functional can be split into two distinct parts: the perimeter of the interface due to the Dirichlet energy that concentrates along the jump, and a boundary integral term due to the constraint u ∈ Κ. In the particular case where Κ = H¹ (i.e. when no boundary condition is prescribed), the boundary term vanishes and we recover a previous result due to Ambrosio et al.