diff --git a/docs/docs/surface/geometry/quantities.md b/docs/docs/surface/geometry/quantities.md index aeed7489..b5c287ca 100644 --- a/docs/docs/surface/geometry/quantities.md +++ b/docs/docs/surface/geometry/quantities.md @@ -547,15 +547,13 @@ In graphics and geometry processing, Crouzeix-Raviart elements have been used, f 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!) +The polygon operators are 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; the code is based on [Astrid Bunge and Mario Botsch's implementation of their paper](https://github.com/mbotsch/polygon-laplacian). The method builds local operators whose matrices are assembled per-polygon, so they will work 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)" +??? func "polygon mesh Laplacian" - ##### polygon mesh Laplacian (simple) + ##### polygon mesh Laplacian The discrete Laplace operator acting on polygon meshes, using Bunge et al.'s virtual refinement method in _Polygon Laplacian Made Simple_. @@ -568,9 +566,9 @@ Here are polygon mesh operators from [Bunge et al.'s _Polygon Laplacian Made Sim - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::simplePolygonLaplacian` - **require:** `void EmbeddedGeometryInterface::requireSimplePolygonLaplacian()` -??? func "polygon mesh vertex lumped mass matrix (simple)" +??? func "polygon mesh vertex lumped mass matrix" - ##### polygon mesh vertex lumped mass matrix (simple) + ##### polygon mesh vertex lumped mass matrix 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. @@ -579,9 +577,9 @@ Here are polygon mesh operators from [Bunge et al.'s _Polygon Laplacian Made Sim - **member:** `Eigen::SparseMatrix EmbeddedGeometryInterface::simplePolygonVertexLumpedMassMatrix` - **require:** `void EmbeddedGeometryInterface::requireSimplePolygonVertexLumpedMassMatrix()` -??? func "polygon mesh vertex Galerkin mass matrix (simple)" +??? func "polygon mesh vertex Galerkin mass matrix" - ##### polygon mesh vertex Galerkin mass matrix (simple) + ##### polygon mesh vertex Galerkin mass matrix A $|V| \times |V|$ real matrix, using Bunge et al.'s virtual refinement method in _Polygon Laplacian Made Simple_. @@ -590,100 +588,23 @@ Here are polygon mesh operators from [Bunge et al.'s _Polygon Laplacian Made Sim - **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 discrete connection Laplacian operator, which applies to vector fields defined in vertex tangent spaces; based on the Laplacian from Bunge et al.'s _Polygon Laplacian Made Simple_. + + 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()` + - **member:** `Eigen::SparseMatrix> EmbeddedGeometryInterface::simplePolygonVertexConnectionLaplacian` + - **require:** `void EmbeddedGeometryInterface::requireSimplePolygonVertexConnectionLaplacian()` ## Extrinsic angles diff --git a/include/geometrycentral/surface/embedded_geometry_interface.h b/include/geometrycentral/surface/embedded_geometry_interface.h index 64d2710d..d85ecc41 100644 --- a/include/geometrycentral/surface/embedded_geometry_interface.h +++ b/include/geometrycentral/surface/embedded_geometry_interface.h @@ -1,7 +1,6 @@ #pragma once #include "geometrycentral/surface/extrinsic_geometry_interface.h" -#include "geometrycentral/surface/polygon_mesh_helpers.h" #include "geometrycentral/surface/surface_mesh.h" #include "geometrycentral/utilities/vector3.h" @@ -58,12 +57,34 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface { // == 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 Eigen::SparseMatrix simplePolygonLaplacian; void requireSimplePolygonLaplacian(); void unrequireSimplePolygonLaplacian(); + // Divergence + Eigen::SparseMatrix simplePolygonDivergenceMatrix; + void requireSimplePolygonDivergenceMatrix(); + void unrequireSimplePolygonDivergenceMatrix(); + + // Gradient + Eigen::SparseMatrix simplePolygonGradientMatrix; + void requireSimplePolygonGradientMatrix(); + void unrequireSimplePolygonGradientMatrix(); + + Eigen::SparseMatrix simplePolygonProlongationMatrix; + void requireSimplePolygonProlongationMatrix(); + void unrequireSimplePolygonProlongationMatrix(); + + // connection Laplacian + Eigen::SparseMatrix> simplePolygonVertexConnectionLaplacian; + void requireSimplePolygonVertexConnectionLaplacian(); + void unrequireSimplePolygonVertexConnectionLaplacian(); + // Vertex Galerkin mass matrix (unlumped) Eigen::SparseMatrix simplePolygonVertexGalerkinMassMatrix; void requireSimplePolygonVertexGalerkinMassMatrix(); @@ -74,38 +95,6 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface { 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; @@ -148,6 +137,22 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface { DependentQuantityD> simplePolygonLaplacianQ; virtual void computeSimplePolygonLaplacian(); + // Divergence + DependentQuantityD> simplePolygonDivergenceMatrixQ; + virtual void computeSimplePolygonDivergenceMatrix(); + + // Gradient + DependentQuantityD> simplePolygonGradientMatrixQ; + virtual void computeSimplePolygonGradientMatrix(); + + // Prolongation + DependentQuantityD> simplePolygonProlongationMatrixQ; + virtual void computeSimplePolygonProlongationMatrix(); + + // Connection Laplacian + DependentQuantityD>> simplePolygonVertexConnectionLaplacianQ; + virtual void computeSimplePolygonVertexConnectionLaplacian(); + // Vertex mass matrix (unlumped) DependentQuantityD> simplePolygonVertexGalerkinMassMatrixQ; virtual void computeSimplePolygonVertexGalerkinMassMatrix(); @@ -156,70 +161,22 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface { DependentQuantityD> simplePolygonVertexLumpedMassMatrixQ; virtual void computeSimplePolygonVertexLumpedMassMatrix(); + // helper functions + Eigen::VectorXd simplePolygonVirtualVertex(const Eigen::MatrixXd& poly); + Eigen::Vector3d gradientHatFunction(const Eigen::Vector3d& a, const Eigen::Vector3d& b, + const Eigen::Vector3d& c) const; + // helper functions -- these all depend on quantities in EmbeddedGeometryInterface, which makes them hard to separate // their declarations into a separate file. FaceData virtualRefinementAreaWeights; + FaceData virtualRefinementAreaPoints; DependentQuantityD> virtualRefinementAreaWeightsQ; // affine weights for each virtual node + DependentQuantityD> virtualRefinementAreaPointsQ; virtual void computeVirtualRefinementAreaWeights(); virtual Eigen::MatrixXd simplePolygonMassMatrix(const Face& f); virtual Eigen::MatrixXd simplePolygonStiffnessMatrix(const Face& f); - virtual SparseMatrix simplePolygonProlongationMatrix(); + virtual SparseMatrix simplePolygonGradientMassMatrix(); virtual Eigen::MatrixXd polygonPositionMatrix(const Face& f); - - // = 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(); - - const double polygonLambda = 1.0; - VertexData polygonVertexNormals; - DependentQuantityD> polygonVertexNormalsQ; - // helper functions -- these all depend on quantities in EmbeddedGeometryInterface, which makes them hard to separate - // their declarations into a separate file. - 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 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); }; diff --git a/include/geometrycentral/surface/polygon_mesh_helpers.h b/include/geometrycentral/surface/polygon_mesh_helpers.h deleted file mode 100644 index 2e9b91ff..00000000 --- a/include/geometrycentral/surface/polygon_mesh_helpers.h +++ /dev/null @@ -1,33 +0,0 @@ -#pragma once - -#include "geometrycentral/surface/embedded_geometry_interface.h" - -namespace geometrycentral { -namespace surface { - -// Helper functions for Bunge et al. "Polygon Laplacian Made Simple" (2020). -// 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) - -Eigen::VectorXd simplePolygonVirtualVertex(const Eigen::MatrixXd& poly); - -// Helper functions for de Goes et al. "Discrete Differential Operators on Polygonal Meshes" (2020). -// 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) - -Eigen::MatrixXd polygonAveragingMatrix(const Face& f); - -Eigen::MatrixXd polygonDerivativeMatrix(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. -Eigen::Matrix3d bracket(const Eigen::Vector3d& n); - -Eigen::Vector3d project(const Eigen::Vector3d& u, const Eigen::Vector3d& n); - -Eigen::MatrixXd kroneckerWithI2(const Eigen::MatrixXd& M); - -} // namespace surface -} // namespace geometrycentral \ No newline at end of file diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index b74e243d..d3c0f4bf 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -33,7 +33,6 @@ SET(SRCS surface/heat_method_distance.cpp surface/vector_heat_method.cpp surface/signed_heat_method.cpp - surface/polygon_mesh_helpers.cpp surface/polygon_mesh_heat_solver.cpp surface/geodesic_centroidal_voronoi_tessellation.cpp surface/trace_geodesic.cpp diff --git a/src/surface/embedded_geometry_interface.cpp b/src/surface/embedded_geometry_interface.cpp index a32035b3..e30e280e 100644 --- a/src/surface/embedded_geometry_interface.cpp +++ b/src/surface/embedded_geometry_interface.cpp @@ -19,20 +19,16 @@ EmbeddedGeometryInterface::EmbeddedGeometryInterface(SurfaceMesh& mesh_) : vertexTangentBasisQ (&vertexTangentBasis, std::bind(&EmbeddedGeometryInterface::computeVertexTangentBasis, 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) + simplePolygonLaplacianQ (&simplePolygonLaplacian, std::bind(&EmbeddedGeometryInterface::computeSimplePolygonLaplacian, this), quantities), + simplePolygonDivergenceMatrixQ (&simplePolygonDivergenceMatrix, std::bind(&EmbeddedGeometryInterface::computeSimplePolygonDivergenceMatrix, this), quantities), + simplePolygonGradientMatrixQ (&simplePolygonGradientMatrix, std::bind(&EmbeddedGeometryInterface::computeSimplePolygonGradientMatrix, this), quantities), + simplePolygonProlongationMatrixQ (&simplePolygonProlongationMatrix, std::bind(&EmbeddedGeometryInterface::computeSimplePolygonProlongationMatrix, this), quantities), + simplePolygonVertexConnectionLaplacianQ (&simplePolygonVertexConnectionLaplacian, std::bind(&EmbeddedGeometryInterface::computeSimplePolygonVertexConnectionLaplacian, 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), + virtualRefinementAreaPointsQ (&virtualRefinementAreaPoints, std::bind(&EmbeddedGeometryInterface::computeVirtualRefinementAreaWeights, this), quantities) + {} // clang-format on @@ -389,6 +385,125 @@ void EmbeddedGeometryInterface::computeSimplePolygonLaplacian() { void EmbeddedGeometryInterface::requireSimplePolygonLaplacian() { simplePolygonLaplacianQ.require(); } void EmbeddedGeometryInterface::unrequireSimplePolygonLaplacian() { simplePolygonLaplacianQ.unrequire(); } +void EmbeddedGeometryInterface::computeSimplePolygonDivergenceMatrix() { + + simplePolygonGradientMatrixQ.ensureHave(); + + SparseMatrix M = simplePolygonGradientMassMatrix(); + simplePolygonDivergenceMatrix = -simplePolygonGradientMatrix.transpose() * M; +} +void EmbeddedGeometryInterface::requireSimplePolygonDivergenceMatrix() { simplePolygonDivergenceMatrixQ.require(); } +void EmbeddedGeometryInterface::unrequireSimplePolygonDivergenceMatrix() { simplePolygonDivergenceMatrixQ.unrequire(); } + +// 3|T^f| x |V| gradient matrix +void EmbeddedGeometryInterface::computeSimplePolygonGradientMatrix() { + + faceIndicesQ.ensureHave(); + vertexIndicesQ.ensureHave(); + vertexPositionsQ.ensureHave(); + simplePolygonProlongationMatrixQ.ensureHave(); + + size_t V = mesh.nVertices(); + size_t F = mesh.nFaces(); + SparseMatrix G; + std::vector> triplets; + int nTriangles = 0; + int k = 0; + for (Face f : mesh.faces()) { + size_t fIdx = faceIndices[f]; + nTriangles += f.degree(); + Eigen::Vector3d p = virtualRefinementAreaPoints[f]; + for (Halfedge he : f.adjacentHalfedges()) { + size_t v0 = vertexIndices[he.tailVertex()]; + size_t v1 = vertexIndices[he.tipVertex()]; + const Vector3& pos0 = vertexPositions[v0]; + const Vector3& pos1 = vertexPositions[v1]; + Eigen::Vector3d p0 = {pos0[0], pos0[1], pos0[2]}; + Eigen::Vector3d p1 = {pos1[0], pos1[1], pos1[2]}; + Eigen::Vector3d grad_p = gradientHatFunction(p, p0, p1); + Eigen::Vector3d grad_p0 = gradientHatFunction(p0, p1, p); + Eigen::Vector3d grad_p1 = gradientHatFunction(p1, p, p0); + for (int j = 0; j < 3; j++) { + triplets.emplace_back(3 * k + j, V + fIdx, grad_p(j)); + triplets.emplace_back(3 * k + j, v0, grad_p0(j)); + triplets.emplace_back(3 * k + j, v1, grad_p1(j)); + } + k++; + } + } + G.resize(3 * nTriangles, V + F); + G.setFromTriplets(triplets.begin(), triplets.end()); + simplePolygonGradientMatrix = G * simplePolygonProlongationMatrix; +} +void EmbeddedGeometryInterface::requireSimplePolygonGradientMatrix() { simplePolygonGradientMatrixQ.require(); } +void EmbeddedGeometryInterface::unrequireSimplePolygonGradientMatrix() { simplePolygonGradientMatrixQ.unrequire(); } + +// Vertex Connection Laplacian. Formed by augmenting a scalar Laplacian with rotations between each pair of vertices +// within each polygon. The rotations are based on extrinsic (esimated) tangent planes. +void EmbeddedGeometryInterface::computeSimplePolygonVertexConnectionLaplacian() { + + simplePolygonLaplacianQ.ensureHave(); + vertexNormalsQ.ensureHave(); + vertexTangentBasisQ.ensureHave(); + vertexIndicesQ.ensureHave(); + + // Compute rotation between tangent plane at vertex i and at vertex j. + // Rotations are not just along halfedges, but between any pair of vertices sharing a face. + std::map, std::complex> rotations; // entry at ij = rotation from i to j + for (Face f : mesh.faces()) { + for (Halfedge he : f.adjacentHalfedges()) { + Vertex vi = he.tailVertex(); + Vertex vj = he.tipVertex(); + Vector3 ni = vertexNormals[vi]; + Vector3 nj = vertexNormals[vj]; + Vector3 axis = cross(ni, nj); + Vector3 ei = vertexTangentBasis[vi][0]; + Vector3 ej = vertexTangentBasis[vj][0]; + Vector3 R_ei = ei; + if (axis.norm() > 1e-8) { + double angle = angleInPlane(ni, nj, axis); + R_ei = ei.rotateAround(axis, angle); + } + std::complex r_ij = std::complex(Vector2::fromAngle(angleInPlane(R_ei, ej, nj))); + + size_t viIdx = vertexIndices[vi]; + size_t vjIdx = vertexIndices[vj]; + if (viIdx > vjIdx) { + viIdx = vertexIndices[vj]; + vjIdx = vertexIndices[vi]; + r_ij = std::conj(r_ij); + } + std::pair key = std::make_pair(viIdx, vjIdx); + rotations[key] = r_ij; + } + } + + size_t V = mesh.nVertices(); + const SparseMatrix& L = simplePolygonLaplacian; + simplePolygonVertexConnectionLaplacian = Eigen::SparseMatrix>(V, V); + std::vector>> triplets; + for (int k = 0; k < L.outerSize(); ++k) { + for (SparseMatrix::InnerIterator it(L, k); it; ++it) { + size_t i = it.row(); + size_t j = it.col(); + if (i == j) { + triplets.emplace_back(it.row(), it.col(), it.value()); + continue; + } + bool sgn = (i < j); + std::pair key = sgn ? std::make_pair(i, j) : std::make_pair(j, i); + std::complex rot = sgn ? rotations[key] : std::conj(rotations[key]); + triplets.emplace_back(it.row(), it.col(), it.value() * rot); + } + } + simplePolygonVertexConnectionLaplacian.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requireSimplePolygonVertexConnectionLaplacian() { + simplePolygonVertexConnectionLaplacianQ.require(); +} +void EmbeddedGeometryInterface::unrequireSimplePolygonVertexConnectionLaplacian() { + simplePolygonVertexConnectionLaplacianQ.unrequire(); +} // Vertex Galerkin mass matrix (unlumped) void EmbeddedGeometryInterface::computeSimplePolygonVertexGalerkinMassMatrix() { @@ -455,6 +570,110 @@ void EmbeddedGeometryInterface::unrequireSimplePolygonVertexLumpedMassMatrix() { // Helper functions +Eigen::VectorXd EmbeddedGeometryInterface::simplePolygonVirtualVertex(const Eigen::MatrixXd& poly) { + + // 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; +} + +/* + * Block diagonal matrix whose i-th block consists of the 3×3 identity matrix multiplied by the area of the i-th + * triangle. + */ +SparseMatrix EmbeddedGeometryInterface::simplePolygonGradientMassMatrix() { + + virtualRefinementAreaWeightsQ.ensureHave(); + vertexPositionsQ.ensureHave(); + + SparseMatrix M; + std::vector> triplets; + int c = 0; + for (Face f : mesh.faces()) { + Eigen::Vector3d areaPoint = virtualRefinementAreaPoints[f]; + Vector3 ap = {areaPoint[0], areaPoint[1], areaPoint[2]}; + int i = 0; + for (Halfedge he : f.adjacentHalfedges()) { + Vector3 p0 = vertexPositions[he.tailVertex()]; + Vector3 p1 = vertexPositions[he.tipVertex()]; + double area = 0.5 * cross(p0 - ap, p1 - ap).norm(); + for (int j = 0; j < 3; j++) { + int idx = c + 3 * i + j; + triplets.emplace_back(idx, idx, area); + } + i++; + } + c += f.degree() * 3; + } + M.resize(c, c); + M.setFromTriplets(triplets.begin(), triplets.end()); + return M; +} + Eigen::MatrixXd EmbeddedGeometryInterface::simplePolygonMassMatrix(const Face& f) { virtualRefinementAreaWeightsQ.ensureHave(); @@ -462,7 +681,7 @@ Eigen::MatrixXd EmbeddedGeometryInterface::simplePolygonMassMatrix(const Face& f 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::Vector3d virtualVertex = virtualRefinementAreaPoints[f]; Eigen::VectorXd ln = Eigen::VectorXd::Zero(n + 1); double l[3], l2[3]; // lengths, lengths squared // Build triangle fan mass and cotan matrices @@ -501,7 +720,7 @@ Eigen::MatrixXd EmbeddedGeometryInterface::simplePolygonStiffnessMatrix(const Fa 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::Vector3d virtualVertex = virtualRefinementAreaPoints[f]; Eigen::VectorXd ln = Eigen::VectorXd::Zero(n + 1); double l[3], l2[3]; // lengths, lengths squared // Build triangle fan mass and cotan matrices @@ -540,14 +759,14 @@ Eigen::MatrixXd EmbeddedGeometryInterface::simplePolygonStiffnessMatrix(const Fa return S; } -SparseMatrix EmbeddedGeometryInterface::simplePolygonProlongationMatrix() { +void EmbeddedGeometryInterface::computeSimplePolygonProlongationMatrix() { virtualRefinementAreaWeightsQ.ensureHave(); vertexIndicesQ.ensureHave(); size_t V = mesh.nVertices(); size_t F = mesh.nFaces(); std::vector> triplets; - SparseMatrix P(V + F, V); + simplePolygonProlongationMatrix.resize(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()) { @@ -560,18 +779,23 @@ SparseMatrix EmbeddedGeometryInterface::simplePolygonProlongationMatrix( } j++; } - P.setFromTriplets(triplets.begin(), triplets.end()); - return P; + simplePolygonProlongationMatrix.setFromTriplets(triplets.begin(), triplets.end()); +} +void EmbeddedGeometryInterface::requireSimplePolygonProlongationMatrix() { simplePolygonProlongationMatrixQ.require(); } +void EmbeddedGeometryInterface::unrequireSimplePolygonProlongationMatrix() { + simplePolygonProlongationMatrixQ.unrequire(); } void EmbeddedGeometryInterface::computeVirtualRefinementAreaWeights() { vertexPositionsQ.ensureHave(); virtualRefinementAreaWeights = FaceData(mesh); + virtualRefinementAreaPoints = FaceData(mesh); for (Face f : mesh.faces()) { Eigen::MatrixXd poly = polygonPositionMatrix(f); Eigen::VectorXd weights = simplePolygonVirtualVertex(poly); + virtualRefinementAreaPoints[f] = poly.transpose() * weights; virtualRefinementAreaWeights[f] = weights; } } @@ -590,409 +814,18 @@ Eigen::MatrixXd EmbeddedGeometryInterface::polygonPositionMatrix(const Face& f) return poly; } +Eigen::Vector3d EmbeddedGeometryInterface::gradientHatFunction(const Eigen::Vector3d& a, const Eigen::Vector3d& b, + const Eigen::Vector3d& c) 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 -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::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::Vector3d gradient; + Eigen::Vector3d site = a - b; + Eigen::Vector3d base = c - b; + double area = 0.5 * (site.cross(base)).norm(); + double baseNorm = base.norm(); + Eigen::Vector3d grad = site - (site.dot(base) / baseNorm) * base / baseNorm; + grad = baseNorm * grad / grad.norm(); + gradient = Eigen::Vector3d(grad[0], grad[1], grad[2]) / (2.0 * area); + return gradient; } } // namespace surface diff --git a/src/surface/polygon_mesh_heat_solver.cpp b/src/surface/polygon_mesh_heat_solver.cpp index a75f5c0a..86d26813 100644 --- a/src/surface/polygon_mesh_heat_solver.cpp +++ b/src/surface/polygon_mesh_heat_solver.cpp @@ -16,12 +16,12 @@ PolygonMeshHeatSolver::PolygonMeshHeatSolver(EmbeddedGeometryInterface& geom_, d shortTime = tCoef * meanEdgeLength * meanEdgeLength; geom.unrequireEdgeLengths(); - geom.requirePolygonVertexLumpedMassMatrix(); - geom.requirePolygonLaplacian(); - massMat = geom.polygonVertexLumpedMassMatrix; - laplaceMat = geom.polygonLaplacian; - geom.unrequirePolygonVertexLumpedMassMatrix(); - geom.unrequirePolygonLaplacian(); + geom.requireSimplePolygonVertexLumpedMassMatrix(); + geom.requireSimplePolygonLaplacian(); + massMat = geom.simplePolygonVertexLumpedMassMatrix; + laplaceMat = geom.simplePolygonLaplacian; + geom.unrequireSimplePolygonVertexLumpedMassMatrix(); + geom.unrequireSimplePolygonLaplacian(); } void PolygonMeshHeatSolver::ensureHaveScalarHeatSolver() { @@ -34,14 +34,14 @@ void PolygonMeshHeatSolver::ensureHaveScalarHeatSolver() { void PolygonMeshHeatSolver::ensureHaveVectorHeatSolver() { if (vectorHeatSolver != nullptr) return; - geom.requirePolygonVertexConnectionLaplacian(); + geom.requireSimplePolygonVertexConnectionLaplacian(); - SparseMatrix>& L = geom.polygonVertexConnectionLaplacian; + SparseMatrix>& L = geom.simplePolygonVertexConnectionLaplacian; SparseMatrix> vectorOp = massMat.cast>() + shortTime * L; vectorHeatSolver.reset(new PositiveDefiniteSolver>(vectorOp)); - geom.unrequirePolygonVertexConnectionLaplacian(); + geom.unrequireSimplePolygonVertexConnectionLaplacian(); } @@ -59,8 +59,8 @@ VertexData PolygonMeshHeatSolver::computeDistance(const Vertex& sourceVe VertexData PolygonMeshHeatSolver::computeDistance(const std::vector& sourceVerts) { GC_SAFETY_ASSERT(sourceVerts.size() != 0, "must have at least one source"); - geom.requirePolygonGradientMatrix(); - geom.requirePolygonDivergenceMatrix(); + geom.requireSimplePolygonGradientMatrix(); + geom.requireSimplePolygonDivergenceMatrix(); geom.requireVertexIndices(); // Flow heat. @@ -76,23 +76,24 @@ VertexData PolygonMeshHeatSolver::computeDistance(const std::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]}; + Vector Y = + geom.simplePolygonGradientMatrix * X; // 3|T^f| x 1 (where |T^f| = # of triangles in the refinement) + for (int i = 0; i < Y.rows(); i += 3) { + Vector3 g = {Y[i], Y[i + 1], Y[i + 2]}; g /= g.norm(); - for (int j = 0; j < 3; j++) Y[3 * i + j] = g[j]; + for (int j = 0; j < 3; j++) Y[i + j] = g[j]; } // Integrate. ensureHavePoissonSolver(); - Vector div = -geom.polygonDivergenceMatrix * Y; + Vector div = geom.simplePolygonDivergenceMatrix * Y; Vector distances = poissonSolver->solve(div); // Shift solution. distances -= distances.minCoeff() * Vector::Ones(V); - geom.unrequirePolygonGradientMatrix(); - geom.unrequirePolygonDivergenceMatrix(); + geom.unrequireSimplePolygonGradientMatrix(); + geom.unrequireSimplePolygonDivergenceMatrix(); geom.unrequireVertexIndices(); return VertexData(mesh, distances); @@ -189,29 +190,47 @@ VertexData PolygonMeshHeatSolver::computeSignedDistance(const std::vecto ensureHaveVectorHeatSolver(); Vector> Xt = vectorHeatSolver->solve(X0); - // Average onto faces, and normalize. - size_t F = mesh.nFaces(); - Vector Y(3 * F); + // Express normalized diffused vectors extrinsically. + Eigen::MatrixXd Yt(V, 3); + for (size_t i = 0; i < V; i++) { + Vector3 g = geom.vertexTangentBasis[i][0] * std::real(Xt[i]) + geom.vertexTangentBasis[i][1] * std::imag(Xt[i]); + g / g.norm(); + for (int j = 0; j < 3; j++) { + Yt(i, j) = g[j]; + } + } + + // In order to apply the divergence operator as-is, average onto faces of refinement. geom.requireFaceIndices(); + geom.requireSimplePolygonDivergenceMatrix(); + geom.requireSimplePolygonProlongationMatrix(); + SparseMatrix D = geom.simplePolygonDivergenceMatrix; + SparseMatrix P = geom.simplePolygonProlongationMatrix; // (|V| + |F|) x |V| + Eigen::MatrixXd PY = P * Yt; // (|V| + |F|) x 3 + Vector Y(D.cols()); // 3|T^f| x 1 + int c = 0; 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(); + int i = 0; size_t fIdx = geom.faceIndices[f]; - for (int j = 0; j < 3; j++) { - Y[3 * fIdx + j] = Yf[j]; + Eigen::Vector3d valM = PY.row(V + fIdx); + for (Halfedge he : f.adjacentHalfedges()) { + size_t v0 = geom.vertexIndices[he.tailVertex()]; + size_t v1 = geom.vertexIndices[he.tipVertex()]; + Eigen::Vector3d val1 = PY.row(v0); + Eigen::Vector3d val2 = PY.row(v1); + Eigen::Vector3d avg = (valM + val1 + val2) / 3.; + avg /= avg.norm(); + for (int j = 0; j < 3; j++) { + int idx = c + 3 * i + j; + Y[idx] = avg[j]; + } + i++; } + c += f.degree() * 3; } - geom.unrequireFaceIndices(); - geom.unrequireVertexTangentBasis(); - geom.requirePolygonDivergenceMatrix(); - Vector divYt = geom.polygonDivergenceMatrix * Y; - geom.unrequirePolygonDivergenceMatrix(); + Vector divYt = D * Y; + geom.unrequireSimplePolygonDivergenceMatrix(); Vector phi; if (levelSetConstraint == LevelSetConstraint::None) { diff --git a/src/surface/polygon_mesh_helpers.cpp b/src/surface/polygon_mesh_helpers.cpp deleted file mode 100644 index 682ffb80..00000000 --- a/src/surface/polygon_mesh_helpers.cpp +++ /dev/null @@ -1,131 +0,0 @@ -#include "geometrycentral/surface/polygon_mesh_helpers.h" - -namespace geometrycentral { -namespace surface { - -// Helper functions for Bunge et al. "Polygon Laplacian Made Simple" (2020). -// 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) - -Eigen::VectorXd simplePolygonVirtualVertex(const Eigen::MatrixXd& poly) { - - // 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; -} - -// Helper functions for de Goes et al. "Discrete Differential Operators on Polygonal Meshes" (2020). -// 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) - -Eigen::MatrixXd polygonAveragingMatrix(const Face& f) { - 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 polygonDerivativeMatrix(const Face& f) { - 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::Matrix3d bracket(const Eigen::Vector3d& n) { - Eigen::Matrix3d B; - B << 0., -n[2], n[1], n[2], 0., -n[0], -n[1], n[0], 0.; - return B; -} - -Eigen::Vector3d project(const Eigen::Vector3d& u, const Eigen::Vector3d& n) { - // Project u on the orthgonal of n. - return u - (u.dot(n) / n.squaredNorm()) * n; -} - -Eigen::MatrixXd kroneckerWithI2(const Eigen::MatrixXd& M) { - 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 \ No newline at end of file diff --git a/test/src/polygon_operators_test.cpp b/test/src/polygon_operators_test.cpp index a9e9756c..4d7a892e 100644 --- a/test/src/polygon_operators_test.cpp +++ b/test/src/polygon_operators_test.cpp @@ -56,19 +56,15 @@ TEST_F(PolygonMeshSuite, TriangularTest) { 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(); @@ -76,41 +72,5 @@ TEST_F(PolygonMeshSuite, TriangularTest) { 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