Clean up polygon DEC

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<double> 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<double> 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<double> 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<double> 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<double> 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<double> 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<double> 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<std::complex<double>> 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<double> EmbeddedGeometryInterface::polygonHodge0` A $|V| \times |V|$ diagonal matrix
-    - `Eigen::SparseMatrix<double> EmbeddedGeometryInterface::polygonHodge0Inverse` A $|V| \times |V|$ diagonal matrix
-    - `Eigen::SparseMatrix<double> EmbeddedGeometryInterface::polygonHodge1` An $|H| \times |H|$ matrix, not necessarily diagonal
-    - `Eigen::SparseMatrix<double> EmbeddedGeometryInterface::polygonHodge2` An $|F| \times |F|$ diagonal matrix
-    - `Eigen::SparseMatrix<double> EmbeddedGeometryInterface::polygonHodge2Inverse` An $|F| \times |F|$ diagonal matrix
-    - `Eigen::SparseMatrix<double> EmbeddedGeometryInterface::polygonD0` An $|H| \times |V|$ matrix with $\{-1, 0, 1\}$ entries
-    - `Eigen::SparseMatrix<double> EmbeddedGeometryInterface::polygonD1` An $|F| \times |H|$ matrix with $\{-1, 0, 1\}$ entries
-
     Only valid on an `EmbeddedGeometryInterface`.
 
-    - **require:** `void EmbeddedGeometryInterface::requirePolygonDECOperators()`
+    - **member:** `Eigen::SparseMatrix<std::complex<double>> EmbeddedGeometryInterface::simplePolygonVertexConnectionLaplacian`
+    - **require:** `void EmbeddedGeometryInterface::requireSimplePolygonVertexConnectionLaplacian()`
 
 ## Extrinsic angles
 

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<double> simplePolygonLaplacian;
   void requireSimplePolygonLaplacian();
   void unrequireSimplePolygonLaplacian();
 
+  // Divergence
+  Eigen::SparseMatrix<double> simplePolygonDivergenceMatrix;
+  void requireSimplePolygonDivergenceMatrix();
+  void unrequireSimplePolygonDivergenceMatrix();
+
+  // Gradient
+  Eigen::SparseMatrix<double> simplePolygonGradientMatrix;
+  void requireSimplePolygonGradientMatrix();
+  void unrequireSimplePolygonGradientMatrix();
+
+  Eigen::SparseMatrix<double> simplePolygonProlongationMatrix;
+  void requireSimplePolygonProlongationMatrix();
+  void unrequireSimplePolygonProlongationMatrix();
+
+  // connection Laplacian
+  Eigen::SparseMatrix<std::complex<double>> simplePolygonVertexConnectionLaplacian;
+  void requireSimplePolygonVertexConnectionLaplacian();
+  void unrequireSimplePolygonVertexConnectionLaplacian();
+
   // Vertex Galerkin mass matrix (unlumped)
   Eigen::SparseMatrix<double> 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<double> polygonLaplacian;
-  void requirePolygonLaplacian();
-  void unrequirePolygonLaplacian();
-
-  // Gradient matrix
-  Eigen::SparseMatrix<double> polygonGradientMatrix;
-  void requirePolygonGradientMatrix();
-  void unrequirePolygonGradientMatrix();
-
-  // Divergence matrix
-  Eigen::SparseMatrix<double> polygonDivergenceMatrix;
-  void requirePolygonDivergenceMatrix();
-  void unrequirePolygonDivergenceMatrix();
-
-  // Vertex mass matrix (lumped)
-  Eigen::SparseMatrix<double> polygonVertexLumpedMassMatrix;
-  void requirePolygonVertexLumpedMassMatrix();
-  void unrequirePolygonVertexLumpedMassMatrix();
-
-  // Vertex connection Laplacian
-  Eigen::SparseMatrix<std::complex<double>> polygonVertexConnectionLaplacian;
-  void requirePolygonVertexConnectionLaplacian();
-  void unrequirePolygonVertexConnectionLaplacian();
-
-  Eigen::SparseMatrix<double> 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<Eigen::SparseMatrix<double>> simplePolygonLaplacianQ;
   virtual void computeSimplePolygonLaplacian();
 
+  // Divergence
+  DependentQuantityD<Eigen::SparseMatrix<double>> simplePolygonDivergenceMatrixQ;
+  virtual void computeSimplePolygonDivergenceMatrix();
+
+  // Gradient
+  DependentQuantityD<Eigen::SparseMatrix<double>> simplePolygonGradientMatrixQ;
+  virtual void computeSimplePolygonGradientMatrix();
+
+  // Prolongation
+  DependentQuantityD<Eigen::SparseMatrix<double>> simplePolygonProlongationMatrixQ;
+  virtual void computeSimplePolygonProlongationMatrix();
+
+  // Connection Laplacian
+  DependentQuantityD<Eigen::SparseMatrix<std::complex<double>>> simplePolygonVertexConnectionLaplacianQ;
+  virtual void computeSimplePolygonVertexConnectionLaplacian();
+
   // Vertex mass matrix (unlumped)
   DependentQuantityD<Eigen::SparseMatrix<double>> simplePolygonVertexGalerkinMassMatrixQ;
   virtual void computeSimplePolygonVertexGalerkinMassMatrix();
@@ -156,70 +161,22 @@ class EmbeddedGeometryInterface : public ExtrinsicGeometryInterface {
   DependentQuantityD<Eigen::SparseMatrix<double>> 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<Eigen::VectorXd> virtualRefinementAreaWeights;
+  FaceData<Eigen::Vector3d> virtualRefinementAreaPoints;
   DependentQuantityD<FaceData<Eigen::VectorXd>> virtualRefinementAreaWeightsQ; // affine weights for each virtual node
+  DependentQuantityD<FaceData<Eigen::Vector3d>> virtualRefinementAreaPointsQ;
   virtual void computeVirtualRefinementAreaWeights();
   virtual Eigen::MatrixXd simplePolygonMassMatrix(const Face& f);
   virtual Eigen::MatrixXd simplePolygonStiffnessMatrix(const Face& f);
-  virtual SparseMatrix<double> simplePolygonProlongationMatrix();
+  virtual SparseMatrix<double> 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<Eigen::SparseMatrix<double>> polygonLaplacianQ;
-  virtual void computePolygonLaplacian();
-
-  DependentQuantityD<Eigen::SparseMatrix<double>> polygonGradientMatrixQ;
-  virtual void computePolygonGradientMatrix();
-
-  DependentQuantityD<Eigen::SparseMatrix<double>> polygonDivergenceMatrixQ;
-  virtual void computePolygonDivergenceMatrix();
-
-  // Vertex mass matrix (lumped)
-  DependentQuantityD<Eigen::SparseMatrix<double>> polygonVertexLumpedMassMatrixQ;
-  virtual void computePolygonVertexLumpedMassMatrix();
-
-  // Vertex connection Laplacian
-  DependentQuantityD<Eigen::SparseMatrix<std::complex<double>>> polygonVertexConnectionLaplacianQ;
-  virtual void computePolygonVertexConnectionLaplacian();
-
-  // DEC Operators
-  std::array<Eigen::SparseMatrix<double>*, 7> polygonDECOperatorArray;
-  DependentQuantityD<std::array<Eigen::SparseMatrix<double>*, 7>> polygonDECOperatorsQ;
-  virtual void computePolygonDECOperators();
-
-  const double polygonLambda = 1.0;
-  VertexData<Eigen::VectorXd> polygonVertexNormals;
-  DependentQuantityD<VertexData<Eigen::VectorXd>> 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);
 };