I have a Poisson equation with mixed boundary conditions:
\begin{equation} \begin{alignedat}{3} -\Delta u(\vec{x}) &= f(\vec{x}), &&\quad \vec{x}\in\Omega \setminus \mathcal{K} \\ \vec{n} \cdot \nabla u(\vec{x}) &= 0, &&\quad \vec{x} \in \partial \Omega \\ u(\vec{x}) &= g(\vec{x}), &&\quad \vec{x} \in \mathcal{K}. \end{alignedat} \end{equation}
The caveat is that the set $\mathcal{K}$ is made up of a collection of points: $\mathcal{K} = \ \{y_1,\ldots, y_m\}$. This is covered for example in "Multivariate Interpolation at Arbitrary Points Made Simple" by Meinguet. When I discretise this using the finite element method I have no issue as long as I have $y_1, \ldots, y_m$ be vertices in my mesh, i.e. I do not even need to approximate the boundary conditions. In the finite volume method one uses the divergence formulation instead which results in:
\begin{align} -\int_{\partial \mathcal{V}}\vec{n} \cdot \nabla u \,dS = -\int_{\mathcal{V}}\nabla \cdot (\nabla u) \,dV &= \int_{\mathcal{V}} f \,dV, \,\forall \mathcal{V} \subset \Omega. \end{align}
A standard finite volume scheme is the two point flux approximation scheme (TPFA) which requires constructing a Delaunay mesh and its corresponding Voronoi diagram. Then it is trivial to prescribe boundary conditions on some of the edges of the triangulation. However I have Dirichlet boundary conditions prescribed at the triangle vertices which have measure zero with respect to $dS$. Is there a clean way to handle this?
