lapy.Solver

class lapy.Solver(geometry, lump=False, aniso=None, aniso_smooth=10, use_cholmod=False, dtype=<class 'numpy.float64'>)[source]

FEM solver for Laplace Eigenvalue and Poisson Equation.

Inputs can be geometry classes, which have vertices and elements. Currently, TriaMesh and TetMesh are implemented. FEM matrices (stiffness (or A) and mass matrix (or B)) are computed during the construction. After that the Eigenvalue solver (lapy.Solver.eigs) or Poisson Solver (lapy.Solver.poisson) can be called.

Parameters:
geometryTriaMesh | TetMesh

Mesh geometry.

lumpbool

If True, lump the mass matrix (diagonal).

anisofloat | tuple of shape (2,)

Anisotropy for curvature-based anisotopic Laplace. If a tuple (a_min, a_max), differentially affects the minimum and maximum curvature directions. e.g. ``(0, 50) will set scaling to 1 into the minimum curvature direction, even if the maximum curvature is large in those regions (i.e. isotropic in regions with large maximum curvature and minimum curvature close to 0, i.e. a concave cylinder).

aniso_smoothint | None

Number of smoothing iterations for curvature computation on vertices.

use_cholmodbool, default: False

If True, attempts to use the Cholesky decomposition for improved execution speed. Requires the scikit-sparse library. If it can not be found, an error will be thrown. If False, will use slower LU decomposition.

dtypenumpy.dtype, default=np.float64

Data type for the stiffness and mass matrices. Defaults to np.float64 for numerical stability. Use np.float32 explicitly if memory is constrained and precision is less critical.

Parameters:

Methods

eigs([k, sigma, which, v0, mode, rng])

Compute the linear finite-element method Laplace-Beltrami spectrum.

fem_tria_mass(tria[, lump, dtype])

Compute the sparse symmetric mass matrix of the Laplace-Beltrami operator for a given triangle mesh.

poisson([h, dtup, ntup])

Solver for the Poisson equation with boundary conditions.

Notes

The class has a static member to create the mass matrix of TriaMesh for external functions that do not need stiffness.

static fem_tria_mass(tria, lump=False, dtype=<class 'numpy.float64'>)[source]

Compute the sparse symmetric mass matrix of the Laplace-Beltrami operator for a given triangle mesh.

The sparse symmetric matrix is computed for a given triangle mesh using the linear finite element method (assuming a closed mesh or Neumann boundary condition). This function is faster than the constructor above and can be used when only a mass matrix is needed.

Parameters:
triaTriaMesh

Triangle mesh.

lumpbool, default=False

If True, B should be lumped (diagonal).

dtypenumpy.dtype, default=np.float64

Data type for the output matrix. Defaults to np.float64 for numerical stability.

Returns:
scipy.sparse.csc_matrix

Sparse symmetric positive definite mass matrix B, shape (n, n).

Parameters:
Return type:

csc_matrix

Notes

This only returns the mass matrix B of the Eigenvalue problem: A x = lambda B x. The area of the surface mesh can be obtained via B.sum().

eigs(k=10, sigma=-0.01, which='LM', v0=None, mode='normal', rng=None)[source]

Compute the linear finite-element method Laplace-Beltrami spectrum.

Parameters:
kint, default=10

The number of eigenvalues and eigenvectors desired. k must be smaller than N (number of vertices). It is not possible to compute all eigenvectors of a matrix.

sigmafloat, default=-0.01

Shift value for the shift-invert mode. The solver finds eigenvalues near sigma. The default small negative value works well for finding the smallest non-negative Laplacian eigenvalues. Adjust if convergence issues occur. The returned eigenvalues are always sorted in ascending order regardless of sigma.

whichstr, default=”LM”

Which k eigenvalues to find. With shift-invert (sigma set), "LM" selects the k eigenvalues of the original problem closest to sigma — this is the standard choice for computing the smallest Laplacian eigenvalues. Other options accepted by scipy.sparse.linalg.eigsh are "SM", "LA", "SA", and "BE".

v0np.ndarray of shape (n_vertices,), default=None

Starting vector for the ARPACK iteration. Providing a fixed vector makes results reproducible. If None, ARPACK uses a random starting vector seeded by rng.

modestr, default=”normal”

Shift-invert mode passed to scipy.sparse.linalg.eigsh. "normal" is the standard choice for Laplacian spectra. Other valid values are "buckling" and "cayley".

rngint or numpy.random.Generator, default=None

Seed or generator used to produce a reproducible starting vector when v0 is None. An integer seed is the most convenient choice, e.g. rng=0. When scipy.sparse.linalg.eigsh exposes an rng parameter, both v0 and rng are forwarded directly to it, which owns all consistency checking. Otherwise the behaviour is replicated locally: numpy.random.default_rng(rng) is used to generate a starting vector with uniform(-1, 1) (matching scipy’s own initialisation); rng is ignored when v0 is provided.

Returns:
eigenvaluesnp.ndarray

Array of k eigenvalues, shape (k,), sorted in ascending order. For closed meshes or Neumann boundary condition with the default small negative sigma, 0 will be the first eigenvalue (with constant eigenvector).

eigenvectorsnp.ndarray

Array representing the k eigenvectors, shape (n_vertices, k). The column eigenvectors[:, i] is the eigenvector corresponding to eigenvalues[i].

Parameters:
Return type:

tuple[ndarray, ndarray]

poisson(h=0.0, dtup=(), ntup=())[source]

Solver for the Poisson equation with boundary conditions.

This solver is based on the A and B Laplace matrices where A x = B h and A is a sparse symmetric positive semi-definite matrix of shape (n, n) and B is a sparse symmetric positive definite matrix of shape (n, n).

Parameters:
hfloat or np.ndarray, default=0.0

Right hand side, can be a scalar, a 1-D array of shape (n_vertices,), or a 2-D array of shape (n_vertices, n_functions) to solve for multiple functions simultaneously. The default 0.0 corresponds to the Laplace equation A x = 0.

dtuptuple, default=()

Dirichlet boundary condition as a tuple containing the index and data arrays of same length. The default, an empty tuple, corresponds to no Dirichlet condition.

ntuptuple, default=()

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.

Returns:
np.ndarray

Array with vertex values of the solution, shape (n_vertices,) for a scalar or 1-D input, or (n_vertices, n_functions) for a 2-D input.

Raises:
ValueError

If input matrices are not square or have different dimensions. If h is not scalar or array matching matrix dimensions. If dtup or ntup don’t contain exactly 2 arrays. If dtup indices are not unique. If dtup or ntup arrays have mismatched lengths. If matrix format is not CSC or CSR.

Parameters:
Return type:

ndarray

Notes

A and B are obtained via computeAB for either triangle or tetrahedral mesh.