Coverage for geometric_kernels/spaces/mesh.py: 97%
69 statements
« prev ^ index » next coverage.py v7.11.3, created at 2025-11-16 21:43 +0000
1"""
2This module provides the :class:`Mesh` space.
3"""
5import lab as B
6import numpy as np
7import potpourri3d as pp3d
8import robust_laplacian
9import scipy.sparse.linalg as sla
10from beartype.typing import Dict, Tuple
11from scipy.linalg import eigh
13from geometric_kernels.lab_extras import dtype_integer
14from geometric_kernels.spaces.base import DiscreteSpectrumSpace
15from geometric_kernels.spaces.eigenfunctions import (
16 Eigenfunctions,
17 EigenfunctionsFromEigenvectors,
18)
21class Mesh(DiscreteSpectrumSpace):
22 """
23 The GeometricKernels space representing the node set of any
24 user-provided mesh.
26 The elements of this space are represented by node indices, integer values
27 from 0 to Nv-1, where Nv is the number of nodes in the user-provided mesh.
29 Each individual eigenfunction constitutes a *level*.
31 .. note::
32 We only support the commonly used 2-dimensional meshes (discrete
33 counterparts of surfaces, 2-dimensional manifolds in a 3-dimensional
34 ambient space).
36 .. note::
37 A tutorial on how to use this space is available in the
38 :doc:`Mesh.ipynb </examples/Mesh>` notebook.
40 .. note::
41 We use `potpourri3d <https://github.com/nmwsharp/potpourri3d>`_ to
42 load meshes and mimic the interface of
43 `PyMesh <https://github.com/PyMesh/PyMesh>`_.
45 :param vertices:
46 A [Nv, 3] array of vertex coordinates, Nv is the number of vertices.
47 :param faces:
48 A [Nf, 3] array of vertex indices that represents a
49 generalized array of faces, where Nf is the number of faces.
51 .. Note:
52 Only 3 vertex indices per face are supported, i.e. mesh must be
53 triangulated.
55 .. admonition:: Citation
57 If you use this GeometricKernels space in your research, please consider
58 citing :cite:t:`borovitskiy2020`.
59 """
61 def __init__(self, vertices: np.ndarray, faces: np.ndarray):
62 self._vertices = vertices
63 if B.shape(self._vertices)[1] != 3:
64 # make sure we are in R^3.
65 raise ValueError("The last dimension (axis) of `_vertices` must be 3.")
67 self._faces = faces
68 self._eigenvalues = None
69 self._eigenfunctions = None
70 self.cache: Dict[int, Tuple[np.ndarray, np.ndarray]] = {}
72 def __str__(self):
73 return f"Mesh({self.num_vertices})"
75 def get_eigensystem(self, num: int) -> Tuple[np.ndarray, np.ndarray]:
76 """
77 Returns the first `num` eigenvalues and eigenvectors of the `robust
78 Laplacian <https://github.com/nmwsharp/nonmanifold-laplacian>`_.
79 Caches the solution to prevent re-computing the same values.
81 .. note::
82 If the `adjacency_matrix` was a sparse SciPy array, requesting
83 **all** eigenpairs will lead to a conversion of the sparse matrix
84 to a dense one due to scipy.sparse.linalg.eigsh limitations.
86 .. warning::
87 Always uses SciPy (thus CPU) for internal computations. We will
88 need to fix this in the future.
90 .. todo::
91 See warning above.
93 :param num:
94 Number of eigenpairs to return. Performs the computation at the
95 first call. Afterwards, fetches the result from cache.
97 :return:
98 A tuple of eigenvectors [nv, num], eigenvalues [num, 1].
99 """
100 if num not in self.cache:
101 L, M = robust_laplacian.mesh_laplacian(self.vertices, self.faces)
102 if L.shape[0] == num:
103 evals, evecs = eigh(L.toarray(), M.toarray())
104 else:
105 evals, evecs = sla.eigsh(L, num, M, sigma=1e-8)
106 evecs, _ = np.linalg.qr(evecs)
107 evecs *= np.sqrt(self.num_vertices)
108 evals = np.clip(
109 evals, a_min=0.0, a_max=None
110 ) # prevent small negative values
111 self.cache[num] = (evecs, evals.reshape(-1, 1))
113 return self.cache[num]
115 def get_eigenvectors(self, num: int) -> B.Numeric:
116 """
117 :param num:
118 Number of eigenvectors to return.
120 :return:
121 Array of eigenvectors, with shape [Nv, num].
122 """
123 return self.get_eigensystem(num)[0]
125 def get_eigenvalues(self, num: int) -> B.Numeric:
126 """
127 :param num:
128 Number of eigenvalues to return.
130 :return:
131 Array of eigenvalues, with shape [num, 1].
132 """
133 return self.get_eigensystem(num)[1]
135 def get_repeated_eigenvalues(self, num: int) -> B.Numeric:
136 """
137 Same as :meth:`get_eigenvalues`.
139 :param num:
140 Same as :meth:`get_eigenvalues`.
141 """
142 return self.get_eigenvalues(num)
144 def get_eigenfunctions(self, num: int) -> Eigenfunctions:
145 """
146 Returns the :class:`~.EigenfunctionsFromEigenvectors` object with
147 `num` levels (i.e., in this case, `num` eigenpairs).
149 :param num:
150 Number of levels.
151 """
152 eigenfunctions = EigenfunctionsFromEigenvectors(self.get_eigenvectors(num))
153 return eigenfunctions
155 @property
156 def num_vertices(self) -> int:
157 """
158 Number of vertices in the mesh, Nv.
159 """
160 return len(self._vertices)
162 @property
163 def num_faces(self) -> int:
164 """
165 Number of faces in the mesh, Nf.
166 """
167 return len(self._faces)
169 @property
170 def dimension(self) -> int:
171 """
172 :return:
173 2.
174 """
175 return 2
177 @property
178 def vertices(self) -> np.ndarray:
179 """
180 A [Nv, 3] array of vertex coordinates, Nv is the number of vertices.
181 """
182 return self._vertices
184 @property
185 def faces(self) -> np.ndarray:
186 """
187 A [Nf, 3] array of vertex indices that represents a generalized array of
188 faces, where Nf is the number of faces.
189 """
190 return self._faces
192 @classmethod
193 def load_mesh(cls, filename: str) -> "Mesh":
194 """
195 Construct :class:`Mesh` by loading a mesh from the file at `filename`.
197 :param filename:
198 Path to read the file from. Supported formats: `obj`,
199 `ply`, `off`, and `stl`. Format inferred automatically from the
200 file extension.
202 :return:
203 And object of class :class:`Mesh` representing the loaded mesh.
204 """
205 # load vertices and faces using potpourri3d
206 vertices, faces = pp3d.read_mesh(filename)
207 # return Mesh
208 return cls(vertices, faces)
210 def random(self, key, number):
211 key, random_vertices = B.randint(
212 key, dtype_integer(key), number, 1, lower=0, upper=self.num_vertices
213 )
214 return key, random_vertices
216 @property
217 def element_shape(self):
218 """
219 :return:
220 [1].
221 """
222 return [1]
224 @property
225 def element_dtype(self):
226 """
227 :return:
228 B.Int.
229 """
230 return B.Int