diff --git a/lapy/diffgeo.py b/lapy/diffgeo.py
index 06e0990..31be8be 100644
--- a/lapy/diffgeo.py
+++ b/lapy/diffgeo.py
@@ -146,21 +146,21 @@ def compute_geodesic_f(
 
     gradf = compute_gradient(geom, vfunc)
     fem = Solver(geom, lump=True, use_cholmod=use_cholmod)
-    fem.mass = sparse.eye(fem.stiffness.shape[0], dtype=fem.stiffness.dtype)
+    # fem.mass = sparse.eye(fem.stiffness.shape[0], dtype=fem.stiffness.dtype)
 
     if scalar_input:
         # gradf: (n_elements, 3)
         gradnorm = gradf / np.sqrt((gradf**2).sum(1))[:, np.newaxis]
         gradnorm = np.nan_to_num(gradnorm)
         divf = compute_divergence(geom, gradnorm)
-        vf = fem.poisson(divf)
+        vf = fem.poisson(divf, integrate=False)
         vf -= vf.min()
     else:
         # gradf: (n_elements, n_functions, 3) — norm along last axis
         gradnorm = gradf / np.sqrt((gradf**2).sum(-1))[:, :, np.newaxis]
         gradnorm = np.nan_to_num(gradnorm)
         divf = compute_divergence(geom, gradnorm)  # (n_vertices, n_functions)
-        vf = fem.poisson(divf)                     # (n_vertices, n_functions)
+        vf = fem.poisson(divf, integrate=False)     # (n_vertices, n_functions)
         vf -= vf.min(axis=0)
     return vf
 
@@ -200,21 +200,21 @@ def tria_compute_geodesic_f(
     fem = Solver(tria, lump=True, use_cholmod=use_cholmod)
     # div is the integrated divergence (so it is already B*div);
     # pass identity instead of B here
-    fem.mass = sparse.eye(fem.stiffness.shape[0])
+    # fem.mass = sparse.eye(fem.stiffness.shape[0])
 
     if scalar_input:
         # gradf: (n_triangles, 3)
         gradnorm = gradf / np.sqrt((gradf**2).sum(1))[:, np.newaxis]
         gradnorm = np.nan_to_num(gradnorm)
         divf = tria_compute_divergence(tria, gradnorm)
-        vf = fem.poisson(divf)
+        vf = fem.poisson(divf, integrate=False)
         vf -= vf.min()
     else:
         # gradf: (n_triangles, n_functions, 3) — norm along last axis
         gradnorm = gradf / np.sqrt((gradf**2).sum(-1))[:, :, np.newaxis]
         gradnorm = np.nan_to_num(gradnorm)
         divf = tria_compute_divergence(tria, gradnorm)  # (n_vertices, n_functions)
-        vf = fem.poisson(divf)                          # (n_vertices, n_functions)
+        vf = fem.poisson(divf, integrate=False)          # (n_vertices, n_functions)
         vf -= vf.min(axis=0)
     return vf
 
@@ -503,20 +503,20 @@ def tria_compute_rotated_f(
     gradf = tria_compute_gradient(tria, vfunc)
     tn = tria.tria_normals()
     fem = Solver(tria, lump=True, use_cholmod=use_cholmod)
-    fem.mass = sparse.eye(fem.stiffness.shape[0], dtype=vfunc.dtype)
+    # fem.mass = sparse.eye(fem.stiffness.shape[0], dtype=vfunc.dtype)
     dtup = (np.array([0]), np.array([0.0]))
 
     if scalar_input:
         # gradf: (n_triangles, 3)
         gradf = np.cross(tn, gradf)
         divf = tria_compute_divergence(tria, gradf)
-        vf = fem.poisson(divf, dtup)
+        vf = fem.poisson(divf, dtup, integrate=False)
     else:
         # gradf: (n_triangles, n_functions, 3)
         # tn: (n_triangles, 3) -> broadcast via (n_triangles, 1, 3)
         gradf = np.cross(tn[:, np.newaxis, :], gradf)  # (n_triangles, n_functions, 3)
         divf = tria_compute_divergence(tria, gradf)    # (n_vertices, n_functions)
-        vf = fem.poisson(divf, dtup)                   # (n_vertices, n_functions)
+        vf = fem.poisson(divf, dtup, integrate=False)   # (n_vertices, n_functions)
     return vf
 
 
diff --git a/lapy/solver.py b/lapy/solver.py
index daa61ca..e5a45ea 100644
--- a/lapy/solver.py
+++ b/lapy/solver.py
@@ -784,6 +784,7 @@ def poisson(
         h: float | np.ndarray = 0.0,
         dtup: tuple = (),
         ntup: tuple = (),
+        integrate: bool = True,
     ) -> np.ndarray:
         """Solver for the Poisson equation with boundary conditions.
 
@@ -808,6 +809,11 @@ def poisson(
             Neumann boundary condition as a tuple containing the index and
             data arrays of same length. The default, an empty tuple,
             corresponds to Neumann on all boundaries.
+        integrate: bool, default=True
+            Whether to integrate the right hand side over the surface using 
+            the mass matrix. If True, the right hand side is effectively 
+            replaced by ``B h``. If False, the right hand side is used as is, 
+            which corresponds to ``A x = h``. 
 
         Returns
         -------
@@ -904,8 +910,9 @@ def poisson(
                 (dim, n_rhs), dtype=dtype
             )
         # compute right hand side
-        mass = self.mass.astype(dtype, copy=False)
-        b = mass * (h - nvec)
+        b = h - nvec
+        if integrate:
+            b = self.mass.astype(dtype, copy=False) * b
         if len(didx) > 0:
             b = b - self.stiffness * dvec
         # remove Dirichlet Nodes