Skip to content

Solver.poisson support pre-integrated vfuncs #115

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
magnesium2400 wants to merge 1 commit into
base: main
Choose a base branch
from
+18 −11
Open

Solver.poisson support pre-integrated vfuncs #115

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
18 changes: 9 additions & 9 deletions lapy/diffgeo.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Copy link

Copilot AI Apr 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There’s now commented-out dead code (# fem.mass = sparse.eye(...)). Since integrate=False is the intended mechanism, consider removing the commented assignment and leaving a short explanatory comment (e.g., that compute_divergence already returns an integrated divergence so Poisson should not re-apply B).

Suggested change
# fem.mass = sparse.eye(fem.stiffness.shape[0], dtype=fem.stiffness.dtype)
# compute_divergence already returns an integrated divergence, so Poisson
# should solve with integrate=False and must not re-apply the mass matrix.

Copilot uses AI. Check for mistakes.

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

Expand Down Expand Up @@ -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])
Comment on lines 201 to +203
Copy link

Copilot AI Apr 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The comment here still describes the old workaround (“pass identity instead of B here”) and keeps the identity assignment commented out. Please update the comment to reflect the new approach (integrate=False) and remove the commented-out fem.mass = sparse.eye(...) line to avoid confusing future readers.

Copilot uses AI. Check for mistakes.

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

Expand Down Expand Up @@ -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)
Copy link

Copilot AI Apr 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

There’s leftover commented-out code (# fem.mass = sparse.eye(...)) from the pre-integrate workaround. Since this function now calls poisson(..., integrate=False), please remove the commented assignment and (if needed) replace it with a brief note explaining why integration is disabled here.

Suggested change
# fem.mass = sparse.eye(fem.stiffness.shape[0], dtype=vfunc.dtype)
# The divergence has already been assembled, so disable integration in the
# Poisson solve instead of using the old mass-matrix workaround.

Copilot uses AI. Check for mistakes.
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


Expand Down
11 changes: 9 additions & 2 deletions lapy/solver.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand All @@ -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``.
Comment on lines +812 to +816
Copy link

Copilot AI Apr 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new parameter doc header uses integrate: bool but the rest of the docstring follows NumPy style (name : type). Also, integrate=False changes the equation being solved and skips mass-matrix scaling for the Neumann term as well (since b = h - nvec is no longer multiplied by B). Please adjust the docstring (including the earlier text that implies A x = B h) and clarify the expected ntup convention when integrate=False.

Copilot uses AI. Check for mistakes.

Returns
-------
Expand Down Expand Up @@ -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
Comment on lines 912 to +915
Copy link

Copilot AI Apr 8, 2026

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new integrate flag introduces a second execution path for RHS construction. Please add unit tests covering integrate=False (including at least one case with Dirichlet/Neumann terms if supported) and asserting it matches the previous workaround (e.g., setting mass to identity or pre-multiplying h appropriately), so behavior stays stable across refactors.

Copilot uses AI. Check for mistakes.
if len(didx) > 0:
b = b - self.stiffness * dvec
# remove Dirichlet Nodes
Expand Down
Loading