A Hands-On Introduction to Discrete Differential Operators on Polygon Meshes

Other Laplacians on Polygon Meshes

Photo of Sven D. Wagner

Sven D. Wagner

TU Dortmund

Photo of Astrid Bunge

Astrid Bunge

AutoForm Engineering

Photo of Mario Botsch

Mario Botsch

TU Dortmund

🚀 by Decker

Outline

  • Other Polygon Operators
    • Martin et al., Polyhedral finite elements using harmonic basis functions, SGP 2008
    • Alexa & Wardetzky, Discrete Laplacians on general polygonal meshes, SIG 2011
    • de Goes et al., Discrete differential operators on polygonal meshes, SIG 2020
    • Bunge et al., The Diamond Laplace for polygonal and polyhedral meshes, SGP 2021
  • Comparisons & Conclusion

Barycentric Coordinates

  • Piecewise linear functions
  • Interpolate nodal data: \(\varphi_i\of{\vec{x}_j} = \delta_{ij}\)
  • Partition of unity: \(\sum_i \varphi_i\of{\vec{x}} = 1\)
  • Barycentric property: \(\sum_i \varphi_i\of{\vec{x}} \vec{x}_i = \vec{x}\)
  • \(C^0\) across elements, \(C^1\) within elements

Generalized Barycentric Coordinates

  • Piecewise linear functions
  • Interpolate nodal data: \(\varphi_i\of{\vec{x}_j} = \delta_{ij}\)
  • Partition of unity: \(\sum_i \varphi_i\of{\vec{x}} = 1\)
  • Barycentric property: \(\sum_i \varphi_i\of{\vec{x}} \vec{x}_i = \vec{x}\)
  • \(C^0\) across elements, \(C^1\) within elements

Generalized Barycentric Coordinates

  • Wachspress Coordinates
    • Wachspress, Warren
  • Mean Value Coordinates
    • Floater, Hormann, Ju, Wicke
  • Maximum Entropy Coordinates
    • Sukumar, Hormann
  • Harmonic Coordinates
    • Joshi, Martin

images/harmonic-coordinates-1.png

Hormann and Sukumar 2017; Martin et al. 2008

Laplace Matrix & Mass Matrix

\[\begin{eqnarray*} \mat{L}_{ij} &=& \int \grad \varphi_i \cdot \grad \varphi_j &=& \color{red}{\xcancel{ \begin{cases} -\frac{\cot\alpha_{ij}+\cot\beta_{ij}}{2} & \text{if } j\in\set{N}\of{i} \,, \\[0.5em] \displaystyle -\sum_{k\in\set{N}\of{i}} \mat{L}_{ik} & \text{if } j=i \,, \\[0.3em] 0 & \text{otherwise}. \end{cases} }} \\[1em] \mat{M}_{ij} &=& \int \varphi_i \, \varphi_j &=& \color{red}{\xcancel{ \begin{cases} \frac{\abs{t_{ijk}} + \abs{t_{jih}}}{12} & \text{if } j\in\set{N}\of{i}\,, \\[0.5em] \displaystyle \sum_{k\in\set{N}\of{i}} \mat{M}_{ik} & \text{if }j=i \,,\\[0.3em] 0 & \text{otherwise}. \end{cases} }} \end{eqnarray*}\]

On triangles, harmonic coordinates reproduce barycentric coordinates

Demo: Parameterization

Properties

Poisson System on 2D Voronoi Meshes

13.5053,	28.5194,	58.8386,	119.78,	241.153
Standard LVR, 0.00768401,0.00193955,0.00057802,0.000152567,3.41754E-05
Harmonic,	0.00951959,	0.00236579,	0.000699382,	0.000188646,	5.46139E-05
Robust LVR ,	0.00795767, 0.00198379 , 0.000572063, 0.000147802,3.33787E-05


<!--
{
  "options": {
    "scales": {
      "x": {
        "title": {
          "text": "inverse mean edge length",
          "display": true
        },
        "type": "logarithmic"
      },
      "y": { 
        "title": {
          "text": "L2 error",
          "display": true
        },
        "type": "logarithmic"
      }
    }
  }
}
-->

Expected convergence rate 😄 Expensive to compute, evaluate, and integrate 😢

images/VoronoiPlane.svg

Harmonic Finite Elements

Outline

  • Other Polygon Operators
    • Martin et al., Polyhedral finite elements using harmonic basis functions, SGP 2008
    • Alexa & Wardetzky, Discrete Laplacians on general polygonal meshes, SIG 2011
    • de Goes et al., Discrete differential operators on polygonal meshes, SIG 2020
    • Bunge et al., The Diamond Laplace for polygonal and polyhedral meshes, SGP 2021
  • Comparisons & Conclusion

DEC Polygon Laplacians

images/alexa-laplace.png
  • Discretization of \(\laplace f = \delta d f\) (Laplace-de Rham)
    • The exterior derivative \(d\) becomes the difference operator \(\mat{d}\)
    • The co-differential \(\delta\) can be discretized as \(\mat{M}_0^{-1} \mat{d}\T \mat{M}_1\)
    • Putting it together: \(\laplace \vec{u} \;\approx\; -\mat{M}_0^{-1} \mat{d}\T \mat{M}_1 \mat{d}\, \vec{u}\)

Difference operator

  • The difference operator maps

Difference operator

  • The difference operator maps
    • From values on vertices (0-forms)
images/0-form.svg
0-simplex

Difference operator

  • The difference operator maps
    • From values on vertices (0-forms)
    • To differece values on (half-)edges (1-forms)
images/1-form.svg
1-simplex

Difference operator

  • Construction for face \(f\) with \(n_f\) vertices/edges \[\mat{d}_f = \begin{pmatrix} 1 & 0 & 0 & \cdots & 0 & -1\\ -1 & 1 & 0 & \cdots & 0 & 0\\ 0 & -1 & 1 & \cdots & 0 & 0\\ & & & \vdots \\ 0 & 0 & 0 & \dots & -1 & 1\end{pmatrix}\in \R^{n_f \times n_f}\]
  • Note that \(\mat{d}\) is rectangular, mapping from all vertices to all halfedges

Properties of DEC Laplacian

\(\laplace \vec{u} \;\approx\; -\mat{M}_0^{-1} \underbrace{\mat{d}\T \mat{M}_1 \mat{d}}_{-\mat{L}} \, \vec{u}\)

  • The weak Laplacian \(\mat{L}\) appears as \(-\mat{d}\T \mat{M}_1 \mat{d}\)
    • Constants are in kernel of \(\mat{L}\) because \(\mat{d}\) annihilates them 👍
    • If \(\mat{M}_1\) is symmetric then \(\mat{L}\) is symmetric 👍
    • If \(\mat{M}_1\) is regular then only the constants are in the kernel 👍
  • Linear precision of \(\mat{L}\)?
    • Strategy: construction per face as \(\mat{M}_f\)
    • Require that \(\mat{L}_f = \mat{d}_f\T \mat{M}_f \mat{d}_f\) yields the area gradient as \(\left(\mat{L}_f\mat{X}_f\right)_i\) for vertex \(i\) 👍
  • Problem: linear precision + \(\mat{M}_f\) regular
    • Solved by regularization that introduces a parameter that affects the results 👎

Regularizer

  • Add regularizer matrix \(\mat{R}_f\) to fill up kernel \[\mat{M}_{f} = \mat{\tilde{M}}_f + \lambda\mat{R}_f\]
    • Both methods are inspired by Brezzi et al. for inner product matrix
    • Use different motivations for regularizer
    • Hyperparameter \(\lambda\) necessary for both discretizations

Brezzi et al. 2005

Laplace and Mass Matrices

  • Local Laplace matrix per polygon \[ {\mat{L}}_f = \mat{d}\T \mat{M}_{f} \mat{d} \]
    • Assemble local matrices into global Laplace matrix \(\mat{L}\)
  • Global mass matrix \(\mat{M}\) with
    \[\mat{M}_{ii} = \sum_{f \ni v_i} \frac{|\bar{f}|}{n_f}\]

Properties

Poisson System on 2D Voronoi Meshes

13.5053,	28.5194,	58.8386,	119.78,	241.153
Alexa11; λ=2,	0.00838748,	0.00209687,	0.000476346,	0.000127937,	3.30444E-05
Alexa11; λ=1,	0.00456036,	0.0010483,	0.000382806,	8.53325E-05,	1.84683E-05
Alexa11; λ=0.5,	0.00800984,	0.00193734,	0.000621653,	0.000143936,	3.19937E-05
Alexa11; λ=0.1,	0.0144414,	0.00337097,	0.000947331,	0.000221091,	5.14271E-05
deGoes20; λ=2,	  0.00965471,	0.00245388,	  0.000542276,	0.000145297,	3.74297E-05
deGoes20; λ=1,	  0.00435302,	0.000996788,	0.000354402,	7.95537E-05,	1.73643E-05
deGoes20; λ=0.5,	0.00747841,	0.00182759,	  0.00059557,	  0.000138136,	3.05869E-05
deGoes20; λ=0.1,	0.0139834,	0.00327528,	  0.00092799,	  0.000217473,	5.03785E-05
<!--
{
  "options": {
    "scales": {
      "x": {
        "title": {
          "text": "inverse mean edge length",
          "display": true
        },
        "type": "logarithmic"
      },
      "y": { 
        "title": {
          "text": "L2 error",
          "display": true
        },
        "type": "logarithmic"
      }
    }
  }
}
-->

images/VoronoiPlane.svg

Poisson System on 2D Voronoi Meshes

13.5053,	28.5194,	58.8386,	119.78,	241.153
Standard LVR, 0.00768401,0.00193955,0.00057802,0.000152567,3.41754E-05
Harmonic,	0.00951959,	0.00236579,	0.000699382,	0.000188646,	5.46139E-05
Alexa11; λ=1,	0.00456036,	0.0010483,	0.000382806,	8.53325E-05,	1.84683E-05,
Alexa11; λ=0.1,	0.0144414,	0.00337097,	0.000947331,	0.000221091,	5.14271E-05
deGoes20; λ=1,	  0.00435302,	0.000996788,	0.000354402,	7.95537E-05,	1.73643E-05
deGoes20; λ=0.1,	0.0139834,	0.00327528,	  0.00092799,	  0.000217473,	5.03785E-05

<!--
{
  "options": {
    "scales": {
      "x": {
        "title": {
          "text": "inverse mean edge length",
          "display": true
        },
        "type": "logarithmic"
      },
      "y": { 
        "title": {
          "text": "L2 error",
          "display": true
        },
        "type": "logarithmic"
      }
    }
  }
}
-->

images/VoronoiPlane.svg

Demo: Curvature

Outline

  • Other Polygon Operators
    • Martin et al., Polyhedral finite elements using harmonic basis functions, SGP 2008
    • Alexa & Wardetzky, Discrete Laplacians on general polygonal meshes, SIG 2011
    • de Goes et al., Discrete differential operators on polygonal meshes, SIG 2020
    • Bunge et al., The Diamond Laplace for polygonal and polyhedral meshes, SGP 2021
  • Comparisons & Conclusion

Finite Volumes

Main idea: Consider the integral of differential in a small region

\[ \text{Divergence: } \quad \iint_{\Omega} \text{div} \vec{u}\ \d{A} = \oint_{\partial \Omega} \vec{u}\tp\vec{n} \d{s}\] \[ \text{Gradient: } \quad \iint_\Omega \nabla u \d{A} = \oint_{\partial \Omega} u\, \mat{n} \d{s}\] \[ \text{Laplacian: } \quad \iint_\Omega \Delta u \d{A} = \oint_{\partial \Omega} \nabla u\tp \mat{n} \d{s}\]

Discrete Duality Finite Volumes (DDFV)

images/mesh.svg images/dual_mesh.svg images/dual_controlvolume.svg images/primal_controlvolume.svg images/dual_mesh_diamond.svg images/diamond_mesh.svg

Hermeline 2000

Quiz

How do we get diamond cells on an arbitrary polygon mesh?

images/poly_mesh_dual.svg images/poly_mesh_diamond.svg

2D DDFV

images/our_diamond.png

\[ \grad u|_D \;=\; \frac{1}{2 \abs{D}} \sum_{(i,j) \in \partial D} \vec{e}_{ij}^\perp \frac{u_i+u_j}{2} \]

Virtual Dual Vertices

images/poly11.svg \[\downarrow\] images/poly31.svg \[\downarrow\] images/poly41.svg

  • Combine ideas from DDFV and virtual refinement
  • Use virtual vertices as dual mesh
  • Compute DDFV operator on the “virtual diamonds”
    • Express diamond in an intrinsic 2D coordinate system
  • Restrict values back to the original mesh

Diamond Gradient Operator

  • The \(i\)-th column of the diamonds gradient matrix \(\mat{G}_D \in \R^{2 \times 4}\) \[ \mat{G}_D(:,i) = \frac{1}{2 \abs{D}} \sum_{(i,j) \in \partial D} \vec{e}_{ij}^\perp \]
  • Mesh gradient \(\hat{\mat{G}}\) is combined with prolongation \(\mat{P}\) to form global operator \[ \mat{G} = \hat{\mat{G}} \, \mat{P} \,\in\R^{2\abs{\mathcal{E}} \times \abs{\mathcal{V}}} \]

Diamond Divergence Operator

  • The diamond divergence matrix is defined as \[ \mat{D} = -\mat{P}\tp \hat{\mat{G}}\tp \hat{\mat{M}}_\diamond \]
  • \(\hat{\mat{M}}_\diamond\) is diagonal matrix of diamond areas

Laplace Operator

  • Laplacian matrix is defined as \[ \mat{L} \;=\; \underbrace{-\mat{P}\tp \hat{\mat{G}}\tp \hat{\mat{M}}_\diamond}_{\mat{D}} \cdot \underbrace{\hat{\mat{G}} \mat{P}}_{\mat{G}} \]
  • Mass matrix derived from standard DDFV mass matrix \(\mat{M}_\diamond\) \[ \mat{M} \;=\; \mat{P}\tp \mat{M}_\diamond \mat{P} \]

Demo: Geodesics

Properties

Poisson System on 2D Voronoi Meshes

13.5053,	28.5194,	58.8386,	119.78,	241.153
Standard LVR, 0.00768401,0.00193955,0.00057802,0.000152567,3.41754E-05
Harmonic,	0.00951959,	0.00236579,	0.000699382,	0.000188646,	5.46139E-05
Alexa11; λ=1,	0.00456036,	0.0010483,	0.000382806,	8.53325E-05,	1.84683E-05
deGoes20; λ=1,	  0.00435302,	0.000996788,	0.000354402,	7.95537E-05,	1.73643E-05
Diamond, 0.00610683,0.00155889,0.000413685,0.000105977,2.46064E-05
<!--
{
  "options": {
    "scales": {
      "x": {
        "title": {
          "text": "inverse mean edge length",
          "display": true
        },
        "type": "logarithmic"
      },
      "y": { 
        "title": {
          "text": "L2 error",
          "display": true
        },
        "type": "logarithmic"
      }
    }
  }
}
-->

images/VoronoiPlane.svg

Outline

  • Other Polygon Operators
    • Martin et al., Polyhedral finite elements using harmonic basis functions, SGP 2008
    • Alexa & Wardetzky, Discrete Laplacians on general polygonal meshes, SIG 2011
    • de Goes et al., Discrete differential operators on polygonal meshes, SIG 2020
    • Bunge et al., The Diamond Laplace for polygonal and polyhedral meshes, SGP 2021
  • Comparisons & Conclusion

Quantitative Comparisons

  • Ensure fair comparisons
  • All methods (re)implemented in C++
  • Same meshes & test for all methods
  • Code checked by original authors (Thanks 👍)
  • Source code for all Laplace operators, quantitative tests, and interactive demos is available on github

Poisson System on 2D Meshes

13.5053,	18.8048, 24.0398 , 37.072 ,50.2679
Robust LVR ,	0.00745003, 0.0036901, 0.00222217 , 0.000924393 , 0.000501055
Standard LVR,	0.0385146, 0.0226358, 0.0147572, 0.00647153, 0.00376134
Alexa11; λ=2, 0.0611674, 0.041198, 0.0291081, 0.0144778, 0.00876606
Alexa11; λ=0.5, 0.01545 , 0.00859024 , 0.00546532 , 0.0023337,0.00132841
deGoes20; λ=1, 0.0363144 , 0.0221452, 0.0148121, 0.00684791, 0.00400674
deGoes20; λ=0.1, 0.0252731, 0.0122202, 0.00700861, 0.00268186, 0.00139388
Diamond, 0.0131201, 0.00576263 , 0.0029424, 0.00113579 , 0.00071783
<!--
{
  "options": {
    "scales": {
      "x": {
        "title": {
          "text": "inverse mean edge length",
          "display": true
        },
        "type": "logarithmic"
      },
      "y": { 
        "title": {
          "text": "L2 error",
          "display": true
        },
        "type": "logarithmic"
      }
    }
  }
}
-->

images/HexStretchedPlane.svg

Poisson System on Sphere Meshes

2.76125, 6.02378, 12.5599, 25.6378, 51.7971
Robust LVR , 0.0195403, 0.00757413, 0.00353687, 0.00173384, 0.000859737
Standard LVR,	 0.0292485, 0.0119197, 0.00558098, 0.00270908, 0.00133463
Alexa11; λ=2, 0.0893085, 0.0490686, 0.0258716, 0.0131541, 0.00660642
Alexa11; λ=0.5, 0.0138582, 0.00877119, 0.00489306, 0.00251059, 0.00126177
deGoes20; λ=1, 0.0500151, 0.0278916, 0.0147208, 0.00751131, 0.00378149
deGoes20; λ=0.1, 0.0267659, 0.0115076 , 0.00583512 , 0.00289191, 0.00143387
Diamond, 0.024761, 0.00765354, 0.00354862, 0.0017351, 0.000857455
<!--
{
  "options": {
    "scales": {
      "x": {
        "title": {
          "text": "inverse mean edge length",
          "display": true
        },
        "type": "logarithmic"
      },
      "y": { 
        "title": {
          "text": "L2 error",
          "display": true
        },
        "type": "logarithmic"
      }
    }
  }
}
-->

images/ConcaveSphere.png

Computational Performance

08_polygon_comparison-deck-code-12300b3a.tex.svg
Time for building and solving a Poisson system

Conclusion

Recommendations

  • DEC Approaches
  • Diamond Laplacian
  • Virtual Refinement Method
  • Harmonic shape functions

💾   Source code of all polygon methods available on github

Conclusion

  • Gave a hands on introduction to Cotan Laplacian
    • Discussed its derivation and several applications
    • Showed methods to avoid numerical pitfalls on bad triangles
  • Presented recent progress for Polygon Laplacians
    • Discuss individual discretization strategies
    • Analyze strength and weaknesses

images/qr-code.svg

  • Acknowledgments
    • Collaborators: Marc Alexa, Philipp Herholz, Misha Kazhdan, Olga Sorkine-Hornung
    • Code-Checkers: Fernando de Goes, Max Wardetzky

References

Alexa, Marc, and Max Wardetzky. 2011. “Discrete Laplacians on General Polygonal Meshes.” ACM Transactions on Graphics 30 (4).
Brezzi, Franco, Konstantin Lipnikov, and Valeria Simoncini. 2005. “A Family of Mimetic Finite Difference Methods on Polygonal and Polyhedral Meshes.” Mathematical Models and Methods in Applied Sciences 15 (10).
deGoes, Fernando, Andrew Butts, and Mathieu Desbrun. 2020. “Discrete Differential Operators on Polygonal Meshes.” ACM Transactions on Graphics 39 (4).
Hermeline, F. 2000. “A Finite Volume Method for the Approximation of Diffusion Operators on Distorted Meshes.” J. Comput. Phys. 160 (2).
Hormann, Kai, and N. Sukumar. 2017. Generalized Barycentric Coordinates in Computer Graphics and Computational Mechanics. Taylor & Francis.
Martin, Sebastian, Peter Kaufmann, Mario Botsch, Martin Wicke, and Markus Gross. 2008. “Polyhedral Finite Elements Using Harmonic Basis Functions.” Computer Graphics Forum 27 (5).