Thanks for making this project available!

The latest version on github will not compile using Microsoft Visual Studio Community 2019 version 16.6.3.
The error is on line 108 of element.h
// Friends
friend std::ostream& std::operator<<<>(std::ostream& output, const Element<T, M>& e);

And the error message is:
Error C2063 'std::operator <<': not a function

Because of the way it's named, Vector::normalize() yields unexpected behavior: the verb-forming suffix "-ize" implies a transitive verb that acts on the vector itself, as in, "We then normalize the vector." The expected behavior of v.normalize() is hence that v itself gets modified, whereas the current behavior is to return a normalized version of v. [The behavior is spelled out in the documentation, but causes confusion nonetheless.]

Proposal: Add a method Vector::unit() that returns the unit version of the given vector. Deprecate Vector::normalize() (by removing it from the documentation) to preserve compatibility with existing code.

Though in a perfect world normalize() would mutate the vector itself, making that change now would likely wreak havok.

Would there be any interest in an implementation of this algorithm within geometry central?
http://mesh.brown.edu/dgp/pdfs/Jones-sg03.pdf
I have a naive implemetation that I would be happy to try to clean up and submit as a PR.
But, maybe a better question is this: Is there a plan already in the works for a suite of smoothing algorithms within the library? If so, what would that look like? Or is it already there somewhere and I've just missed it?

This may fall under the heading of "too much sugar", but it could be nice to have accessors for the two endpoints of an edge. Currently there are neither accessors like tip() and tail(), nor is there an iterator Edge::adjacentVertices(). The latter is of course useful when you need to do the same thing at both endpoints. The former is natural if you want to, say, average values at the two endpoints, or (more exotic example) perform a conformal scaling. Of course tip and tail make sense only for an oriented edge, but perhaps Edge::firstVertex() and Edge::secondVertex() (or similar) could be used.

A fairly common stencil for edges in a (manifold) triangle mesh is the boundary of the "diamond" made by the two incident triangles. Maybe this is getting too specialized. But it might not be so bad to have a method Edge::diamondBoundary() that iterates over e.halfedge().next(), e.halfedge.next().next(), e.halfedge.twin().next(), and e.halfedge().twin().next().next(). This list is conceptually simple, but currently translates to some pretty ugly code when written inline.

One simple implementation, perhaps, is to just construct and return this list (rather than defining a new iterator class):

array<Edge,4> diamondBoundary( Edge e )
{
   Halfedge h = e.halfedge();
   Halfedge hN = h.next();
   Halfedge hT = h.twin();
   Halfedge hTN = hT.next();
   return { hN.edge(), hN.next().edge(), hTN.edge(), hTN.next().edge() };
}

Perhaps, again, this is getting too specialized and it is better to just let someone write this code themselves. But adding additional "convenience" neighborhoods doesn't seem to add too much glut to the code/documentation, and will only make other code leaner, cleaner, and less error-prone.

Hi I might just not have understood things correctly but I was wondering if there was a way or a plan to implement a method to convert the Vector2 vectors from the Vector field algorithms to real world Vector3 vectors without using polyscope? Or if there is a way to get the real world Vector3 from polyscope ?

Dear nmwsharp,
I am currently trying to use the heat method distance algorithm in GC for some .off datasets. Here is the link of datasets I used:
https://github.com/ItakEjgo/TerrainDataset
However, I met some problems. I follow your gc tutorial for the heat method distance algorithm through the following code:

HeatMethodDistanceSolver heatSolver(*geometry);
for (decltype(mesh->nVertices()) i = 0; i < mesh->nVertices(); i++){
    Vertex v = mesh->vertex(i);
    VertexData<double> distToSource = heatSolver.computeDistance(v);
}

This algorithm wells really well for some datasets and performs really well, e.g. small_terrain.off and BH_low.off.
But, it failed in the other .off datasets with the following exception:

terminate called after throwing an instance of 'std::logic_error'
  what():  checkFinite() failure. Non-finite row vector entry [39] = -nan
Aborted (core dumped)

The value of the vector became -nan and thus the checkFinite() function throws the previous exception for these datasets. I tried to fix this problem and found that the problem is happened in the following line in heat_method_distance.cpp line 174:
Vector<double> distVec = poissonSolver->solve(divergenceVec);
But I have no idea about solving this problem...
All these dataset could be successfully visulized through the mesh lab in https://github.com/cnr-isti-vclab/meshlab.
Could you kindly give me some adivce about how to solve this problem？Are the datasets I used are vaild or invalid?
Thanks in advance for your help!

In the function SimplePolygonMesh::readMeshFromBinaryStlFile(std::istream& in) we need to seek back to the beginning of the file before processing the 80 byte header and triangle count. In the calling function, the stream is examined to see if the first line contains the token "solid", which flags the file as an ASCII STL file. However, once we determine that the file represents a binary STL, the current stream position is at the end of the first "delimited line" in the file. So, inserting a in.seekg(0, in.beg); before the header read fixes the issue.

Hi, Nick!

Thank you for sharing your library!
[Question]
Could you please explain how can I initialize SurfaceMesh with Eigen::MatrixXd (vertices and faces); and further use SurfaceMesh for geodesic distance calculation. I've looked through docs but didn't find it.

Currently the default behavior of void IntrinsicGeometryInterface::requireLaplacian() is to build the zero-Neumann Laplace matrix. Although there are of course many boundary conditions that could arise, two are common enough in practice that they may be worth supporting:

• Laplace with (possibly non-zero) Dirichlet boundary conditions
• Laplace with (possibly non-zero) Neumann boundary conditions

The exact interface for requesting these operators may require some thought. If all you care about is the matrix itself, then it might be as easy as adding a new method/member requireDirichletLaplacian()/dirichletLaplacian, and letting the current method requireLaplacian()/laplacian handle the Neumann case. Or the latter could be renamed to requireNeumannLaplace()/neumannLaplacian for clarity.

A more helpful, but also more complicated strategy would be to also build the vector of boundary data. For instance, to solve the Dirichlet problem

Δu = f on M,
u = g on 𝝏M

one needs to "move the boundary data to the right-hand side," i.e., any time a boundary degree of freedom appears in the expression for the discrete Laplacian, it gets multiplied by the cotan weight, negated, and moved to a vector b that becomes the right-hand side of the final linear system. A GC user can of course build this vector themselves, independent of the construction of the Laplacian, but it then feels like no real convenience is provided by the GC routine (since the user has to write some code of comparable complexity just to get the RHS vector).

Another possible interface could therefore look something like

DenseMatrix<double> IntrinsicGeometryInterface::requireLaplacian( BCType boundaryConditions = BCNeumann, MeshData<double>* boundaryData )

The default behavior would be to just build the zero-Neumann Laplace matrix. If BCDirichlet is passed in, it would build the (smaller) zero-Dirichlet matrix. If in addition the boundaryData pointer is non-null, it would return the vector that should appear on the right-hand side (which in the case of Neumann BC's is, IIRC, just a negated copy of the input vector).

There may be an even easier way of doing all this: the Dirichlet Laplacian is a submatrix of the Neumann Laplacian, corresponding to the interior vertices. So, the usage pattern could be to call requiresLaplacian() which builds the full matrix, and then call getDirichletLaplacian() or getNeumannLaplacian() which return the appropriate submatrices. These latter two methods could again take a vector of boundary data, and apply the right transformation to turn it into data on the RHS.

Finally, it's worth thinking about how to handle multiple RHSs. You don't want to re-build (or even extract) the same matrix many times just because you solve a sequence of many problems.

It's easy to introduce the following bug:

VertexData<double> x;
x[v] = 1.23;

Two possible resolutions:

• Don't allow a data array to be declared without a reference to a mesh (i.e., don't provide an empty constructor).
• Include assertions in debug mode that check for a mesh reference.

Pairwise conversions between:

• 1 forms on edges
• vector fields in face basis
• vector fields in vertex basis

including symmetric fields.

Note that most of the conversions are over determined or undetermined, so some kind of solve/projection will be required.

Sometimes when I require quantities inside of functions geometry-central doesn't fill the associated buffer. For example, if I try to run the following code, I segfault at the std::cout line.

#include "geometrycentral/surface/halfedge_mesh.h"
#include "geometrycentral/surface/meshio.h"
#include "geometrycentral/surface/vertex_position_geometry.h"

#include "args/args.hxx"
#include "imgui.h"

using namespace geometrycentral;
using namespace geometrycentral::surface;

// == Geometry-central data
std::unique_ptr<HalfedgeMesh> mesh;
std::unique_ptr<VertexPositionGeometry> geometry;

void badFn(VertexPositionGeometry geo) {
  geo.requireVertexAngleSums();
  std::cout << geo.vertexAngleSums[0] << std::endl;
}

int main(int argc, char **argv) {

  // Configure the argument parser
  args::ArgumentParser parser("geometry-central & Polyscope example project");
  args::Positional<std::string> inputFilename(parser, "mesh", "A mesh file.");

  // Parse args
  try {
    parser.ParseCLI(argc, argv);
  } catch (args::Help) {
    std::cout << parser;
    return 0;
  } catch (args::ParseError e) {
    std::cerr << e.what() << std::endl;
    std::cerr << parser;
    return 1;
  }

  // Make sure a mesh name was given
  if (!inputFilename) {
    std::cerr << "Please specify a mesh file as argument" << std::endl;
    return EXIT_FAILURE;
  }

  // Load mesh
  std::tie(mesh, geometry) = loadMesh(args::get(inputFilename));

  // Uncomment to fix
  // geometry->requireVertexAngleSums();
  badFn(*geometry);

  return EXIT_SUCCESS;
}

Error message:

Assertion failed: (i < size() && "Attempted to access MeshData with out of bounds index"), function operator[], file /Users/mark/Documents/Code/bad-gc-mesh/deps/geometry-central/src/../include/geometrycentral/surface/halfedge_containers.ipp, line 209.

If I call requireVertexAngleSums before calling badFn, then everything works fine.

You can find this example in a repo here: https://github.com/MarkGillespie/bad-gc-mesh.

It's currently possible to iterate over a MeshData instance x like this

for( size_t i = 0; i < x.size(); i++ ) {
   // do science
}

or like this

for( auto i : x.getMesh()->vertices() ) {
   // do science
}

but it would be nice to be able to just do this:

for( auto i : x ) {
   // do science
}

Hi,

It is not clear from your documentation if Suitesparse is actually a requirement, or just a nice to have. That is, if Suitesparse is not available on the system, will Eigen be used instead?

At this stage I would like to check out your geodesic distance features that use the different heat methods, but I will need to run it on Windows. Since Suitesparse is quite difficult to install on Windows, I was hoping to avoid that if not absolutely necessary.

Thanks for your time,
Michael

Container types such as VertexData<T> f( M ) effectively describe a map f: M -> T. If the type T supports arithmetic operations, it would therefore be natural to support the same arithmetic on the container type. E.g.,

VertexData<double> f(M), g(M), h(M);
double c;
h = f + c*g;

would be equivalent to

VertexData<T> f(M), g(M), h(M);
double c;
for( Vertex i : M.vertices() )
{
   h[i] = f[i] + c*g[i];
}

This syntax would greatly improve readability of code for, e.g., linear algebra and optimization, and bring geometry-central in better parity with packages that provide array-style manipulation.

The Eigen vector interoperability does not provide an attractive solution because:

• it requires the user to be familiar with Eigen,
• it still results in messy code, because of all the extra declarations and conversions,
• extra time is spent making copies, and
• most importantly, it can cause data to go out of synch, since methods like MeshData::toVector() make copies by value rather than by reference

This last point is particularly annoying for code that wants to intermingle container-wise and element-wise operations. E.g., a very common use case is: I want to treat a whole container like a vector for doing some kind of gradient descent or line search, but I want to use the nice Geometry methods (such as computing the area of a triangle, say) for expressing the objective function. Currently the only feasible usage pattern is:

• copy data from containers into Eigen vectors
• perform vector-wise operations on Eigen vectors
• copy Eigen values back into containers
• evaluate objective function using idiosyncratic geometry-central code
• repeat

or

• copy data from containers into Eigen vectors
• perform vector-wise operations on Eigen vectors
• copy Eigen values back into containers
• evaluate objective function by writing custom code that operates directly on Eigen vectors
• repeat

(But at this point, why bother using geometry-central? Just do everything in Eigen...)

Right now if you try to trace intrinsic edges on a mesh with completely degenerate (e.g., zero-area) triangles, it can cause a hard exit. Would be nice (e.g., for visualization purposes) if there's a "forgiving" version of this routine that just returns, say, an empty list for edges that are mathematically impossible to trace. For instance

signpostTri::traceEdges( bool beForgiving = false );

Currently one has to grab vertex positions from a VertexPositionGeometry object via

Vector3 p = geometry->InputVertexPositions[p];

Would be nice to provide access via a method VertexPositionGeometry::operator[], so that one can just write

Vector3 p = geometry[v];

One could also provide an analogous accessor for EdgeLengthGeometry, i.e.,

Edge ij;
double Lij = geometry[ij];

It looks like ManifoldSurfaceMesh::splitEdgeTriangular causes something strange to happen to the inputVertexPositions buffer. Here's the simplest example on which I see it:

void myCallback()
{
  if (ImGui::Button("Split one edge"))
  {
    std::cout << "Before: " << geometry->inputVertexPositions.size() << " positions in inputVertexPositions" << std::endl;
    std::cout << "Before: " << mesh->nVertices() << " vertices according to nVertices()" << std::endl;

    // Split an edge and set the new vertex to be at the midpoint
    Edge e = mesh->edge(1);
    Vector3 p1 = geometry->inputVertexPositions[e.halfedge().vertex()];
    Vector3 p2 = geometry->inputVertexPositions[e.halfedge().twin().vertex()];
    Vector3 newPos = (p1 + p2) / 2;
    Halfedge he = mesh->splitEdgeTriangular(e);
    Vertex newV = he.vertex();
    geometry->inputVertexPositions[newV] = newPos;

    std::cout << "After: " << geometry->inputVertexPositions.size() << " positions in inputVertexPositions" << std::endl;
    std::cout << "After: " << mesh->nVertices() << " vertices according to nVertices()" << std::endl;

    geometry->refreshQuantities();

    // Register the mesh with polyscope
    psMesh = polyscope::registerSurfaceMesh(
        "meshWithSplit",
        geometry->inputVertexPositions, mesh->getFaceVertexList(),
        polyscopePermutations(*mesh));
  }
}

I replaced the sample myCallback() function in the current version of the template project (https://github.com/nmwsharp/gc-polyscope-project-template) with this one. When I run the project on any mesh and click the resulting "Split one edge" button, the following console output results:

The inputVertexPositions buffer has now doubled in size (not capacity, but actual reported size), but nVertices() indicates that only one vertex was actually added. The resulting mismatch produces an error when attempting to reregister the mesh in polyscope.

It looks like the polyscope error only started being raised after this commit, but the underlying behavior has been silently occurring before then. Prior to that commit, this was what would happen:


Polyscope would think that there were twice as many vertices, but it would allow the registration to go through without an error. So the issue would only be noticed if one was paying attention to the vertex count on the UI (which apparently I hadn't been).

There's a super old comment in the intrinsic geometry implementation

  // ONEDAY try these for better accuracy in near-degenerate triangles?
  // "Miscalculating Area and Angles of a Needle-like Triangle" https://www.cs.unc.edu/~snoeyink/c/c205/Triangle.pdf

Just turning this comment into an issue, in the hopes that some enterprising young student might take a stab at it!

Hi,

A while back I came across this paper : Discrete Differential Forms for Computational Modeling
The paper led me to study smooth manifolds and some fundamentals of Riemannian Geometry. Although the paper explains how to implement discrete differential forms I still can't quite match it with the theory.

Can you point in the code where is the implementation of the differential forms, hodge dual etc in the code so I can actually dig into it and see how this is done in practice.

I also would like if formalizing the algorithms in such language actually make it easier the implementation or does it simplify maybe the formalization of the algorithm itself?

This is a minor feature request which I've found myself wanting a few times. I wish mesh elements had a reinterpretTo function similar to the mesh data function. For example, I often have a pair of meshes with the same vertex set, and it would be convenient to say

Vertex v1 = // some vertex from mesh1
Vertex v2 = v1.reinterpretTo(mesh2);

It seems like the functions equivalentPointOnIntrinsic(const SurfacePoint& pointOnInput) and equivalentPointOnInput(const SurfacePoint& pointOnIntrinsic) may be declared but not defined for SignpostIntrinsicTriangulation? When I grep the master branch I seem to find their declarations in signpost_intrinisc_triangulation.h, but I can't seem to locate a definition. Obviously, my first clue was unresolved external link errors when I trie to use them. Just curious if their implementation might have existed at one time?

I am trying to run the pointcloud examples from http://geometry-central.net/pointcloud/basics. I am unable to compile the examples, though, because it cannot find the pointcloud headers. I am able to compile with surface or utility headers, so it is able to access the geometry-central library.

I noticed that the geometry-central src CMakeLists does not include pointcloud headers. Even after adding them, the pointcloud examples still don't compile.

In geometry processing, it's common to solve problems like Δu = f, where u and f are vector-valued. In the discrete setting this problem translates to a matrix equation where the entries are vector-valued. There are various ways to translate these problems into the scalar problems already supported by GeometryCentral; would be nice to encapsulate this translation in a friendly wrapper.

It would be nice if it's easier to write "simple" C-style code in geometry-central without thinking about e.g. std::unqiue_ptr<>. Some changes that would push us in that direction without breaking existing code:

• Add copy/move constructors to the heavy classes like SurfaceMesh and Geometry subtypes. There are some pitfalls here, but that's outweighted by the benefit of avoiding pointer contortions.
• Add versions of the factory functions (loaders / copiers) which return either raw pointers or values

Hi,

First of all, great work on this library !

I'm trying to do some simple mutation, like adding a vertex along an edge. I have this geometry in the image below which could be seen as 2 rectangles merged together, and I have different CornerData (let's call them UVs) for each corner of the rectangles.

When I try to add a vertex along the shared edge, it seems that some CornerData are lost (showed with the ?).

Is this normal behavior ? How could I change this besides copying the CornerData before the mutation and rewrite them later ?

Right now I do :

//Store UVs
Vector2 previousUVF1 = uvData[edge.halfedge().corner()];

//Insert new vertex, and get the halfedge in the center point, then rewrite UVs
auto newhe = mesh->insertVertexAlongEdge(edge);
uvData[newhe.prevOrbitFace().corner()] = previousUVF1;

(And this is not related to this issue, but will mutation like adding a vertex along an edge be available for non manifold meshes ?)

Thanks again for you work.

-Edit : after further investigation, it seems that only one corner lose his data. Update the image to reflect that

Currently, the smoothest vertex direction fields algorithms do not precisely implement the energy from the Globally Optimal Direction Fields paper; they use a slightly different, simpler algorithm.

The Globally Optimal Direction Fields paper carefully derives a FEM expression which incorporates curvature (Eqn 18). The geometry-central implementation instead just uses a simple connection Laplacian which does not really incorporate curvature properly. This is mostly fine most of the time, but some examples expose the distinction. For instance, try computing a curvature-aligned cross field on a uniform sphere.

Fix: implement the proper energy.

When linking to geometry-central on macOS (e.g. using https://github.com/nmwsharp/gc-polyscope-project-template), cmake/make is complaining about -lm and -dl targets:

[ 95%] Built target geometry-central
make[2]: *** No rule to make target `../deps/geometry-central/deps/-lm', needed by `bin/gc_project'.
make[2]: *** No rule to make target `../deps/geometry-central/deps/-ldl', needed by `bin/gc_project'.
make[2]: Target `CMakeFiles/gc_project.dir/build' not remade because of errors.

The problem comes from path variables in cmake and when I comment the line 168 of geometry-central/deps/CMakeLists.txt:

   convert_paths_to_absolute(GC_DEP_LIBRARIES)

==>

  #convert_paths_to_absolute(GC_DEP_LIBRARIES)

the problem disappears.

(and no problem on linux)

Attached below is the code snippet that reproduces the problem.

#include <geometrycentral/surface/manifold_surface_mesh.h>
#include <geometrycentral/surface/vertex_position_geometry.h>
#include <geometrycentral/surface/simple_polygon_mesh.h>
#include <geometrycentral/surface/surface_mesh_factories.h>
#include <geometrycentral/utilities/vector3.h>

#include "polyscope/polyscope.h"
#include "polyscope/surface_mesh.h"

namespace gc = ::geometrycentral;
namespace gcs = ::geometrycentral::surface;

std::tuple<std::unique_ptr<gcs::ManifoldSurfaceMesh>,
           std::unique_ptr<gcs::VertexPositionGeometry>>
diamond(double dihedral) {
  std::vector<gc::Vector3> coords;
  std::vector<std::vector<std::size_t>> polygons;

  // Initialize vertex coordinates
  auto makeVertex = [](double x, double y, double z) -> gc::Vector3 {
    return gc::Vector3{std::move(x), std::move(y), std::move(z)};
  };

  coords.emplace_back(makeVertex(0, 0.5, 0));
  coords.emplace_back(makeVertex(0, -0.5, 0));
  coords.emplace_back(makeVertex(cos(dihedral) * std::sqrt(3) / 2, 0,
                                 -sin(dihedral) * std::sqrt(3) / 2));
  coords.emplace_back(makeVertex(-std::sqrt(3) / 2, 0, 0));
  // Initialize Faces
  polygons.emplace_back(std::vector<std::size_t>{0, 1, 2});
  polygons.emplace_back(std::vector<std::size_t>{0, 3, 1});

  gcs::SimplePolygonMesh soup(polygons, coords);
  soup.mergeIdenticalVertices();
  return gcs::makeManifoldSurfaceMeshAndGeometry(soup.polygons,
                                                 soup.vertexCoordinates);
}

int main() {
  // Initialize mesh and geometry
  std::unique_ptr<gcs::ManifoldSurfaceMesh> mesh;
  std::unique_ptr<gcs::VertexPositionGeometry> vpg;
  std::tie(mesh, vpg) = diamond(1);
  vpg->requireEdgeDihedralAngles();

  // split edge
  for (gcs::Edge e : mesh->edges()) {
    if (!e.isBoundary()) {
      gcs::Halfedge newhe = mesh->splitEdgeTriangular(e);
    }
  }

  // update edge dihedral angles
  mesh->compress();
  vpg->refreshQuantities();

  /// Initiate in polyscope
  polyscope::init();
  polyscope::SurfaceMesh *polyMesh = polyscope::registerSurfaceMesh(
      "Membrane", vpg->inputVertexPositions, mesh->getFaceVertexList());
  polyscope::getSurfaceMesh("mesh")
      ->addEdgeScalarQuantity("edge_dihedral", vpg->edgeDihedralAngles);
  polyscope::show();

  return 0;
}

This leads to the visualization in Polysocpe:
Before mutation:

After mutation:

I am also tagging Chris Here @ctlee.

We should add support for the following formats (both reading and writing):

• spmat (MATLAB) — 1-based indexing, just a list of (row, column, value) triples. No matrix size. No special assumptions about ordering of triples. Redundant triples are allowed on read. For matrices with an empty last row or column, we should be sure to write a "dummy" zero in the bottom-right entry.
• npz (SciPy) — see https://docs.scipy.org/doc/scipy-1.2.3/reference/generated/scipy.sparse.load_npz.html
• dmat (libigl) — see https://libigl.github.io/libigl-python-bindings/igl_docs/#read_dmat First line gives matrix dimensions; remaining lines give (row, column, value) entries. Not clear if entries are 0-or 1-based. Seems a special (sorted) ordering is assumed, probably to provide fast load into an Eigen matrix without sorting.

Hi,

I've been searching for an implementation of the heat method that would be easily usable on intrinsic geometry and your library seems exactly what I was looking for.

However, I'm wondering whether the data I'm using may not be suitable.
Does the implicit geometry need to be watertight?

The real data I'm looking to apply geodesic distance computations only has local intrinsic information (notably edge lengths), and no global embedding of the vertices.
I have the whole global topology, but it's made of multiple groups of local 2D embeddings, and thus I was hoping your library would do the trick (especially since it seems pretty easy to compile it to Wasm and I'm planning to use it in Javascript).

However, here's a test which I was hoping would pass but currently ends up crashing:

#include "geometrycentral/surface/surface_mesh.h"
#include "geometrycentral/utilities/mesh_data.h"
#include "geometrycentral/surface/edge_length_geometry.h"
#include "geometrycentral/surface/heat_method_distance.h"
#include <Eigen/Core>
#include <stdio.h>
#include <iostream>
#include <math.h>

using namespace geometrycentral;
using namespace geometrycentral::surface;

extern "C" {

  int main(){
    // create quad with two triangles in CCW order
    //
    //   0--3
    //   | /|
    //   |/ |
    //   1--2
    //

    const size_t F = 2;
    Eigen::MatrixX3i faces(F, 3);
    faces << 0, 1, 3,
             1, 2, 3;
    
    Eigen::MatrixX3d edges(F, 3);
    edges << 1, M_SQRT2, 1,
             1, 1, M_SQRT2;

    std::cout << "Faces:\n" << faces << "\n";
    std::cout << "Edges:\n" << edges << "\n";

    // create surface mesh
    std::unique_ptr<SurfaceMesh> mesh;
    mesh.reset(new SurfaceMesh(faces));
    mesh->compress();
    mesh->printStatistics();

    // create geometry data
    EdgeData<double> edgeLengths(*mesh);
    for(size_t i = 0; i < faces.rows(); ++i){
      Face f = mesh->face(i);
      Halfedge he = f.halfedge(); edgeLengths[he.edge()] = edges(i, 0);
      he = he.next(); edgeLengths[he.edge()] = edges(i, 1);
      he = he.next(); edgeLengths[he.edge()] = edges(i, 2);
    }
    std::cout << "Edges:\n" << edgeLengths.raw() << "\n";
    for(Edge e : mesh->edges()){
      printf("Edge #%zu = %g\n", e.getIndex(), edgeLengths[e]);
    }
    std::unique_ptr<EdgeLengthGeometry> geometry;
    geometry.reset(new EdgeLengthGeometry(*mesh, edgeLengths));

    std::cout << "Precomputation\n";

    // create heat method distance solver (precomputation happens here)
    std::unique_ptr<HeatMethodDistanceSolver> heatSolver;
    heatSolver.reset(new HeatMethodDistanceSolver(*geometry));

    std::cout << "Solve\n";

    // solve for distance from first vertex
    const Vertex v = mesh->vertex(0);
    VertexData<double> distToSource = heatSolver->computeDistance(v);
    std::cout << "Distances:\n" << distToSource.raw() << "\n";

    return 0;
  }

}

My output (through node using Wasm) shows this:

Faces:
0 1 3
1 2 3
Edges:
      1 1.41421       1
      1       1 1.41421
Halfedge mesh with: 
    # verts =  4
    # edges =  5
    # faces =  2
    # halfedges =  6  (6 interior, 0 exterior)
      and 0 boundary components. 
Edges:
      1
1.41421
      1
      1
      1
      0
      0
      0
Edge #0 = 1
Edge #1 = 1.41421
Edge #2 = 1
Edge #3 = 1
Edge #4 = 1
Precomputation
Solve
exception thrown: 5310560

Thus the face and edge data seems appropriate (except possibly for the edges data having a few zeros for additional edges, but I guess that's due to some memory allocation strategy).
The precomputation seems also fine (initially, I erroneously used M_SQRT1_2, which leads to issues during the precomputation already, but I expect that's due to the math behind since such shorter edge may lead to trouble during factorization).

My current expectation is that the distances on the square are going to be 1 from the source, except for the diagonal that ends up being again 1.41421.
Am I making wrong assumptions on the data structure or expectations as to what the implicit geometry requires?

Thanks for any clarifications.

Hi, I found this project really interesting.

I want to use halfedge mesh for Finite element and Jacobian calculation on 2D domain, where $\Delta \sigma \Delta f = 0$ and a boundary valued problem is solved. The Jacobian is defined as the first order derivative with respect to the conductivity $\sigma$.

Can geometry central represent 2D meshes (triangles) efficiently? as well as the calculation on 2D meshes?

Best.

Have you considered providing python wrapper of geometry-central?
Yes - I realize it would be a significant amount of work :)

Could a method be added to the heat solver so that it can use an existing VertexData<double> object? When calling this function many times it would be useful to reuse an existing array. Something like:

void HeatMethodDistanceSolver::computeDistance(Vertex v, VertexData<double>& out)

As far as I can tell, there are not yet implemented standard mesh processing algorithms, like:

• decimation (quadratic-error based simplification)
• hole filling
• components detection
• self-intersecting faces detection (preferably correction)
• degenerate face removal
• non-manifold correction (though can be easily implemented as far as I understand using a succession of non-manifold edge and vertex fixing)

Do you have plans to add them at some point?

Would there be any interest extending the edge flip method to provide a gift wrappers algorithm over surfaces? There don't appear to be any implementations online I can find - would throw in a vote for it being a neat feature to have :)

Can possibly submit a PR at a later date myself if it's of interest!

Hi,

Can geometry central given source and destination point on a 3D mesh represent somehow the geodesic path?
I can see this : https://github.com/nmwsharp/geometry-central/blob/master/include/geometrycentral/surface/exact_polyhedral_geodesics.h

But the description isn't to me 100% clear of what it does, it sounds however it might do the job.
Thank you

After #52 SurfaceMesh::flipEdge() now looks like:

  bool flip(Edge e, bool preventSelfEdges = true);

This gets at a design question: how should we expose behavior options?

Two possible strategies are:

• Options per-method to control behavior (as above)
• Stateful (mutable or immutable) flags on the class, which set a consistent behavior for all such function calls.

We should think about which of these strategies we want to pursue, and what the consequences would be.

Some motivating questions to consider:

• What is the most natural for a user? What minimizes verbosity-creep?
• What is least likely to lead to bugs, and easiest to debug?
• What scales well as these classes grow in complexity and functionality?
• What is not-too-painful to implement internally?

gcc5 loves to erroneously reject code with templates/tuples/initializers. However, a decent number of people still use it, as it is the default in Ubuntu 16.04. We should make sure our code passes on gcc5.

TODOs

• Add CI testing to catch gcc5 issues, similar to nmwsharp/polyscope@d5ca9fa
• Will probably also need to clean up some small failing code snippets

Thanks for the great library! I'm just starting to work through it. I noticed along the way that the documentation page lists:

size_t HalfedgeMesh::nBoundaryVertices()
Returns the number of vertices incident on the boundary.

The only reference I can find to this method though is in the docs.

Perhaps the method was unimplemented, renamed, or removed?

Problem:
Right now, Eigen is a submodule. This is bad for two reasons:

• The Eigen git repo is pretty giant: 100MB, which dwarfs everything else in this project
• Some users might very reasonably want to use a version of Eigen already in their codebase

I'm against removing the Eigen from the repo and requiring that users manually acquire Eigen. It's true that for many C++ geometry programmers Eigen is very familiar, but nonetheless preserving a download-build-run workflow for new novice users is very important.

Unfortunately, there's no such thing as an "optional" submodule in git, selectively downloading only some submodule is not really feasible.

Solution:
Remove Eigen as a submodule, and download Eigen via CMake script at configure time. We can set options so Eigen is only downloaded if not automatically resolved, and since we don't have the whole git history the source size should be smaller.

After installation on Ubuntu 18.04.3 following the tutorial from http://geometry-central.net/build/building/, there is no geometry.h header in geometrycentral/surface directory as required by example code. Here is the contents of geometrycentral/surface:

barycentric_coordinate_helpers.h    halfedge_element_types.h	    ply_halfedge_mesh_data.h
barycentric_coordinate_helpers.ipp  halfedge_element_types.ipp	    ply_halfedge_mesh_data.ipp
base_geometry_interface.h	    halfedge_factories.h	    polygon_soup_mesh.h
detect_symmetry.h		    halfedge_iterators.h	    signpost_intrinsic_triangulation.h
direction_fields.h		    halfedge_iterators.ipp	    signpost_intrinsic_triangulation.ipp
edge_length_geometry.h		    halfedge_logic_templates.ipp    surface_centers.h
edge_length_geometry.ipp	    halfedge_mesh.h		    surface_point.h
embedded_geometry_interface.h	    halfedge_mesh.ipp		    surface_point.ipp
exact_polyhedral_geodesics.h	    heat_method_distance.h	    trace_geodesic.h
extrinsic_geometry_interface.h	    intrinsic_geometry_interface.h  vector_heat_method.h
fast_marching_method.h		    mesh_graph_algorithms.h	    vertex_position_geometry.h
halfedge_containers.h		    meshio.h			    vertex_position_geometry.ipp
halfedge_containers.ipp		    mesh_ray_tracer.h

It's fairly common (e.g., in GLSL) to support "swizzled" access to vector components. E.g.,

u.yz = v.xy;

will copy the x and y components of v into the y and z components of u, respectively. A common use case is translating between 2D and 3D vectors, e.g., mapping UV coordinates into three-dimensional space, or translating 3D coordinates of a planar mesh into UV coordinates. Etc.

There are all sorts of hacks to achieve GLSL-like syntax, but it's probably just as good to simply add some member functions. E.g., the line above might become

u.yz() = v.xy();

The right-hand side is easy enough to implement: just return a Vector2 build from the appropriate components. The left-hand side would require a bit of hacking. For instance, you could return a Swizzle type that internally holds references to components of the calling vector (in this case, u), and an assignment operator that copies values into these referenced components.

Even just implementing the simpler right-hand sides would be useful. Plenty of times you just want to do stuff like

Vector2 u;
Vector3 p;
u = p.xy();

(The place where this really matters is when you're extracting components inline in some larger computation. Otherwise, sure, you can just write u = Vector2{ p.x, p.y };)

Hi Folks,
Would it make sense to add stl iterator_categories (and all the other sundry stl iterator typedefs) for the RangeSetBase class? There are times when I would like to use stl algorithms on some entitie of the mesh, like:

std::transform(begin(mesh.faces()), end(mesh.faces()),...)

But, it seems we need category support for that. My guess is that std::forward_iterator_tag would be the most likely candidate, but maybe std::bidirectional_iterator_tag would work too?

Ultimately, I would like to be able to use the parallel version of the stl algorithms, which I currently do in some code. But, I end up first copying the entity handles to a std::vector with a range based for loop, and then use that sequence as my iteration range for the algorithm. From a performance standpoint, the penalty of copying is probably in the noise for cases where parallelizing things is a benefit, but it still seems wasteful. I only bring that up because I don't quite know all the details for how the parallel versions of the std algorithms partition things, etc... Therefore, I don't know how smart/safe/etc... a given iterator/container pair need to actually be to work correctly in the parallel versions of the algorithms. A vector of handles seems to work fine, but obviously that is the least restrictive container/iterator type.

As a last thought, this may also play into enhancing general const-correctness for the overall api as well. However, I seem to gather from #66 it is not as trival as sprinkling const everywhere.

I am trying to do arithmetic on vertices, such as finding the distance between two vertices, but cannot find a method that is supported. geometrycentral::surface::Vertex does not seem to define vertex coordinates, nor are arithmetic operations defined. I also tried converting to geomectrycentral::Vector3, but that does not seem possible either. Is there a way to do this?

Currently VertexPositionGeometry contains both vertexPositions and inputVertexPositions, which arose from a well-intentioned design but turns out to be pretty annoying/confusing in most common use cases.

Proposal: Have just one copy of the vertex positions, vertexPositions.

Probably there are quite a few consequences to work through here... but it feels like the right direction moving forward.