diff --git a/docs/docs/media/polygon_heat_solvers.png b/docs/docs/media/polygon_heat_solvers.png new file mode 100644 index 00000000..0b8eee94 Binary files /dev/null and b/docs/docs/media/polygon_heat_solvers.png differ diff --git a/docs/docs/pointcloud/algorithms/heat_solver.md b/docs/docs/pointcloud/algorithms/heat_solver.md index f8f90544..4aae61e3 100644 --- a/docs/docs/pointcloud/algorithms/heat_solver.md +++ b/docs/docs/pointcloud/algorithms/heat_solver.md @@ -101,7 +101,7 @@ The `SignedHeatOptions` struct allows several options: | Field | Default value |Meaning| |---|---|---| | `#!cpp bool preserveSourceNormals`| `false` | If `true`, preserve the initial curve normals at the source curve during vector diffusion. | -| `#!cpp LevelSetConstraint levelSetConstraint`| `LevelSetConstraint::ZeroSet` | Specifies how/if level sets should be preserved. Can be set to `LevelSetConstraint::ZeroSet`, `LevelSetConstraint::Multiple`, or `LevelSetConstraint::None`, corresponding to preserving the zero set, mulitple level sets (one for each curve component), or no level sets, respectively. | +| `#!cpp LevelSetConstraint levelSetConstraint`| `LevelSetConstraint::ZeroSet` | Specifies how/if level sets should be preserved. Can be set to `LevelSetConstraint::ZeroSet`, `LevelSetConstraint::Multiple`, or `LevelSetConstraint::None`, corresponding to preservation of `curves` as the zero set, as multiple level sets (one for each curve component), or no constraint, respectively. | | `#!cpp double softLevelSetWeight`| `-1` | If greater than 0, gives the weight with which the given level set constraint is "softly" enforced. | ## Scalar extension diff --git a/docs/docs/surface/algorithms/geodesic_distance.md b/docs/docs/surface/algorithms/geodesic_distance.md index e6e456f9..d42a41fc 100644 --- a/docs/docs/surface/algorithms/geodesic_distance.md +++ b/docs/docs/surface/algorithms/geodesic_distance.md @@ -182,6 +182,8 @@ This class also supports any (possibly-nonmanifold) triangle mesh as input, and `#include "geometrycentral/surface/heat_method_distance.h"` +For the polygon mesh version, see the [polygon mesh heat solver](/surface/algorithms/polygon_heat_solver); for point clouds, see the [point cloud heat solver](/pointcloud/algorithms/heat_solver/). + ### Single Solves A one-off utility function is provided which computes the distance from a source vertex using the heat method. Repeated solves or more general source data should use the stateful version below. diff --git a/docs/docs/surface/algorithms/polygon_heat_solver.md b/docs/docs/surface/algorithms/polygon_heat_solver.md new file mode 100644 index 00000000..e8d48ead --- /dev/null +++ b/docs/docs/surface/algorithms/polygon_heat_solver.md @@ -0,0 +1,115 @@ +# Heat distance and transport on general polygon meshes + +Compute signed and unsigned geodesic distance, and transport tangent vectors using fast solvers based on short-time heat flow. + +![polygon mesh heat solve results](/media/polygon_heat_solvers.png) + +These routines implement polygon mesh versions of the algorithms from: + +- [The Heat Method for Distance Computation](http://www.cs.cmu.edu/~kmcrane/Projects/HeatMethod/index.html) (distance) +- [The Vector Heat Method](https://nmwsharp.com/research/vector-heat-method) (parallel transport) +- [A Heat Method for Generalized Signed Distance](https://nzfeng.github.io/research/SignedHeatMethod/index.html) (signed distance) + +All computation is encapsulated by the `PolygonMeshHeatSolver` class, which maintains prefactored linear systems for the various methods. Some setup work is performed both on construction, and after the first query. Subsequent queries, even with different source points, will be fast. + +`#include "geometrycentral/surface/polygon_mesh_heat_solver.h"` + +For the original algorithms on triangle meshes, see the documentation for the [Unsigned Heat Method](/surface/algorithms/geodesic_distance/#heat-method-for-distance), [Vector Heat Method](/surface/algorithms/vector_heat_method), and [Signed Heat Method](/surface/algorithms/signed_heat_method). + +**Example:** Basic usage + +```cpp +#include "geometrycentral/surface/polygon_mesh_heat_solver.h" + +using namespace geometrycentral; +using namespace geometrycentral::surface; + +// Read in a polygon mesh +std::unique_ptr mesh; +std::unique_ptr geom; +std::tie(mesh, geom) = readSurfaceMesh("my_mesh.obj"); + +// Create the solver +PolygonMeshHeatSolver solver(*geom); + +// Pick a source point or two +Vertex vSource = mesh->vertex(7); +Vertex vSource2 = mesh->vertex(8); + +// Pick some source curves +std::vector> curves; +curves.push_back({mesh->vertex(11), mesh->vertex(12), mesh->vertex(15), mesh->vertex(14), mesh->vertex(13)}); +curves.push_back({mesh->vertex(17), mesh->vertex(18), mesh->vertex(19)}); + +// Compute geodesic distance +VertexData distance = solver.computeDistance(vSource); + +// Compute signed distance to a set of curves. +VertexData signedDistance = solver.computeSignedDistance(curves); + +// Compute scalar extension +VertexData extended = solver.extendScalars({{vSource, 3.}, + {vSource2, -5.}}); + +// Compute parallel transport +Vector2 sourceVec{1, 2}; +VertexData transport = solver.transportTangentVector(vSource, sourceVec); +``` + +### Constructor + +??? func "`#!cpp PolygonMeshHeatSolver::PolygonMeshHeatSolver(EmbeddedGeometryInterface& geom, double tCoef = 1.0)`" + + Create a new solver for the heat methods. Precomputation is performed at startup and lazily as needed. + + - `geom` is the geometry (and hence mesh) on which to compute. Any embedded geometry object (`VertexPositionGeometry`, etc) can be passed here. + + - `tCoef` is the time to use for short time heat flow, as a factor `m * h^2`, where `h` is the maximum between-point spacing. The default value of `1.0` is almost always sufficient. + + Algorithm options (like `tCoef`) cannot be changed after construction; create a new solver object with the new settings. + + +## Geodesic distance + +_Geodesic distance_ is the distance from a given source along the surface represented by the point cloud. Specifying multiple source points yields the distance to the nearest source. + +??? func "`#!cpp VertexData PolygonMeshHeatSolver::computeDistance(const Vertex& sourceVert)`" + + Compute the geodesic distance from `sourceVert` to all other points. + +??? func "`#!cpp VertexData PolygonMeshHeatSolver::computeDistance(const std::vector& sourceVerts)`" + + Like above, but for multiple source points. + +## Signed geodesic distance + +??? func "`#!cpp VertexData PolygonMeshHeatSolver::computeSignedDistance(const std::vector>& curves, const LevelSetConstraint& levelSetConstraint = LevelSetConstraint::ZeroSet)`" + + Compute the signed geodesic distance from a set of curves `curves` to all other points. Each curve in `curves` is a single connected component specified as an ordered sequence of points; curve orientations are derived from this order. The argument `levelSetConstraint` can be set to `LevelSetConstraint::ZeroSet`, `LevelSetConstraint::Multiple`, or `LevelSetConstraint::None`, corresponding to preservation of `curves` as the zero set, as multiple level sets (one for each curve component), or no constraint, respectively. + +## Scalar extension + +Given scalar values defined at isolated vertices in the domain, extend it to a scalar field on all vertices. Each point will take the value from the nearest source point, in the sense of geodesic distance. Note that the fast diffusion algorithm means the result is a slightly smoothed-out field. + +??? func "`#!cpp VertexData PolygonMeshHeatSolver::extendScalars(const std::vector>& sources)`" + + Given a collection of source vertices and scalars at those vertices, extends the scalar field to the whole mesh as a nearest-geodesic-neighbor interpolant. + + +## Vector extension + +Given tangent vectors defined at one or more isolated source locations on a surface, extend transport the vectors across the entire domain according to parallel transport. Each point on the domain will take the value of the nearest source point. Note that the fast diffusion algorithm means the result is a slightly smoothed-out field. + +??? func "`#!cpp VertexData PolygonMeshHeatSolver::transportTangentVector(const Vertex& sourceVert, const Vector2& sourceVector)`" + + Shorthand for the general version below when there is just one vector. + + Computes parallel transport of the given vector along shortest geodesics to the rest of the domain. + + Polar directions are defined in each vertex's tangent space. [See `EmbeddedGeometryInterface::vertexTangentBasis`](/surface/geometry#vertex-tangent-basis). + +??? func "`#!cpp VertexData PolygonMeshHeatSolver::transportTangentVectors(const std::vector>& sources)`" + + Given a collection of source vertices and tangent vectors at those vertices, extends the vector field to the whole domain as a nearest-geodesic-neighbor interpolant via parallel transport along shortest geodesics. + + Polar directions are defined in each vertex's tangent space. [See `EmbeddedGeometryInterface::vertexTangentBasis`](/surface/geometry#vertex-tangent-basis). diff --git a/docs/docs/surface/algorithms/signed_heat_method.md b/docs/docs/surface/algorithms/signed_heat_method.md index 6cf799a3..afb8a366 100644 --- a/docs/docs/surface/algorithms/signed_heat_method.md +++ b/docs/docs/surface/algorithms/signed_heat_method.md @@ -2,10 +2,11 @@ This section describes the _Signed Heat Method_ in geometry-central, which compu Note that these quantities all depend on the _intrinsic_ geometry of a surface (via the `IntrinsicGeometryInterface`). Therefore, these routines can be run on abstract geometric domains as well as traditional surfaces in 3D. -These algorithms are described in [A Heat Method for Generalized Signed Distance](https://nzfeng.github.io/research/SignedHeatMethod/SignedDistance.pdf). +This algorithm is described in [A Heat Method for Generalized Signed Distance](https://nzfeng.github.io/research/SignedHeatMethod/SignedDistance.pdf). `#include "geometrycentral/surface/signed_heat_method.h"` +For the polygon mesh version, see the [polygon mesh heat solver](/surface/algorithms/polygon_heat_solver); for point clouds, see the [point cloud heat solver](/pointcloud/algorithms/heat_solver/). ## Signed Heat Solver @@ -78,7 +79,7 @@ Options are passed in to `computeDistance` via a `SignedHeatOptions` struct, whi | Field | Default value |Meaning| |---|---|---| | `#!cpp bool preserveSourceNormals`| `false` | If `true`, preserve the initial curve normals at the source curve during vector diffusion. | -| `#!cpp LevelSetConstraint levelSetConstraint`| `LevelSetConstraint::ZeroSet` | Specifies how/if level sets should be preserved. Can be set to `LevelSetConstraint::ZeroSet`, `LevelSetConstraint::Multiple`, or `LevelSetConstraint::None`, corresponding to preserving the zero set, mulitple level sets (one for each curve component), or no level sets, respectively. | +| `#!cpp LevelSetConstraint levelSetConstraint`| `LevelSetConstraint::ZeroSet` | Specifies how/if level sets should be preserved. Can be set to `LevelSetConstraint::ZeroSet`, `LevelSetConstraint::Multiple`, or `LevelSetConstraint::None`, corresponding to preservation of the input curves as the zero set, as multiple level sets (one for each curve component), or no constraint, respectively. | | `#!cpp double softLevelSetWeight`| `-1` | If greater than 0, gives the weight with which the given level set constraint is "softly" enforced. | ## Citation diff --git a/docs/docs/surface/algorithms/vector_heat_method.md b/docs/docs/surface/algorithms/vector_heat_method.md index 1779177f..aebfca55 100644 --- a/docs/docs/surface/algorithms/vector_heat_method.md +++ b/docs/docs/surface/algorithms/vector_heat_method.md @@ -6,6 +6,7 @@ These algorithms are described in [The Vector Heat Method](http://www.cs.cmu.edu `#include "geometrycentral/surface/vector_heat_method.h"` +For the polygon mesh version, see the [polygon mesh heat solver](/surface/algorithms/polygon_heat_solver); for point clouds, see the [point cloud heat solver](/pointcloud/algorithms/heat_solver/). ## Vector Heat Solver diff --git a/docs/docs/surface/geometry/quantities.md b/docs/docs/surface/geometry/quantities.md index a1caf3d0..01ddbc75 100644 --- a/docs/docs/surface/geometry/quantities.md +++ b/docs/docs/surface/geometry/quantities.md @@ -426,7 +426,7 @@ All operators are indexed over mesh elements according to the natural iteration A $|V| \times |V|$ real matrix. Always symmetric and positive semi-definite. If and only the underlying geometry is _Delaunay_, the matrix will furthermore have all negative off-diagonal entries, satisfy a maximum principle, and be an _M-matrix_. - This is the _weak_ Laplace operator, if we use it to evalutae $\mathsf{y} \leftarrow \mathsf{L} \mathsf{x}$, $\mathsf{x}$ should hold _pointwise_ quantities at vertices, and the result $\mathsf{y}$ will contain _integrated_ values of the result in the neighborhood of each vertex. If used to solve a Poisson problem, a mass matrix (such as the lumped or Galerkin mass matrices below) are likely necessary on the right hand side. + This is the _weak_ Laplace operator, if we use it to evaluate $\mathsf{y} \leftarrow \mathsf{L} \mathsf{x}$, $\mathsf{x}$ should hold _pointwise_ quantities at vertices, and the result $\mathsf{y}$ will contain _integrated_ values of the result in the neighborhood of each vertex. If used to solve a Poisson problem, a mass matrix (such as the lumped or Galerkin mass matrices below) are likely necessary on the right hand side. Only valid on triangular meshes. @@ -543,6 +543,147 @@ In graphics and geometry processing, Crouzeix-Raviart elements have been used, f - **member:** `Eigen::SparseMatrix IntrinsicGeometryInterface::crouzeixRaviartConnectionLaplacian` - **require:** `void IntrinsicGeometryInterface::requireCrouzeixRaviartConnectionLaplacian()` +## Polygon mesh operators + +The following quantities are designed for general polygon meshes, and are defined for any `EmbeddedGeometryInterface`. On triangle meshes, they will reduce to the classical discrete exterior calculus & finite element operators. + +Two classes of polygon operators are provided: those based on [Bunge et al.'s _Polygon Laplacian Made Simple_](https://www.cs.jhu.edu/~misha/MyPapers/EUROG20.pdf), whose discretization is based on virtual refinement of the polygon mesh; and those based on [de Goes et al.'s _Discrete Differential Operators on Polygonal Meshes_](https://graphics.pixar.com/library/PolyDDG/paper.pdf), whose discretization is based on an adaptation of the virtual element method. The former uses [Astrid Bunge and Mario Botsch's implementation of their paper](https://github.com/mbotsch/polygon-laplacian), while the latter uses [David Coeurjolly, Jacques-Olivier Lachaud, and Baptiste Genest's DGtal implementation of de Goes et al.'s paper](https://www.dgtal.org/doc/stable/modulePolygonalCalculus.html). Both methods build local operators whose matrices are assembled per-polygon, so they will run out-of-the-box on non-manifold meshes (but no guarantees are provided!) + +All operators are indexed over mesh elements according to the natural iteration order of the elements, or equivalently the indices from `SurfaceMesh::getVertexIndices()` (etc). + +Here are polygon mesh operators from [Bunge et al.'s _Polygon Laplacian Made Simple_](https://www.cs.jhu.edu/~misha/MyPapers/EUROG20.pdf). + +??? func "polygon mesh Laplacian (simple)" + + ##### polygon mesh Laplacian (simple) + + The discrete Laplace operator acting on polygon meshes, using Bunge et al.'s virtual refinement method in _Polygon Laplacian Made Simple_. + + A $|V| \times |V|$ real matrix. Always symmetric and positive semi-definite. On triangle meshes, this polygon Laplacian becomes the standard cotan Laplacian. + + This is the _weak_ Laplace operator. + + Only valid on an `EmbeddedGeometryInterface`. + + - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::simplePolygonLaplacian` + - **require:** `void EmbeddedGeometryInterface::requireSimplePolygonLaplacian()` + +??? func "polygon mesh vertex lumped mass matrix (simple)" + + ##### polygon mesh vertex lumped mass matrix (simple) + + A $|V| \times |V|$ real diagonal matrix, using Bunge et al.'s virtual refinement method in _Polygon Laplacian Made Simple_. Obtained by setting each diagonal entry to the row sum in the Galerkin mass matrix. Bunge et al. note that the lumped mass matrix gives better results than the unlumped Galerkin mass matrix for most applications. + + Only valid on an `EmbeddedGeometryInterface`. + + - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::simplePolygonVertexLumpedMassMatrix` + - **require:** `void EmbeddedGeometryInterface::requireSimplePolygonVertexLumpedMassMatrix()` + +??? func "polygon mesh vertex Galerkin mass matrix (simple)" + + ##### polygon mesh vertex Galerkin mass matrix (simple) + + A $|V| \times |V|$ real matrix, using Bunge et al.'s virtual refinement method in _Polygon Laplacian Made Simple_. + + Only valid on an `EmbeddedGeometryInterface`. + + - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::simplePolygonVertexGalerkinMassMatrix` + - **require:** `void EmbeddedGeometryInterface::requireSimplePolygonVertexGalerkinMassMatrix()` + +And here are polygon mesh operators from [de Goes et al.'s _Discrete Differential Operators on Polygonal Meshes_](https://graphics.pixar.com/library/PolyDDG/paper.pdf). + +??? func "polygon mesh Laplacian" + + ##### polygon mesh Laplacian + + The discrete Laplace operator acting on polygon meshes, using de Goes et al.'s _Discrete Differential Operators on Polygonal Meshes_. + + A $|V| \times |V|$ real matrix. Always symmetric and positive semi-definite. Uses an additional parameter $\lambda$ whose default value is $1$. On triangle meshes, this polygon Laplacian becomes the standard cotan Laplacian. + + This is the _weak_ Laplace operator. + + Only valid on an `EmbeddedGeometryInterface`. + + - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonLaplacian` + - **require:** `void EmbeddedGeometryInterface::requirePolygonLaplacian()` + +??? func "polygon mesh gradient matrix" + + ##### polygon mesh gradient matrix + + The discrete gradient operator acting on polygon meshes, using de Goes et al.'s _Discrete Differential Operators on Polygonal Meshes_. + + A $3|F| \times |V|$ real matrix $\mathsf{G}$. If $\mathsf{x}\in\mathbb{R}^{|V|}$ holds pointwise quantities at vertices, then $\mathsf{y} \leftarrow \mathsf{G}\mathsf{x}\in\mathbb{R}^{3|F|}$ gives the face-wise constant gradient of $\mathsf{x}$, where the gradient in face $f$ is the vector $[\mathsf{y}_{3f}, \mathsf{y}_{3f+1}, \mathsf{y}_{3f+2}]^\top$. + + Only valid on an `EmbeddedGeometryInterface`. + + - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonGradientMatrix` + - **require:** `void EmbeddedGeometryInterface::requirePolygonGradientMatrix()` + +??? func "polygon mesh divergence matrix" + + ##### polygon mesh divergence matrix + + The discrete divergence operator acting on polygon meshes, using de Goes et al.'s _Discrete Differential Operators on Polygonal Meshes_. + + A $|V| \times 3|F|$ real matrix $\mathsf{D}$. If $\mathsf{x}\in\mathbb{R}^{3|F|}$ defines a face-wise constant vector field whose value in face $f$ is the vector $[\mathsf{x}_{3f}, \mathsf{x}_{3f+1}, \mathsf{x}_{3f+2}]^\top$, then $\mathsf{y} \leftarrow \mathsf{D}\mathsf{x}\in\mathbb{R}^{|V|}$ gives the divergence of $\mathsf{x}$ as _integrated_ quantities at vertices. To obtain pointwise quantities, one would compute $\mathsf{A}^{-1}\mathsf{D}$, where $\mathsf{A}\in\mathbb{R}^{|V|\times|V|}$ is a diagonal mass matrix of local areas at vertices. + + The divergence matrix $\mathsf{D}$ is related to the gradient matrix $\mathsf{G}$ as $\mathsf{D} = \mathsf{G}^\top\mathsf{M}$, where $\mathsf{M}\in\mathbb{R}^{3|F|}$ is a diagonal mass matrix containing face areas. Note that this assumes the convention that inflow corresponds to positive divergence, corresponding to the convention that $\mathsf{D}\mathsf{G}$ yields a positive-semidefinite Laplacian. + + Only valid on an `EmbeddedGeometryInterface`. + + - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonGradientMatrix` + - **require:** `void EmbeddedGeometryInterface::requirePolygonGradientMatrix()` + +??? func "polygon mesh vertex lumped mass matrix" + + ##### polygon mesh vertex lumped mass matrix + + A $|V| \times |V|$ real diagonal matrix, using de Goes et al.'s _Discrete Differential Operators on Polygonal Meshes_. + + Only valid on an `EmbeddedGeometryInterface`. + + - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonVertexLumpedMassMatrix` + - **require:** `void EmbeddedGeometryInterface::requirePolygonVertexLumpedMassMatrix()` + +??? func "polygon mesh vertex connection Laplacian" + + ##### polygon mesh vertex connection Laplacian + + A discrete connection Laplacian operator, which applies to vector fields defined in vertex tangent spaces; defined in de Goes et al.'s _Discrete Differential Operators on Polygonal Meshes_. Always symmetric and positive-definite. + + A $|V| \times |V|$ complex matrix. + + Given a complex vector $\mathsf{x}$ of tangent vectors at vertices, apply the operator by multiplying $\mathsf{L} * \mathsf{x}$. + + Vertex tangent spaces are defined in a similar manner to [the convention taken on triangle meshes](#vertex-tangent-spaces), namely the local $x$-axis is taken to be in the direction of `vertex.halfedge()`. On polygon meshes, the local $y$-axis is defined to be 90 degrees counterclockwise relative to the $x$-axis, rotated about the vertex normal. + + Only valid on an `EmbeddedGeometryInterface`. + + - **member:** `Eigen::SparseMatrix> EmbeddedGeometryInterface::polygonVertexConnectionLaplacian` + - **require:** `void EmbeddedGeometryInterface::requirePolygonVertexConnectionLaplacian()` + +??? func "polygon mesh DEC operators" + + ##### polygon mesh DEC operators + + These operators are the basic building blocks for _discrete exterior calculus_ on polygon meshes, using de Goes et al.'s _Discrete Differential Operators on Polygonal Meshes_. Takes in an additional parameter $\lambda$ defining a stabilization term to ensure inner products of discrete 1-forms remain positive-definite on non-triangular faces. + + **Note:** These quantities slightly deviate from the usual naming scheme for quantities. Rather than `requireD0()`, `requireD1()`, etc, there is a single `requirePolygonDECOperators()` function which manages all 7 of the members listed below. There is no `polygonHodge1Inverse`, since although `polygonHodge1` will be diagonal on triangle meshes, on general polygon meshes it will not be diagonal. Also note that the coboundary operators `polygonD0` and `polygonD1` have row/column dimension equal to the number of halfedges $|H|$ rather than $|E|$. + + The following members are constructed: + + - `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonHodge0` A $|V| \times |V|$ diagonal matrix + - `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonHodge0Inverse` A $|V| \times |V|$ diagonal matrix + - `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonHodge1` An $|H| \times |H|$ matrix, not necessarily diagonal + - `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonHodge2` An $|F| \times |F|$ diagonal matrix + - `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonHodge2Inverse` An $|F| \times |F|$ diagonal matrix + - `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonD0` An $|H| \times |V|$ matrix with $\{-1, 0, 1\}$ entries + - `Eigen::SparseMatrix EmbeddedGeometryInterface::polygonD1` An $|F| \times |H|$ matrix with $\{-1, 0, 1\}$ entries + + Only valid on an `EmbeddedGeometryInterface`. + + - **require:** `void EmbeddedGeometryInterface::requirePolygonDECOperators()` ## Extrinsic angles diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml index ebae77ae..29f17e70 100644 --- a/docs/mkdocs.yml +++ b/docs/mkdocs.yml @@ -80,6 +80,7 @@ nav: - 'Geodesic Centroidal Voronoi Tessellations' : 'surface/algorithms/geodesic_voronoi_tessellations.md' - 'Intersection' : 'surface/algorithms/intersection.md' - 'Parameterization' : 'surface/algorithms/parameterization.md' + - 'Polygon Mesh Heat Distance & Transport' : 'surface/algorithms/polygon_heat_solver.md' - 'Remeshing' : 'surface/algorithms/remeshing.md' - 'Robust Geometry' : 'surface/algorithms/robust_geometry.md' - 'Signed Heat Method' : 'surface/algorithms/signed_heat_method.md' diff --git a/include/geometrycentral/pointcloud/point_cloud_heat_solver.h b/include/geometrycentral/pointcloud/point_cloud_heat_solver.h index e1d91308..9dcd29a7 100644 --- a/include/geometrycentral/pointcloud/point_cloud_heat_solver.h +++ b/include/geometrycentral/pointcloud/point_cloud_heat_solver.h @@ -14,6 +14,7 @@ namespace geometrycentral { #ifndef SHM_H #define SHM_H +// Note: Duplicated in surface/signed_heat_method enum class LevelSetConstraint { None = 0, ZeroSet, Multiple }; struct SignedHeatOptions { diff --git a/include/geometrycentral/surface/embedded_geometry_interface.h b/include/geometrycentral/surface/embedded_geometry_interface.h index 956c8001..711d970c 100644 --- a/include/geometrycentral/surface/embedded_geometry_interface.h +++ b/include/geometrycentral/surface/embedded_geometry_interface.h @@ -54,6 +54,57 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface { void requireVertexDualMeanCurvatureNormals(); void unrequireVertexDualMeanCurvatureNormals(); + // == Polygon Operators + + // = Bunge et al. "Polygon Laplacian Made Simple" (2020), based on virtual refinement (virtual node method). + + // Laplacian + Eigen::SparseMatrix simplePolygonLaplacian; + void requireSimplePolygonLaplacian(); + void unrequireSimplePolygonLaplacian(); + + // Vertex Galerkin mass matrix (unlumped) + Eigen::SparseMatrix simplePolygonVertexGalerkinMassMatrix; + void requireSimplePolygonVertexGalerkinMassMatrix(); + void unrequireSimplePolygonVertexGalerkinMassMatrix(); + + // Vertex mass matrix (lumped) + Eigen::SparseMatrix simplePolygonVertexLumpedMassMatrix; + void requireSimplePolygonVertexLumpedMassMatrix(); + void unrequireSimplePolygonVertexLumpedMassMatrix(); + + // = de Goes et al. "Discrete Differential Operators on Polygonal Meshes" (2020), based on the virtual element method. + + // Laplacian + Eigen::SparseMatrix polygonLaplacian; + void requirePolygonLaplacian(); + void unrequirePolygonLaplacian(); + + // Gradient matrix + Eigen::SparseMatrix polygonGradientMatrix; + void requirePolygonGradientMatrix(); + void unrequirePolygonGradientMatrix(); + + // Divergence matrix + Eigen::SparseMatrix polygonDivergenceMatrix; + void requirePolygonDivergenceMatrix(); + void unrequirePolygonDivergenceMatrix(); + + // Vertex mass matrix (lumped) + Eigen::SparseMatrix polygonVertexLumpedMassMatrix; + void requirePolygonVertexLumpedMassMatrix(); + void unrequirePolygonVertexLumpedMassMatrix(); + + // Vertex connection Laplacian + Eigen::SparseMatrix> polygonVertexConnectionLaplacian; + void requirePolygonVertexConnectionLaplacian(); + void unrequirePolygonVertexConnectionLaplacian(); + + Eigen::SparseMatrix polygonHodge0, polygonHodge0Inverse, polygonHodge1, polygonHodge2, polygonHodge2Inverse, + polygonD0, polygonD1; + void requirePolygonDECOperators(); + void unrequirePolygonDECOperators(); + protected: // == Implmentations of quantities from base classes virtual void computeEdgeLengths() override; @@ -84,6 +135,96 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface { virtual void computeCornerAngles() override; virtual void computeHalfedgeCotanWeights() override; virtual void computeEdgeCotanWeights() override; + + // == Polygon Operators + + // = Bunge et al. "Polygon Laplacian Made Simple" (2020), based on virtual refinement (virtual node method). + // Copyright (C) 2020 Astrid Bunge, Philipp Herholz, Misha Kazhdan, Mario Botsch, MIT license + // (Modified to work in geometry-central. Original code can be found here: + // https://github.com/mbotsch/polygon-laplacian) + + // Laplacian + DependentQuantityD> simplePolygonLaplacianQ; + virtual void computeSimplePolygonLaplacian(); + + // Vertex mass matrix (unlumped) + DependentQuantityD> simplePolygonVertexGalerkinMassMatrixQ; + virtual void computeSimplePolygonVertexGalerkinMassMatrix(); + + // Vertex mass matrix (lumped) + DependentQuantityD> simplePolygonVertexLumpedMassMatrixQ; + virtual void computeSimplePolygonVertexLumpedMassMatrix(); + + // helper functions + FaceData virtualRefinementAreaWeights; + DependentQuantityD> virtualRefinementAreaWeightsQ; // affine weights for each virtual node + virtual void computeVirtualRefinementAreaWeights(); + virtual Eigen::MatrixXd simplePolygonMassMatrix(const Face& f); + virtual Eigen::MatrixXd simplePolygonStiffnessMatrix(const Face& f); + virtual SparseMatrix simplePolygonProlongationMatrix(); + virtual Eigen::MatrixXd polygonPositionMatrix(const Face& f); + virtual Eigen::VectorXd simplePolygonVirtualVertex(const Eigen::MatrixXd& poly) const; + + // = de Goes et al. "Discrete Differential Operators on Polygonal Meshes" (2020), based on the virtual element method. + // Use of this source code is governed by a LGPL-3.0 license. + // (Modified to work in geometry-central. Original code can be found here: + // https://github.com/DGtal-team/DGtal/blob/master/src/DGtal/dec/PolygonalCalculus.h) + + // Laplacian + DependentQuantityD> polygonLaplacianQ; + virtual void computePolygonLaplacian(); + + DependentQuantityD> polygonGradientMatrixQ; + virtual void computePolygonGradientMatrix(); + + DependentQuantityD> polygonDivergenceMatrixQ; + virtual void computePolygonDivergenceMatrix(); + + // Vertex mass matrix (lumped) + DependentQuantityD> polygonVertexLumpedMassMatrixQ; + virtual void computePolygonVertexLumpedMassMatrix(); + + // Vertex connection Laplacian + DependentQuantityD>> polygonVertexConnectionLaplacianQ; + virtual void computePolygonVertexConnectionLaplacian(); + + // DEC Operators + std::array*, 7> polygonDECOperatorArray; + DependentQuantityD*, 7>> polygonDECOperatorsQ; + virtual void computePolygonDECOperators(); + + // helper functions + const double polygonLambda = 1.0; + VertexData polygonVertexNormals; + DependentQuantityD> polygonVertexNormalsQ; + virtual void computePolygonVertexNormals(); // area-weighted normals + virtual Eigen::MatrixXd polygonPerFaceLaplacian(const Face& f); + virtual Eigen::MatrixXd polygonPerFaceInnerProductMatrix(const Face& f); + virtual Eigen::MatrixXd polygonProjectionMatrix(const Face& f); + virtual Eigen::MatrixXd polygonPerFaceGradientMatrix(const Face& f); + virtual Eigen::MatrixXd polygonCoGradientMatrix(const Face& f); + virtual Eigen::MatrixXd polygonAveragingMatrix(const Face& f) const; + virtual Eigen::MatrixXd polygonDerivativeMatrix(const Face& f) const; + virtual Eigen::MatrixXd polygonEdgeVectorMatrix(const Face& f); + virtual Eigen::MatrixXd polygonEdgeMidpointMatrix(const Face& f); + virtual Eigen::MatrixXd polygonFlat(const Face& f); + virtual Eigen::MatrixXd polygonSharp(const Face& f); + virtual Eigen::Vector3d polygonCentroid(const Face& f); + // connections + virtual Eigen::MatrixXd polygonPerFaceConnectionLaplacian(const Face& f); + virtual Eigen::MatrixXd polygonBlockConnection(const Face& f); + virtual Eigen::MatrixXd polygonCovariantGradient(const Face& f); + virtual Eigen::MatrixXd polygonCovariantProjection(const Face& f); + // tangent space helpers + virtual Eigen::MatrixXd Tv(const Vertex& v); + virtual Eigen::MatrixXd Tf(const Face& f); + virtual Eigen::Matrix2d Rvf(const Vertex& v, const Face& f); + virtual Eigen::Matrix3d Qvf(const Vertex& v, const Face& f); + // helpers to the helper functions: generic linear algebra stuff, though probably wouldn't find much use elsewhere + // so keeping them here -- also they use Eigen::Vectors here for matrix-multiply compatibility. + virtual Eigen::Matrix3d bracket(const Eigen::Vector3d& n) const; + virtual Eigen::Vector3d project(const Eigen::Vector3d& u, const Eigen::Vector3d& n) const; + virtual Eigen::MatrixXd kroneckerWithI2(const Eigen::MatrixXd& M) const; }; diff --git a/include/geometrycentral/surface/polygon_mesh_heat_solver.h b/include/geometrycentral/surface/polygon_mesh_heat_solver.h new file mode 100644 index 00000000..bc951ae3 --- /dev/null +++ b/include/geometrycentral/surface/polygon_mesh_heat_solver.h @@ -0,0 +1,63 @@ +#pragma once + +#include "geometrycentral/numerical/linear_solvers.h" +#include "geometrycentral/surface/embedded_geometry_interface.h" +#include "geometrycentral/surface/signed_heat_method.h" +#include "geometrycentral/surface/surface_mesh.h" + +namespace geometrycentral { +namespace surface { + +class PolygonMeshHeatSolver { + +public: + // === Constructor + PolygonMeshHeatSolver(EmbeddedGeometryInterface& geom, double tCoef = 1.0); + + // === Methods + + // Solve for distance from a single vertex (or collection of points) + VertexData computeDistance(const Vertex& sourceVert); + VertexData computeDistance(const std::vector& sourceVerts); + + // Scalar Extension + VertexData extendScalars(const std::vector>& sources); + + // Compute parallel transport along shortest geodesics from sources at points + VertexData transportTangentVector(const Vertex& sourceVert, const Vector2& sourceVector); + VertexData transportTangentVectors(const std::vector>& sources); + + // Solve for signed distance from a curve comprising a sequence of vertices. + VertexData computeSignedDistance(const std::vector>& curves, + const LevelSetConstraint& levelSetConstraint = LevelSetConstraint::ZeroSet); + + // === Options and parameters + + const double tCoef; // the time parameter used for heat flow, measured as time = tCoef * maxDiagonalLength + // default: 1.0 + +private: + // === Members + + // Input mesh and geometry + SurfaceMesh& mesh; + EmbeddedGeometryInterface& geom; + + // Parameters + double shortTime; + + // Solvers + void ensureHaveScalarHeatSolver(); + void ensureHaveVectorHeatSolver(); + void ensureHavePoissonSolver(); + std::unique_ptr>> vectorHeatSolver; + std::unique_ptr> scalarHeatSolver, poissonSolver; + SparseMatrix massMat, laplaceMat; + + // Helpers + void buildSignedCurveSource(const std::vector& curve, Vector>& X0) const; + double computeAverageValue(const std::vector>& curves, const Vector& u); +}; + +} // namespace surface +} // namespace geometrycentral \ No newline at end of file diff --git a/include/geometrycentral/surface/surface_point.ipp b/include/geometrycentral/surface/surface_point.ipp index 9d324f4d..e3cf28ce 100644 --- a/include/geometrycentral/surface/surface_point.ipp +++ b/include/geometrycentral/surface/surface_point.ipp @@ -263,11 +263,11 @@ inline bool SurfacePoint::operator==(const SurfacePoint& other) const { break; } case SurfacePointType::Edge: { - return edge == other.edge && abs(tEdge - other.tEdge) < eps; + return edge == other.edge && tEdge == other.tEdge; break; } case SurfacePointType::Face: { - return face == other.face && (faceCoords - other.faceCoords).norm() < eps; + return face == other.face && faceCoords == other.faceCoords; break; } } @@ -306,8 +306,6 @@ inline Edge sharedEdge(const SurfacePoint& pA, const SurfacePoint& pB) { if (pA.type == SurfacePointType::Face || pB.type == SurfacePointType::Face) return Edge(); - // This is kind of ugly code, esp. since there's only two options for the switch statement. But a nested-if looks even - // worse... switch (pA.type) { case SurfacePointType::Vertex: { switch (pB.type) { diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index e11380b7..320427c6 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -32,6 +32,7 @@ SET(SRCS surface/heat_method_distance.cpp surface/vector_heat_method.cpp surface/signed_heat_method.cpp + surface/polygon_mesh_heat_solver.cpp surface/geodesic_centroidal_voronoi_tessellation.cpp surface/trace_geodesic.cpp surface/normal_coordinates.cpp @@ -125,6 +126,7 @@ SET(HEADERS ${INCLUDE_ROOT}/surface/mesh_ray_tracer.h ${INCLUDE_ROOT}/surface/parameterize.h ${INCLUDE_ROOT}/surface/poisson_disk_sampler.h + ${INCLUDE_ROOT}/surface/polygon_mesh_heat_solver.h ${INCLUDE_ROOT}/surface/quadric_error_simplification.h ${INCLUDE_ROOT}/surface/remeshing.h ${INCLUDE_ROOT}/surface/rich_surface_mesh_data.h diff --git a/src/surface/embedded_geometry_interface.cpp b/src/surface/embedded_geometry_interface.cpp index 81e11d94..f9f73cfe 100644 --- a/src/surface/embedded_geometry_interface.cpp +++ b/src/surface/embedded_geometry_interface.cpp @@ -17,8 +17,22 @@ EmbeddedGeometryInterface::EmbeddedGeometryInterface(SurfaceMesh& mesh_) : vertexNormalsQ (&vertexNormals, std::bind(&EmbeddedGeometryInterface::computeVertexNormals, this), quantities), faceTangentBasisQ (&faceTangentBasis, std::bind(&EmbeddedGeometryInterface::computeFaceTangentBasis, this), quantities), vertexTangentBasisQ (&vertexTangentBasis, std::bind(&EmbeddedGeometryInterface::computeVertexTangentBasis, this), quantities), - vertexDualMeanCurvatureNormalsQ (&vertexDualMeanCurvatureNormals, std::bind(&EmbeddedGeometryInterface::computeVertexDualMeanCurvatureNormals, this), quantities) + vertexDualMeanCurvatureNormalsQ (&vertexDualMeanCurvatureNormals, std::bind(&EmbeddedGeometryInterface::computeVertexDualMeanCurvatureNormals, this), quantities), + + simplePolygonLaplacianQ (&simplePolygonLaplacian, std::bind(&EmbeddedGeometryInterface::computeSimplePolygonLaplacian, this), quantities), + simplePolygonVertexGalerkinMassMatrixQ (&simplePolygonVertexGalerkinMassMatrix, std::bind(&EmbeddedGeometryInterface::computeSimplePolygonVertexGalerkinMassMatrix, this), quantities), + simplePolygonVertexLumpedMassMatrixQ (&simplePolygonVertexLumpedMassMatrix, std::bind(&EmbeddedGeometryInterface::computeSimplePolygonVertexLumpedMassMatrix, this), quantities), + virtualRefinementAreaWeightsQ (&virtualRefinementAreaWeights, std::bind(&EmbeddedGeometryInterface::computeVirtualRefinementAreaWeights, this), quantities), + + polygonLaplacianQ (&polygonLaplacian, std::bind(&EmbeddedGeometryInterface::computePolygonLaplacian, this), quantities), + polygonGradientMatrixQ (&polygonGradientMatrix, std::bind(&EmbeddedGeometryInterface::computePolygonGradientMatrix, this), quantities), + polygonDivergenceMatrixQ (&polygonDivergenceMatrix, std::bind(&EmbeddedGeometryInterface::computePolygonDivergenceMatrix, this), quantities), + polygonVertexLumpedMassMatrixQ (&polygonVertexLumpedMassMatrix, std::bind(&EmbeddedGeometryInterface::computePolygonVertexLumpedMassMatrix, this), quantities), + polygonVertexConnectionLaplacianQ (&polygonVertexConnectionLaplacian, std::bind(&EmbeddedGeometryInterface::computePolygonVertexConnectionLaplacian, this), quantities), + polygonDECOperatorArray{&polygonHodge0, &polygonHodge0Inverse, &polygonHodge1, &polygonHodge2, &polygonHodge2Inverse, &polygonD0, &polygonD1}, + polygonDECOperatorsQ(&polygonDECOperatorArray, std::bind(&EmbeddedGeometryInterface::computePolygonDECOperators, this), quantities), + polygonVertexNormalsQ (&polygonVertexNormals, std::bind(&EmbeddedGeometryInterface::computePolygonVertexNormals, this), quantities) {} // clang-format on @@ -184,7 +198,6 @@ void EmbeddedGeometryInterface::computeVertexTangentBasis() { } halfedgeVectorsInVertexQ.ensureHave(); - for (Vertex v : mesh.vertices()) { // For general polygons, take the average of each edge vector projected to tangent plane @@ -249,18 +262,15 @@ void EmbeddedGeometryInterface::computeFaceAreas() { faceAreas = FaceData(mesh); for (Face f : mesh.faces()) { - - // WARNING: Logic duplicated between cached and immediate version - Halfedge he = f.halfedge(); - Vector3 pA = vertexPositions[he.vertex()]; - he = he.next(); - Vector3 pB = vertexPositions[he.vertex()]; - he = he.next(); - Vector3 pC = vertexPositions[he.vertex()]; - - GC_SAFETY_ASSERT(he.next() == f.halfedge(), "faces must be triangular"); - - double area = 0.5 * norm(cross(pB - pA, pC - pA)); + Vector3 N = {0, 0, 0}; + for (Halfedge he : f.adjacentHalfedges()) { + Vertex vA = he.vertex(); + Vertex vB = he.next().vertex(); + Vector3 pA = vertexPositions[vA]; + Vector3 pB = vertexPositions[vB]; + N += cross(pA, pB); + } + double area = 0.5 * norm(N); faceAreas[f] = area; } } @@ -278,11 +288,11 @@ void EmbeddedGeometryInterface::computeCornerAngles() { Vector3 pA = vertexPositions[he.vertex()]; he = he.next(); Vector3 pB = vertexPositions[he.vertex()]; - he = he.next(); + do { + he = he.next(); + } while (he.next() != c.halfedge()); Vector3 pC = vertexPositions[he.vertex()]; - GC_SAFETY_ASSERT(he.next() == c.halfedge(), "faces must be triangular"); - double q = dot(unit(pB - pA), unit(pC - pA)); q = clamp(q, -1.0, 1.0); double angle = std::acos(q); @@ -349,6 +359,758 @@ void EmbeddedGeometryInterface::computeEdgeCotanWeights() { } } +// === Polygon Operators + +// = Bunge et al. "Polygon Laplacian Made Simple" (2020), based on virtual refinement (virtual node method). +// Copyright (C) 2020 Astrid Bunge, Philipp Herholz, Misha Kazhdan, Mario Botsch, MIT license +// (Modified to work in geometry-central. Original code can be found here: https://github.com/mbotsch/polygon-laplacian) + +// Laplacian +void EmbeddedGeometryInterface::computeSimplePolygonLaplacian() { + vertexIndicesQ.ensureHave(); + + size_t V = mesh.nVertices(); + simplePolygonLaplacian = Eigen::SparseMatrix(V, V); + std::vector> triplets; + std::vector vIndices; // indices of vertices of polygon face + Eigen::MatrixXd Si; // local per-polygon matrix + for (Face f : mesh.faces()) { + vIndices.clear(); + for (Vertex v : f.adjacentVertices()) vIndices.push_back(vertexIndices[v]); + size_t n = f.degree(); + // Get local stiffness matrix. + Si = simplePolygonStiffnessMatrix(f); + // Add contribution to global mass matrix. + for (size_t j = 0; j < n; j++) { + for (size_t i = 0; i < n; i++) { + triplets.emplace_back(vIndices[i], vIndices[j], Si(i, j)); + } + } + } + simplePolygonLaplacian.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requireSimplePolygonLaplacian() { simplePolygonLaplacianQ.require(); } +void EmbeddedGeometryInterface::unrequireSimplePolygonLaplacian() { simplePolygonLaplacianQ.unrequire(); } + + +// Vertex Galerkin mass matrix (unlumped) +void EmbeddedGeometryInterface::computeSimplePolygonVertexGalerkinMassMatrix() { + vertexIndicesQ.ensureHave(); + + size_t V = mesh.nVertices(); + simplePolygonVertexGalerkinMassMatrix = Eigen::SparseMatrix(V, V); + std::vector> triplets; + std::vector vIndices; // indices of vertices of polygon face + Eigen::MatrixXd Mi; // local per-polygon matrix + for (Face f : mesh.faces()) { + vIndices.clear(); + for (Vertex v : f.adjacentVertices()) vIndices.push_back(vertexIndices[v]); + size_t n = f.degree(); + // Get local mass matrix. + Mi = simplePolygonMassMatrix(f); + // Add contribution to global mass matrix. + for (size_t j = 0; j < n; j++) { + for (size_t i = 0; i < n; i++) { + triplets.emplace_back(vIndices[i], vIndices[j], Mi(i, j)); + } + } + } + simplePolygonVertexGalerkinMassMatrix.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requireSimplePolygonVertexGalerkinMassMatrix() { + simplePolygonVertexGalerkinMassMatrixQ.require(); +} +void EmbeddedGeometryInterface::unrequireSimplePolygonVertexGalerkinMassMatrix() { + simplePolygonVertexGalerkinMassMatrixQ.unrequire(); +} + + +// Vertex mass matrix (lumped) +void EmbeddedGeometryInterface::computeSimplePolygonVertexLumpedMassMatrix() { + vertexIndicesQ.ensureHave(); + + size_t V = mesh.nVertices(); + simplePolygonVertexLumpedMassMatrix = Eigen::SparseMatrix(V, V); + std::vector> triplets; + std::vector vIndices; // indices of vertices of polygon face + Eigen::MatrixXd Mi; // local per-polygon matrix + for (Face f : mesh.faces()) { + vIndices.clear(); + for (Vertex v : f.adjacentVertices()) vIndices.push_back(vertexIndices[v]); + size_t n = f.degree(); + // Get local mass matrix. + Mi = simplePolygonMassMatrix(f); + // Add contribution to global mass matrix. + for (size_t j = 0; j < n; j++) { + for (size_t i = 0; i < n; i++) { + triplets.emplace_back(vIndices[i], vIndices[i], Mi(i, j)); + } + } + } + simplePolygonVertexLumpedMassMatrix.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requireSimplePolygonVertexLumpedMassMatrix() { + simplePolygonVertexLumpedMassMatrixQ.require(); +} +void EmbeddedGeometryInterface::unrequireSimplePolygonVertexLumpedMassMatrix() { + simplePolygonVertexLumpedMassMatrixQ.unrequire(); +} + +// Helper functions + +Eigen::MatrixXd EmbeddedGeometryInterface::simplePolygonMassMatrix(const Face& f) { + virtualRefinementAreaWeightsQ.ensureHave(); + + size_t n = f.degree(); + Eigen::MatrixXd poly = polygonPositionMatrix(f); + Eigen::MatrixXd M = Eigen::MatrixXd::Zero(n, n); + const Eigen::VectorXd& weights = virtualRefinementAreaWeights[f]; + Eigen::Vector3d virtualVertex = poly.transpose() * weights; + Eigen::VectorXd ln = Eigen::VectorXd::Zero(n + 1); + double l[3], l2[3]; // lengths, lengths squared + // Build triangle fan mass and cotan matrices + for (size_t i = 0; i < n; i++) { + const size_t i1 = (i + 1) % n; + l2[2] = (poly.row(i) - poly.row(i1)).squaredNorm(); + l2[0] = (poly.row(i1) - virtualVertex.transpose()).squaredNorm(); + l2[1] = (poly.row(i) - virtualVertex.transpose()).squaredNorm(); + l[0] = std::sqrt(l2[0]); + l[1] = std::sqrt(l2[1]); + l[2] = std::sqrt(l2[2]); + const double arg = + (l[0] + (l[1] + l[2])) * (l[2] - (l[0] - l[1])) * (l[2] + (l[0] - l[1])) * (l[0] + (l[1] - l[2])); + const double area = 0.25 * std::sqrt(arg); + l[0] = 1.0 / 6.0 * area; + l[1] = 1.0 / 12.0 * area; + M(i1, i1) += 1.0 / 6.0 * area; + M(i, i) += 1.0 / 6.0 * area; + M(i1, i) += 1.0 / 12.0 * area; + M(i, i1) += 1.0 / 12.0 * area; + ln(i1) += l[1]; + ln(i) += l[1]; + ln(n) += l[0]; + } + // Apply prolongation + for (size_t j = 0; j < n; ++j) + for (size_t i = 0; i < n; ++i) M(i, j) += weights(i) * ln(j) + weights(j) * ln(i) + weights(i) * weights(j) * ln(n); + + return M; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::simplePolygonStiffnessMatrix(const Face& f) { + virtualRefinementAreaWeightsQ.ensureHave(); + + size_t n = f.degree(); + Eigen::MatrixXd poly = polygonPositionMatrix(f); + Eigen::MatrixXd S = Eigen::MatrixXd::Zero(n, n); + const Eigen::VectorXd& weights = virtualRefinementAreaWeights[f]; + Eigen::Vector3d virtualVertex = poly.transpose() * weights; + Eigen::VectorXd ln = Eigen::VectorXd::Zero(n + 1); + double l[3], l2[3]; // lengths, lengths squared + // Build triangle fan mass and cotan matrices + for (size_t i = 0; i < n; i++) { + const size_t i1 = (i + 1) % n; + l2[2] = (poly.row(i) - poly.row(i1)).squaredNorm(); + l2[0] = (poly.row(i1) - virtualVertex.transpose()).squaredNorm(); + l2[1] = (poly.row(i) - virtualVertex.transpose()).squaredNorm(); + l[0] = std::sqrt(l2[0]); + l[1] = std::sqrt(l2[1]); + l[2] = std::sqrt(l2[2]); + const double arg = + (l[0] + (l[1] + l[2])) * (l[2] - (l[0] - l[1])) * (l[2] + (l[0] - l[1])) * (l[0] + (l[1] - l[2])); + const double area = 0.5 * std::sqrt(arg); + if (area > 1e-7) { + l[0] = 0.25 * (l2[1] + l2[2] - l2[0]) / area; + l[1] = 0.25 * (l2[2] + l2[0] - l2[1]) / area; + l[2] = 0.25 * (l2[0] + l2[1] - l2[2]) / area; + + S(i1, i1) += l[0]; + S(i, i) += l[1]; + S(i1, i) -= l[2]; + S(i, i1) -= l[2]; + S(i, i) += l[2]; + S(i1, i1) += l[2]; + + ln(i1) -= l[0]; + ln(i) -= l[1]; + ln(n) += l[0] + l[1]; + } + } + // Apply prolongation + for (size_t j = 0; j < n; ++j) + for (size_t i = 0; i < n; ++i) S(i, j) += weights(i) * ln(j) + weights(j) * ln(i) + weights(i) * weights(j) * ln(n); + + return S; +} + +SparseMatrix EmbeddedGeometryInterface::simplePolygonProlongationMatrix() { + virtualRefinementAreaWeightsQ.ensureHave(); + vertexIndicesQ.ensureHave(); + + size_t V = mesh.nVertices(); + size_t F = mesh.nFaces(); + std::vector> triplets; + SparseMatrix P(V + F, V); + for (size_t i = 0; i < V; i++) triplets.emplace_back(i, i, 1); + int j = 0; + for (Face f : mesh.faces()) { + Eigen::VectorXd weights = virtualRefinementAreaWeights[f]; + int i = 0; + for (Vertex v : f.adjacentVertices()) { + size_t vIdx = vertexIndices[v]; + triplets.emplace_back(V + j, vIdx, weights[i]); + i++; + } + j++; + } + P.setFromTriplets(triplets.begin(), triplets.end()); + return P; +} + +void EmbeddedGeometryInterface::computeVirtualRefinementAreaWeights() { + vertexPositionsQ.ensureHave(); + + virtualRefinementAreaWeights = FaceData(mesh); + + for (Face f : mesh.faces()) { + Eigen::MatrixXd poly = polygonPositionMatrix(f); + Eigen::VectorXd weights = simplePolygonVirtualVertex(poly); + virtualRefinementAreaWeights[f] = weights; + } +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonPositionMatrix(const Face& f) { + vertexPositionsQ.ensureHave(); + + Eigen::MatrixXd poly(f.degree(), 3); + int i = 0; + for (Vertex v : f.adjacentVertices()) { + for (int j = 0; j < 3; j++) { + poly(i, j) = vertexPositions[v][j]; + } + i++; + } + return poly; +} + +Eigen::VectorXd EmbeddedGeometryInterface::simplePolygonVirtualVertex(const Eigen::MatrixXd& poly) const { + + // Given a polygon face, computes the affine weights that determine the position of the virtual vertex that minimizes + // the sum of the squared areas of the triangles in the induced triangle fan. While the location of this vertex (the + // minimizer) is unique, its expression as an affine combination of the polygon verties may not be -- regularize by + // picking the weights with minimum L_2 norm, which encourages the weights to be as uniform as possible. + + int n = poly.rows(); + Eigen::VectorXd weights(n); + Eigen::MatrixXd J(n, n); + Eigen::VectorXd b(n); + for (int i = 0; i < n; i++) { + Eigen::Vector3d pk = poly.row(i); + + double Bk1_d2 = 0.0; + double Bk1_d1 = 0.0; + + double Bk2_d0 = 0.0; + double Bk2_d2 = 0.0; + + double Bk3_d0 = 0.0; + double Bk3_d1 = 0.0; + + double CBk = 0.0; + Eigen::Vector3d d = Eigen::MatrixXd::Zero(3, 1); + + for (int j = 0; j < n; j++) { + Eigen::Vector3d pi = poly.row(j); + Eigen::Vector3d pj = poly.row((j + 1) % n); + d = pi - pj; + + double Bik1 = d(1) * pk(2) - d(2) * pk(1); + double Bik2 = d(2) * pk(0) - d(0) * pk(2); + double Bik3 = d(0) * pk(1) - d(1) * pk(0); + + double Ci1 = d(1) * pi(2) - d(2) * pi(1); + double Ci2 = d(2) * pi(0) - d(0) * pi(2); + double Ci3 = d(0) * pi(1) - d(1) * pi(0); + + Bk1_d1 += d(1) * Bik1; + Bk1_d2 += d(2) * Bik1; + + Bk2_d0 += d(0) * Bik2; + Bk2_d2 += d(2) * Bik2; + + Bk3_d0 += d(0) * Bik3; + Bk3_d1 += d(1) * Bik3; + + CBk += Ci1 * Bik1 + Ci2 * Bik2 + Ci3 * Bik3; + } + for (int k = 0; k < n; k++) { + Eigen::Vector3d xj = poly.row(k); + J(i, k) = + 0.5 * (xj(2) * Bk1_d1 - xj(1) * Bk1_d2 + xj(0) * Bk2_d2 - xj(2) * Bk2_d0 + xj(1) * Bk3_d0 - xj(0) * Bk3_d1); + } + b(i) = 0.5 * CBk; + } + + Eigen::MatrixXd M(n + 1, n); + M.block(0, 0, n, n) = 4 * J; + M.block(n, 0, 1, n).setOnes(); + + Eigen::VectorXd b_(n + 1); + b_.block(0, 0, n, 1) = 4 * b; + + b_(n) = 1.; + weights = M.completeOrthogonalDecomposition().solve(b_).topRows(n); + + return weights; +} + + +// = de Goes et al. "Discrete Differential Operators on Polygonal Meshes" (2020), based on the virtual element method. +// Use of this source code is governed by a LGPL-3.0 license. +// (Modified to work in geometry-central. Original code can be found here: +// https://github.com/DGtal-team/DGtal/blob/master/src/DGtal/dec/PolygonalCalculus.h) + +// Laplacian +void EmbeddedGeometryInterface::computePolygonLaplacian() { + vertexIndicesQ.ensureHave(); + + size_t V = mesh.nVertices(); + polygonLaplacian = Eigen::SparseMatrix(V, V); + std::vector> triplets; + // Assemble per-face operators. + Eigen::MatrixXd Lf; // local per-polygon matrix + std::vector vIndices; // indices of vertices of polygon face + for (Face f : mesh.faces()) { + vIndices.clear(); + for (Vertex v : f.adjacentVertices()) vIndices.push_back(vertexIndices[v]); + size_t n = f.degree(); + // Add contribution to global matrix. + Lf = polygonPerFaceLaplacian(f); + for (size_t j = 0; j < n; j++) { + for (size_t i = 0; i < n; i++) { + triplets.emplace_back(vIndices[i], vIndices[j], Lf(i, j)); + } + } + } + polygonLaplacian.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requirePolygonLaplacian() { polygonLaplacianQ.require(); } +void EmbeddedGeometryInterface::unrequirePolygonLaplacian() { polygonLaplacianQ.unrequire(); } + + +// Gradient matrix +void EmbeddedGeometryInterface::computePolygonGradientMatrix() { + vertexIndicesQ.ensureHave(); + faceIndicesQ.ensureHave(); + + size_t V = mesh.nVertices(); + size_t F = mesh.nFaces(); + polygonGradientMatrix = Eigen::SparseMatrix(3 * F, V); + std::vector> triplets; + // Assemble per-face operators. + Eigen::MatrixXd Gf; // local per-polygon matrix + std::vector vIndices; // indices of vertices of polygon face + for (Face f : mesh.faces()) { + size_t fIdx = faceIndices[f]; + vIndices.clear(); + for (Vertex v : f.adjacentVertices()) vIndices.push_back(vertexIndices[v]); + size_t n = f.degree(); + // Add contribution to global matrix. + Gf = polygonPerFaceGradientMatrix(f); // 3 x nf + for (size_t i = 0; i < 3; i++) { + for (size_t j = 0; j < n; j++) { + triplets.emplace_back(3 * fIdx + i, vIndices[j], Gf(i, j)); + } + } + } + polygonGradientMatrix.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requirePolygonGradientMatrix() { polygonGradientMatrixQ.require(); } +void EmbeddedGeometryInterface::unrequirePolygonGradientMatrix() { polygonGradientMatrixQ.unrequire(); } + + +// Gradient matrix +void EmbeddedGeometryInterface::computePolygonDivergenceMatrix() { + vertexIndicesQ.ensureHave(); + faceIndicesQ.ensureHave(); + faceAreasQ.ensureHave(); + + size_t V = mesh.nVertices(); + size_t F = mesh.nFaces(); + polygonDivergenceMatrix = Eigen::SparseMatrix(V, 3 * F); + std::vector> triplets; + // Assemble per-face operators. + Eigen::MatrixXd Gf; // local per-polygon matrix + std::vector vIndices; // indices of vertices of polygon face + for (Face f : mesh.faces()) { + size_t fIdx = faceIndices[f]; + double area = faceAreas[f]; + vIndices.clear(); + for (Vertex v : f.adjacentVertices()) vIndices.push_back(vertexIndices[v]); + size_t n = f.degree(); + // Add contribution to global matrix. + Gf = polygonPerFaceGradientMatrix(f); // 3 x nf + for (size_t i = 0; i < 3; i++) { + for (size_t j = 0; j < n; j++) { + triplets.emplace_back(vIndices[j], 3 * fIdx + i, Gf(i, j) * area); + } + } + } + polygonDivergenceMatrix.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requirePolygonDivergenceMatrix() { polygonDivergenceMatrixQ.require(); } +void EmbeddedGeometryInterface::unrequirePolygonDivergenceMatrix() { polygonDivergenceMatrixQ.unrequire(); } + + +// Vertex mass matrix (lumped) +void EmbeddedGeometryInterface::computePolygonVertexLumpedMassMatrix() { + vertexIndicesQ.ensureHave(); + faceAreasQ.ensureHave(); + + size_t V = mesh.nVertices(); + Eigen::VectorXd hodge0 = Eigen::VectorXd::Zero(V); + for (Face f : mesh.faces()) { + double w = faceAreas[f] / f.degree(); + for (Vertex v : f.adjacentVertices()) { + size_t vIdx = vertexIndices[v]; + hodge0[vIdx] += w; + } + } + polygonVertexLumpedMassMatrix = hodge0.asDiagonal(); +} +void EmbeddedGeometryInterface::requirePolygonVertexLumpedMassMatrix() { polygonVertexLumpedMassMatrixQ.require(); } +void EmbeddedGeometryInterface::unrequirePolygonVertexLumpedMassMatrix() { polygonVertexLumpedMassMatrixQ.unrequire(); } + + +// vertex connection Laplacian +void EmbeddedGeometryInterface::computePolygonVertexConnectionLaplacian() { + vertexIndicesQ.ensureHave(); + + size_t V = mesh.nVertices(); + polygonVertexConnectionLaplacian = Eigen::SparseMatrix>(V, V); + std::vector>> triplets, tripletsTest; + Eigen::SparseMatrix> testMat(2 * V, 2 * V); + Eigen::MatrixXd Lf; + std::vector vIndices; + for (Face f : mesh.faces()) { + vIndices.clear(); + for (Vertex v : f.adjacentVertices()) vIndices.push_back(vertexIndices[v]); + size_t n = f.degree(); + Lf = polygonPerFaceConnectionLaplacian(f); + for (size_t j = 0; j < n; j++) { + for (size_t i = 0; i < n; i++) { + double re = Lf(2 * i, 2 * j); + double im = Lf(2 * i + 1, 2 * j); + // Split up into symmetric contributions to ensure Hermitian-ness. + triplets.emplace_back(vIndices[i], vIndices[j], 0.5 * std::complex(re, im)); + triplets.emplace_back(vIndices[j], vIndices[i], 0.5 * std::complex(re, -im)); + } + } + } + polygonVertexConnectionLaplacian.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requirePolygonVertexConnectionLaplacian() { + polygonVertexConnectionLaplacianQ.require(); +} +void EmbeddedGeometryInterface::unrequirePolygonVertexConnectionLaplacian() { + polygonVertexConnectionLaplacianQ.unrequire(); +} + + +void EmbeddedGeometryInterface::computePolygonDECOperators() { + vertexIndicesQ.ensureHave(); + edgeIndicesQ.ensureHave(); + faceIndicesQ.ensureHave(); + halfedgeIndicesQ.ensureHave(); + faceAreasQ.ensureHave(); + + size_t V = mesh.nVertices(); + size_t E = mesh.nEdges(); + size_t F = mesh.nFaces(); + size_t H = mesh.nHalfedges(); // technically only needs to be nInteriorHalfedges(), but interior halfedges aren't + // indexed densely + std::vector> tripletsD0, tripletsD1, tripletsH1; + + // exterior derivatives + polygonD0 = Eigen::SparseMatrix(H, V); + for (Face f : mesh.faces()) { + for (Halfedge he : f.adjacentHalfedges()) { + size_t hIdx = halfedgeIndices[he]; + size_t vA = vertexIndices[he.tailVertex()]; + size_t vB = vertexIndices[he.tipVertex()]; + tripletsD0.emplace_back(hIdx, vA, -1.); + tripletsD0.emplace_back(hIdx, vB, 1.); + } + } + polygonD0.setFromTriplets(tripletsD0.begin(), tripletsD0.end()); + + polygonD1 = Eigen::SparseMatrix(F, H); + for (Face f : mesh.faces()) { + size_t fIdx = faceIndices[f]; + for (Halfedge he : f.adjacentHalfedges()) { + size_t hIdx = halfedgeIndices[he]; + tripletsD1.emplace_back(fIdx, hIdx, 1.); + } + } + polygonD1.setFromTriplets(tripletsD1.begin(), tripletsD1.end()); + + // hodge0 + Eigen::VectorXd h0 = Eigen::VectorXd::Zero(V); + for (Face f : mesh.faces()) { + double w = faceAreas[f] / f.degree(); + for (Vertex v : f.adjacentVertices()) { + size_t vIdx = vertexIndices[v]; + h0[vIdx] += w; + } + } + polygonHodge0 = h0.asDiagonal(); + polygonHodge0Inverse = h0.asDiagonal().inverse(); + + // hodge1 (inner product matrix on 1-forms). + polygonHodge1 = Eigen::SparseMatrix(H, H); + Eigen::MatrixXd Mf; // local per-polygon matrix + std::vector hIndices; // indices of halfedges of polygon face + for (Face f : mesh.faces()) { + size_t n = f.degree(); + Mf = polygonPerFaceInnerProductMatrix(f); + hIndices.clear(); + for (Halfedge he : f.adjacentHalfedges()) hIndices.push_back(halfedgeIndices[he]); + for (size_t j = 0; j < n; j++) { + for (size_t i = 0; i < n; i++) { + tripletsH1.emplace_back(hIndices[i], hIndices[j], Mf(i, j)); + } + } + } + polygonHodge1.setFromTriplets(tripletsH1.begin(), tripletsH1.end()); + + // hodge2 + Eigen::VectorXd h2(F); + for (Face f : mesh.faces()) { + size_t fIdx = faceIndices[f]; + h2[fIdx] = 1. / faceAreas[f]; + } + polygonHodge2 = h2.asDiagonal(); + polygonHodge2Inverse = h2.asDiagonal().inverse(); +} +void EmbeddedGeometryInterface::requirePolygonDECOperators() { polygonDECOperatorsQ.require(); } +void EmbeddedGeometryInterface::unrequirePolygonDECOperators() { polygonDECOperatorsQ.unrequire(); } + + +// Helper functions + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonPerFaceLaplacian(const Face& f) { + Eigen::MatrixXd Df = polygonDerivativeMatrix(f); + return Df.transpose() * polygonPerFaceInnerProductMatrix(f) * Df; // build positive-definite +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonPerFaceInnerProductMatrix(const Face& f) { + faceAreasQ.ensureHave(); + + Eigen::MatrixXd Uf = polygonSharp(f); + Eigen::MatrixXd Pf = polygonProjectionMatrix(f); + double A = faceAreas[f]; + return A * Uf.transpose() * Uf + polygonLambda * Pf.transpose() * Pf; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonPerFaceConnectionLaplacian(const Face& f) { + faceAreasQ.ensureHave(); + + Eigen::MatrixXd G = polygonCovariantGradient(f); + Eigen::MatrixXd P = polygonCovariantProjection(f); + double A = faceAreas[f]; + Eigen::MatrixXd L = A * G.transpose() * G + polygonLambda * P.transpose() * P; + return L; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonBlockConnection(const Face& f) { + size_t d = f.degree(); + Eigen::MatrixXd R = Eigen::MatrixXd::Zero(2 * d, 2 * d); + size_t cpt = 0; + for (Vertex v : f.adjacentVertices()) { + Eigen::Matrix2d Rv = Rvf(v, f); + R.block<2, 2>(2 * cpt, 2 * cpt) = Rv; + cpt++; + } + return R; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonCovariantGradient(const Face& f) { + return kroneckerWithI2(Tf(f).transpose() * polygonPerFaceGradientMatrix(f)) * polygonBlockConnection(f); +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonCovariantProjection(const Face& f) { + return kroneckerWithI2(polygonProjectionMatrix(f) * polygonDerivativeMatrix(f)) * polygonBlockConnection(f); +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonProjectionMatrix(const Face& f) { + size_t d = f.degree(); + Eigen::MatrixXd P = Eigen::MatrixXd::Identity(d, d) - polygonFlat(f) * polygonSharp(f); + return P; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonPerFaceGradientMatrix(const Face& f) { + faceNormalsQ.ensureHave(); + faceAreasQ.ensureHave(); + + // equivalent to applying the (per-face) sharp operator to the (per-face) exterior derivative + double A = faceAreas[f]; + Vector3 n = faceNormals[f]; + Eigen::Vector3d N = {n[0], n[1], n[2]}; + return 1. / A * bracket(N) * polygonCoGradientMatrix(f); +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonCoGradientMatrix(const Face& f) { + return polygonEdgeVectorMatrix(f).transpose() * polygonAveragingMatrix(f); +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonAveragingMatrix(const Face& f) const { + size_t d = f.degree(); + Eigen::MatrixXd A = Eigen::MatrixXd::Zero(d, d); + for (size_t i = 0; i < d; i++) { + A(i, (i + 1) % d) = 0.5; + A(i, i) = 0.5; + } + return A; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonDerivativeMatrix(const Face& f) const { + size_t d = f.degree(); + Eigen::MatrixXd D = Eigen::MatrixXd::Zero(d, d); + for (size_t i = 0; i < d; i++) { + D(i, (i + 1) % d) = 1.; + D(i, i) = -1.; + } + return D; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonEdgeVectorMatrix(const Face& f) { + return polygonDerivativeMatrix(f) * polygonPositionMatrix(f); +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonEdgeMidpointMatrix(const Face& f) { + return polygonAveragingMatrix(f) * polygonPositionMatrix(f); +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonFlat(const Face& f) { + faceNormalsQ.ensureHave(); + + Vector3 n = faceNormals[f]; + Eigen::Vector3d N = {n[0], n[1], n[2]}; + return polygonEdgeVectorMatrix(f) * (Eigen::MatrixXd::Identity(3, 3) - N * N.transpose()); +} + +Eigen::MatrixXd EmbeddedGeometryInterface::polygonSharp(const Face& f) { + faceAreasQ.ensureHave(); + faceNormalsQ.ensureHave(); + + size_t d = f.degree(); + double A = faceAreas[f]; + Vector3 n = faceNormals[f]; + Eigen::Vector3d N = {n[0], n[1], n[2]}; + Eigen::Vector3d c = polygonCentroid(f); + return 1. / A * bracket(N) * (polygonEdgeMidpointMatrix(f).transpose() - c * Eigen::VectorXd::Ones(d).transpose()); +} + +void EmbeddedGeometryInterface::computePolygonVertexNormals() { + faceAreasQ.ensureHave(); + faceNormalsQ.ensureHave(); + + polygonVertexNormals = VertexData(mesh); + for (Vertex v : mesh.vertices()) { + Eigen::Vector3d vN(0., 0., 0.); + for (Face f : v.adjacentFaces()) { + Vector3 n = faceNormals[f]; + Eigen::Vector3d N = {n[0], n[1], n[2]}; + vN += N * faceAreas[f]; + } + vN /= vN.norm(); + polygonVertexNormals[v] = vN; + } +} + +Eigen::Vector3d EmbeddedGeometryInterface::polygonCentroid(const Face& f) { + vertexPositionsQ.ensureHave(); + + Vector3 c = {0, 0, 0}; + for (Vertex v : f.adjacentVertices()) { + c += vertexPositions[v]; + } + c /= f.degree(); + return Eigen::Vector3d(c[0], c[1], c[2]); +} + +Eigen::MatrixXd EmbeddedGeometryInterface::Tv(const Vertex& v) { + vertexTangentBasisQ.ensureHave(); + + // Return 3 x 2 matrix defining the tangent space at vertex v, with basis vectors in columns. + Vector3 xVec = vertexTangentBasis[v][0]; + Vector3 yVec = vertexTangentBasis[v][1]; + Eigen::Vector3d uu = {xVec[0], xVec[1], xVec[2]}; + Eigen::Vector3d vv = {yVec[0], yVec[1], yVec[2]}; + Eigen::MatrixXd B(3, 2); + B.col(0) = uu; + B.col(1) = vv; + return B; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::Tf(const Face& f) { + faceTangentBasisQ.ensureHave(); + + // Return 3 x 2 matrix defining the tangent space at face f, with basis vectors in columns. + Vector3 xVec = faceTangentBasis[f][0]; + Vector3 yVec = faceTangentBasis[f][1]; + Eigen::Vector3d uu = {xVec[0], xVec[1], xVec[2]}; + Eigen::Vector3d vv = {yVec[0], yVec[1], yVec[2]}; + Eigen::MatrixXd B(3, 2); + B.col(0) = uu; + B.col(1) = vv; + return B; +} + +Eigen::Matrix2d EmbeddedGeometryInterface::Rvf(const Vertex& v, const Face& f) { + return Tf(f).transpose() * Qvf(v, f) * Tv(v); +} + +Eigen::Matrix3d EmbeddedGeometryInterface::Qvf(const Vertex& v, const Face& f) { + polygonVertexNormalsQ.ensureHave(); + faceNormalsQ.ensureHave(); + + // Return 3 x 3 rotation matrix to align n_v to n_f. + Vector3 n = faceNormals[f]; + Eigen::Vector3d nf = {n[0], n[1], n[2]}; + Eigen::Vector3d nv = polygonVertexNormals[v]; + double c = nv.dot(nf); + + // Special case for opposite nv and nf vectors. + if (std::abs(c + 1.0) < 1e-5) return -Eigen::Matrix3d::Identity(); + + Eigen::Vector3d vv = nv.cross(nf); + Eigen::Matrix3d skew = bracket(vv); + return Eigen::Matrix3d::Identity() + skew + 1.0 / (1.0 + c) * skew * skew; +} + +Eigen::Matrix3d EmbeddedGeometryInterface::bracket(const Eigen::Vector3d& n) const { + Eigen::Matrix3d B; + B << 0., -n[2], n[1], n[2], 0., -n[0], -n[1], n[0], 0.; + return B; +} + +Eigen::Vector3d EmbeddedGeometryInterface::project(const Eigen::Vector3d& u, const Eigen::Vector3d& n) const { + // Project u on the orthgonal of n. + return u - (u.dot(n) / n.squaredNorm()) * n; +} + +Eigen::MatrixXd EmbeddedGeometryInterface::kroneckerWithI2(const Eigen::MatrixXd& M) const { + size_t h = M.rows(); + size_t w = M.cols(); + Eigen::MatrixXd MK = Eigen::MatrixXd::Zero(h * 2, w * 2); + for (size_t j = 0; j < h; j++) + for (size_t i = 0; i < w; i++) { + MK(2 * j, 2 * i) = M(j, i); + MK(2 * j + 1, 2 * i + 1) = M(j, i); + } + return MK; +} } // namespace surface } // namespace geometrycentral diff --git a/src/surface/polygon_mesh_heat_solver.cpp b/src/surface/polygon_mesh_heat_solver.cpp new file mode 100644 index 00000000..a75f5c0a --- /dev/null +++ b/src/surface/polygon_mesh_heat_solver.cpp @@ -0,0 +1,331 @@ +#include "geometrycentral/surface/polygon_mesh_heat_solver.h" + +namespace geometrycentral { +namespace surface { + +PolygonMeshHeatSolver::PolygonMeshHeatSolver(EmbeddedGeometryInterface& geom_, double tCoef_) + : tCoef(tCoef_), mesh(geom_.mesh), geom(geom_) + +{ + // Though de Goes et al. suggest using the maximum length over all polygon diagonals, using the mean length seems more + // accurate. + geom.requireEdgeLengths(); + double meanEdgeLength = 0.; + for (Edge e : mesh.edges()) meanEdgeLength += geom.edgeLengths[e]; + meanEdgeLength /= mesh.nEdges(); + shortTime = tCoef * meanEdgeLength * meanEdgeLength; + geom.unrequireEdgeLengths(); + + geom.requirePolygonVertexLumpedMassMatrix(); + geom.requirePolygonLaplacian(); + massMat = geom.polygonVertexLumpedMassMatrix; + laplaceMat = geom.polygonLaplacian; + geom.unrequirePolygonVertexLumpedMassMatrix(); + geom.unrequirePolygonLaplacian(); +} + +void PolygonMeshHeatSolver::ensureHaveScalarHeatSolver() { + if (scalarHeatSolver != nullptr) return; + + SparseMatrix heatOp = massMat + shortTime * laplaceMat; + scalarHeatSolver.reset(new PositiveDefiniteSolver(heatOp)); +} + +void PolygonMeshHeatSolver::ensureHaveVectorHeatSolver() { + if (vectorHeatSolver != nullptr) return; + + geom.requirePolygonVertexConnectionLaplacian(); + + SparseMatrix>& L = geom.polygonVertexConnectionLaplacian; + SparseMatrix> vectorOp = massMat.cast>() + shortTime * L; + + vectorHeatSolver.reset(new PositiveDefiniteSolver>(vectorOp)); + + geom.unrequirePolygonVertexConnectionLaplacian(); +} + + +void PolygonMeshHeatSolver::ensureHavePoissonSolver() { + if (poissonSolver != nullptr) return; + + poissonSolver.reset(new PositiveDefiniteSolver(laplaceMat)); +} + +VertexData PolygonMeshHeatSolver::computeDistance(const Vertex& sourceVert) { + std::vector v{sourceVert}; + return computeDistance(v); +} + +VertexData PolygonMeshHeatSolver::computeDistance(const std::vector& sourceVerts) { + GC_SAFETY_ASSERT(sourceVerts.size() != 0, "must have at least one source"); + + geom.requirePolygonGradientMatrix(); + geom.requirePolygonDivergenceMatrix(); + geom.requireVertexIndices(); + + // Flow heat. + size_t V = mesh.nVertices(); + size_t F = mesh.nFaces(); + Vector rho = Vector::Zero(V); + for (Vertex v : sourceVerts) { + size_t vIdx = geom.vertexIndices[v]; + rho[vIdx] += 1; + } + + ensureHaveScalarHeatSolver(); + Vector X = scalarHeatSolver->solve(rho); + + // Normalize gradient. + Vector Y = geom.polygonGradientMatrix * X; // 3|F| + for (size_t i = 0; i < F; i++) { + Vector3 g = {Y[3 * i], Y[3 * i + 1], Y[3 * i + 2]}; + g /= g.norm(); + for (int j = 0; j < 3; j++) Y[3 * i + j] = g[j]; + } + + // Integrate. + ensureHavePoissonSolver(); + Vector div = -geom.polygonDivergenceMatrix * Y; + Vector distances = poissonSolver->solve(div); + + // Shift solution. + distances -= distances.minCoeff() * Vector::Ones(V); + + geom.unrequirePolygonGradientMatrix(); + geom.unrequirePolygonDivergenceMatrix(); + geom.unrequireVertexIndices(); + + return VertexData(mesh, distances); +} + +VertexData PolygonMeshHeatSolver::extendScalars(const std::vector>& sources) { + GC_SAFETY_ASSERT(sources.size() != 0, "must have at least one source"); + + ensureHaveScalarHeatSolver(); + + size_t V = mesh.nVertices(); + Vector rhsVals = Vector::Zero(V); + Vector rhsOnes = Vector::Zero(V); + for (size_t i = 0; i < sources.size(); i++) { + size_t ind = std::get<0>(sources[i]).getIndex(); + double val = std::get<1>(sources[i]); + rhsVals[ind] = val; + rhsOnes[ind] = 1.; + } + + Vector interpVals = scalarHeatSolver->solve(rhsVals); + Vector interpOnes = scalarHeatSolver->solve(rhsOnes); + Vector resultArr = (interpVals.array() / interpOnes.array()); + + VertexData result(mesh, resultArr); + return result; +} + +VertexData PolygonMeshHeatSolver::transportTangentVector(const Vertex& sourceVert, + const Vector2& sourceVector) { + + return transportTangentVectors({std::make_tuple(sourceVert, sourceVector)}); +} + +VertexData +PolygonMeshHeatSolver::transportTangentVectors(const std::vector>& sources) { + GC_SAFETY_ASSERT(sources.size() != 0, "must have at least one source"); + + ensureHaveVectorHeatSolver(); + geom.requireVertexIndices(); + + // Construct rhs + size_t V = mesh.nVertices(); + Vector> dirRHS = Vector>::Zero(V); + std::vector> magnitudeSources; + bool normsAllSame = true; + double firstNorm = std::get<1>(sources[0]).norm(); + for (size_t i = 0; i < sources.size(); i++) { + Vertex v = std::get<0>(sources[i]); + size_t ind = geom.vertexIndices[v]; + Vector2 vec = std::get<1>(sources[i]); + dirRHS(ind) += std::complex(vec); + magnitudeSources.emplace_back(v, vec.norm()); + + // Check if all norms same + double thisNorm = vec.norm(); + if (std::abs(firstNorm - thisNorm) > std::fmax(firstNorm, thisNorm) * 1e-10) { + normsAllSame = false; + } + } + geom.unrequireVertexIndices(); + + // Transport + Vector> vecSolution = vectorHeatSolver->solve(dirRHS); + + // Set scale + VertexData solution(mesh); + if (normsAllSame) { + vecSolution = (vecSolution.array() / vecSolution.array().abs()) * firstNorm; + for (size_t i = 0; i < V; i++) solution[i] = Vector2::fromComplex(vecSolution[i]); + } else { + VertexData interpMags = extendScalars(magnitudeSources); + for (size_t i = 0; i < V; i++) { + Vector2 dir = Vector2::fromComplex(vecSolution[i]).normalize(); + solution[i] = dir * interpMags[i]; + } + } + + return solution; +} + +VertexData PolygonMeshHeatSolver::computeSignedDistance(const std::vector>& curves, + const LevelSetConstraint& levelSetConstraint) { + + size_t V = mesh.nVertices(); + Vector> X0 = Vector>::Zero(V); + geom.requireVertexIndices(); + geom.requireVertexTangentBasis(); + geom.requireVertexNormals(); + for (const auto& curve : curves) buildSignedCurveSource(curve, X0); + geom.unrequireVertexNormals(); + GC_SAFETY_ASSERT(X0.norm() != 0, "must have at least one source"); + + ensureHaveVectorHeatSolver(); + Vector> Xt = vectorHeatSolver->solve(X0); + + // Average onto faces, and normalize. + size_t F = mesh.nFaces(); + Vector Y(3 * F); + geom.requireFaceIndices(); + for (Face f : mesh.faces()) { + Vector3 Yf = {0, 0, 0}; + for (Vertex v : f.adjacentVertices()) { + size_t vIdx = geom.vertexIndices[v]; + Yf += std::real(Xt[vIdx]) * geom.vertexTangentBasis[v][0]; + Yf += std::imag(Xt[vIdx]) * geom.vertexTangentBasis[v][1]; + } + Yf /= Yf.norm(); + size_t fIdx = geom.faceIndices[f]; + for (int j = 0; j < 3; j++) { + Y[3 * fIdx + j] = Yf[j]; + } + } + geom.unrequireFaceIndices(); + geom.unrequireVertexTangentBasis(); + + geom.requirePolygonDivergenceMatrix(); + Vector divYt = geom.polygonDivergenceMatrix * Y; + geom.unrequirePolygonDivergenceMatrix(); + + Vector phi; + if (levelSetConstraint == LevelSetConstraint::None) { + ensureHavePoissonSolver(); + phi = poissonSolver->solve(divYt); + // Shift so that average value is zero along input curves. + double shift = computeAverageValue(curves, phi); + phi -= shift * Vector::Ones(V); + } else if (levelSetConstraint == LevelSetConstraint::ZeroSet) { + Vector setAMembership = Vector::Ones(V); + for (const auto& curve : curves) { + for (const Vertex& v : curve) { + setAMembership[geom.vertexIndices[v]] = false; + } + } + int nB = V - setAMembership.cast().sum(); + Vector bcVals = Vector::Zero(nB); + BlockDecompositionResult decomp = blockDecomposeSquare(laplaceMat, setAMembership, true); + Vector rhsValsA, rhsValsB; + decomposeVector(decomp, divYt, rhsValsA, rhsValsB); + Vector combinedRHS = rhsValsA; + Vector Aresult = solvePositiveDefinite(decomp.AA, combinedRHS); + phi = reassembleVector(decomp, Aresult, bcVals); + } else if (levelSetConstraint == LevelSetConstraint::Multiple) { + std::vector> triplets; + SparseMatrix A; + size_t m = 0; + for (const auto& curve : curves) { + size_t nNodes = curve.size(); + if (nNodes < 2) continue; + size_t v0 = geom.vertexIndices[curve[0]]; + for (size_t i = 1; i < nNodes; i++) { + size_t vIdx = geom.vertexIndices[curve[i]]; + triplets.emplace_back(m, v0, 1); + triplets.emplace_back(m, vIdx, -1); + m++; + } + } + A.resize(m, V); + A.setFromTriplets(triplets.begin(), triplets.end()); + SparseMatrix Z(m, m); + SparseMatrix LHS1 = horizontalStack({laplaceMat, A.transpose()}); + SparseMatrix LHS2 = horizontalStack({A, Z}); + SparseMatrix LHS = verticalStack({LHS1, LHS2}); + Vector RHS = Vector::Zero(V + m); + RHS.head(V) = divYt; + Vector soln = solveSquare(LHS, RHS); + phi = soln.head(V); + double shift = computeAverageValue(curves, phi); + phi -= shift * Vector::Ones(V); + } + + geom.unrequireVertexIndices(); + + return VertexData(mesh, phi); +} + +void PolygonMeshHeatSolver::buildSignedCurveSource(const std::vector& curve, + Vector>& X0) const { + + // Encode curve input by expressing curve normals in vertex tangent spaces. + size_t nNodes = curve.size(); + for (size_t i = 0; i < nNodes - 1; i++) { + Vertex vA = curve[i]; + Vertex vB = curve[i + 1]; + size_t vIdxA = geom.vertexIndices[vA]; + size_t vIdxB = geom.vertexIndices[vB]; + Vector3 pA = geom.vertexPositions[vA]; + Vector3 pB = geom.vertexPositions[vB]; + Vector3 tangent = pB - pA; + Vector3 vnA = geom.vertexNormals[vA]; + Vector3 vnB = geom.vertexNormals[vB]; + Vector3 xAxisA = geom.vertexTangentBasis[vA][0]; + Vector3 yAxisA = geom.vertexTangentBasis[vA][1]; + Vector3 xAxisB = geom.vertexTangentBasis[vB][0]; + Vector3 yAxisB = geom.vertexTangentBasis[vB][1]; + Vector3 tangentA = dot(xAxisA, tangent) * xAxisA + dot(yAxisA, tangent) * yAxisA; + Vector3 tangentB = dot(xAxisB, tangent) * xAxisB + dot(yAxisB, tangent) * yAxisB; + Vector3 normalA = cross(vnA, tangentA); + Vector3 normalB = cross(vnB, tangentB); + X0[vIdxA] += std::complex(dot(geom.vertexTangentBasis[vA][0], normalA), + dot(geom.vertexTangentBasis[vA][1], normalA)); + X0[vIdxB] += std::complex(dot(geom.vertexTangentBasis[vB][0], normalB), + dot(geom.vertexTangentBasis[vB][1], normalB)); + } +} + +double PolygonMeshHeatSolver::computeAverageValue(const std::vector>& curves, + const Vector& u) { + + geom.requireVertexIndices(); + + double shift = 0.; + double totalLength = 0.; + for (const auto& curve : curves) { + size_t nNodes = curve.size(); + for (size_t i = 0; i < nNodes - 1; i++) { + Vertex vA = curve[i]; + Vertex vB = curve[i + 1]; + size_t vIdxA = geom.vertexIndices[vA]; + size_t vIdxB = geom.vertexIndices[vB]; + Vector3 pA = geom.vertexPositions[vA]; + Vector3 pB = geom.vertexPositions[vB]; + double length = (pB - pA).norm(); + shift += 0.5 * length * (u[vIdxA] + u[vIdxB]); + totalLength += length; + } + } + shift /= totalLength; + + geom.unrequireVertexIndices(); + + return shift; +} + +} // namespace surface +} // namespace geometrycentral \ No newline at end of file diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 3f13bf1f..12869f92 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -106,6 +106,7 @@ set(TEST_SRCS src/load_test_meshes.cpp src/point_cloud_test.cpp src/poisson_disk_sampler_test.cpp + src/polygon_operators_test.cpp src/stl_reader_test.cpp src/surface_misc_test.cpp ) diff --git a/test/src/polygon_operators_test.cpp b/test/src/polygon_operators_test.cpp new file mode 100644 index 00000000..a9e9756c --- /dev/null +++ b/test/src/polygon_operators_test.cpp @@ -0,0 +1,116 @@ +#include "geometrycentral/surface/embedded_geometry_interface.h" +#include "geometrycentral/surface/meshio.h" +#include "geometrycentral/surface/vertex_position_geometry.h" + +#include "load_test_meshes.h" + +#include "gtest/gtest.h" + +#include +#include + +using namespace geometrycentral; +using namespace geometrycentral::surface; + +class PolygonMeshSuite : public MeshAssetSuite {}; + +/* Test that the lumped mass matrices correspond with their unlumped versions. */ +TEST_F(PolygonMeshSuite, MassLumpingTest) { + + double epsilon = 1e-8; + for (MeshAsset& a : allMeshes()) { + a.printThyName(); + SurfaceMesh& mesh = *a.mesh; + VertexPositionGeometry& geometry = *a.geometry; + + geometry.requireSimplePolygonVertexLumpedMassMatrix(); + geometry.requireSimplePolygonVertexGalerkinMassMatrix(); + + SparseMatrix Mt = geometry.simplePolygonVertexGalerkinMassMatrix.transpose(); + SparseMatrix& M = geometry.simplePolygonVertexLumpedMassMatrix; + for (size_t i = 0; i < mesh.nVertices(); i++) { + double rowSum = 0.; + for (SparseMatrix::InnerIterator it(Mt, i); it; ++it) { + rowSum += it.value(); + } + EXPECT_LT(std::abs(rowSum - M.coeffRef(i, i)), epsilon); + } + + geometry.unrequireSimplePolygonVertexLumpedMassMatrix(); + geometry.unrequireSimplePolygonVertexGalerkinMassMatrix(); + } +} + +/* Check that polygon Laplacians and mass matrices coincide with simplicial versions on a triangle mesh. */ +TEST_F(PolygonMeshSuite, TriangularTest) { + + double epsilon = 1e-8; + for (MeshAsset& a : triangularMeshes()) { + a.printThyName(); + SurfaceMesh& mesh = *a.mesh; + VertexPositionGeometry& geometry = *a.geometry; + + geometry.requireVertexGalerkinMassMatrix(); + geometry.requireVertexLumpedMassMatrix(); + geometry.requireCotanLaplacian(); + geometry.requireSimplePolygonLaplacian(); + geometry.requireSimplePolygonVertexLumpedMassMatrix(); + geometry.requireSimplePolygonVertexGalerkinMassMatrix(); + geometry.requirePolygonLaplacian(); + geometry.requirePolygonVertexLumpedMassMatrix(); + + double L = geometry.cotanLaplacian.norm(); + double Mg = geometry.vertexGalerkinMassMatrix.norm(); + double Ml = geometry.vertexLumpedMassMatrix.norm(); + + EXPECT_LT((geometry.simplePolygonLaplacian - geometry.cotanLaplacian).norm() / L, epsilon); + EXPECT_LT((geometry.polygonLaplacian - geometry.cotanLaplacian).norm() / L, epsilon); + EXPECT_LT((geometry.simplePolygonVertexGalerkinMassMatrix - geometry.vertexGalerkinMassMatrix).norm() / Mg, + epsilon); + EXPECT_LT((geometry.simplePolygonVertexLumpedMassMatrix - geometry.vertexLumpedMassMatrix).norm() / Ml, epsilon); + EXPECT_LT((geometry.polygonVertexLumpedMassMatrix - geometry.vertexLumpedMassMatrix).norm() / Ml, epsilon); + + geometry.unrequireVertexGalerkinMassMatrix(); + geometry.unrequireVertexLumpedMassMatrix(); + geometry.unrequireCotanLaplacian(); + geometry.unrequireSimplePolygonLaplacian(); + geometry.unrequireSimplePolygonVertexLumpedMassMatrix(); + geometry.unrequireSimplePolygonVertexGalerkinMassMatrix(); + geometry.unrequirePolygonLaplacian(); + geometry.unrequirePolygonVertexLumpedMassMatrix(); + } +} + +/* Check that Laplacian can be assembled as expected from DEC operators. */ +TEST_F(PolygonMeshSuite, DECTest) { + + double epsilon = 1e-8; + for (MeshAsset& a : allMeshes()) { + a.printThyName(); + SurfaceMesh& mesh = *a.mesh; + VertexPositionGeometry& geometry = *a.geometry; + + geometry.requirePolygonDECOperators(); + geometry.requirePolygonLaplacian(); + + SparseMatrix& L = geometry.polygonLaplacian; + SparseMatrix& h0 = geometry.polygonHodge0; + SparseMatrix& h0Inv = geometry.polygonHodge0Inverse; + SparseMatrix& h1 = geometry.polygonHodge1; + SparseMatrix& h2 = geometry.polygonHodge2; + SparseMatrix& h2Inv = geometry.polygonHodge2Inverse; + SparseMatrix& d0 = geometry.polygonD0; + SparseMatrix& d1 = geometry.polygonD1; + EXPECT_LT((L - d0.transpose() * h1 * d0).norm(), epsilon); + + if (mesh.isTriangular()) { + geometry.requireDECOperators(); + EXPECT_LT((h0 - geometry.hodge0).norm(), epsilon); + EXPECT_LT((h2 - geometry.hodge2).norm(), epsilon); + geometry.unrequireDECOperators(); + } + + geometry.unrequirePolygonDECOperators(); + geometry.unrequirePolygonLaplacian(); + } +} \ No newline at end of file