Program Listing for File SurfaceMesh.h

↰ Return to documentation for file (include/gamer/SurfaceMesh.h)

// This file is part of the GAMer software.
// Copyright (C) 2016-2021
// by Christopher T. Lee and contributors
//
// This library is free software; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public
// License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, see <http://www.gnu.org/licenses/>
// or write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330,
// Boston, MA 02111-1307 USA

#pragma once

#include <iostream>
#include <memory>
#include <stdexcept>
#include <string>
#include <unordered_set>
#include <utility>

#include <Eigen/Dense>
#include <Eigen/Eigenvalues>
#include <casc/casc>

#include "gamer/Vertex.h"

namespace gamer {

struct SMGlobal {
  int marker;
  float volumeConstraint;
  bool useVolumeConstraint;
  bool ishole;

  SMGlobal(int marker = -1, float volumeConstraint = -1,
           bool useVolumeConstraint = false, bool ishole = false)
      : marker(marker), volumeConstraint(volumeConstraint),
        useVolumeConstraint(useVolumeConstraint), ishole(ishole) {}
};

struct SMVertex : Vertex {
  Vector normal;
  using Vertex::Vertex;
};

struct SMEdge {
  bool selected;

  SMEdge() : SMEdge(0) {}

  SMEdge(bool select) : selected(select) {}
};

struct SMFaceProperties {
  int marker;
  bool selected;
  Vector normal;
  SMFaceProperties(int marker, bool selected)
      : marker(marker), selected(selected) {}
};

struct SMFace : casc::Orientable, SMFaceProperties {
  SMFace() : SMFace(Orientable{0}, SMFaceProperties{-1, false}) {}

  SMFace(int marker, bool selected)
      : SMFace(Orientable{0}, SMFaceProperties{marker, selected}) {}

  SMFace(int orient, int marker, bool selected)
      : SMFace(Orientable{orient}, SMFaceProperties{marker, selected}) {}

  SMFace(Orientable orient, SMFaceProperties prop)
      : Orientable(orient), SMFaceProperties(prop) {}

  friend std::ostream &operator<<(std::ostream &output, const SMFace &f) {
    output << "SMFace("
           << "m:" << f.marker << ";sel:" << std::boolalpha << f.selected
           << ";o:" << f.orientation << ")";
    return output;
  }

  std::string to_string() const {
    std::ostringstream output;
    output << *this;
    return output.str();
  }
};

namespace surfmesh_detail {
struct surfmesh_traits {
  using KeyType = int;
  using NodeTypes = util::type_holder<SMGlobal, SMVertex, SMEdge, SMFace>;
  using EdgeTypes =
      util::type_holder<casc::Orientable, casc::Orientable, casc::Orientable>;
};
} // end namespace surfmesh_detail

using SurfaceMesh = casc::simplicial_complex<surfmesh_detail::surfmesh_traits>;

tensor<double, 3, 2> getTangent(const SurfaceMesh &mesh,
                                SurfaceMesh::SimplexID<1> vertexID);

tensor<double, 3, 2> getTangent(const SurfaceMesh &mesh,
                                SurfaceMesh::SimplexID<3> faceID);

Vector getNormalFromTangent(const tensor<double, 3, 2> tangent);

Vector getNormal(const SurfaceMesh &mesh, SurfaceMesh::SimplexID<1> vertexID);

Vector getNormal(const SurfaceMesh &mesh, SurfaceMesh::SimplexID<3> faceID);

namespace surfacemesh_detail {
void decimateVertex(SurfaceMesh &mesh, SurfaceMesh::SimplexID<1> vertexID,
                    std::size_t rings = 2);

tensor<double, 3, 2>
computeLocalStructureTensor(const SurfaceMesh &mesh,
                            const SurfaceMesh::SimplexID<1> vertexID,
                            const int rings);

tensor<double, 3, 2>
computeLSTFromCache(const SurfaceMesh &mesh,
                    const SurfaceMesh::SimplexID<1> vertexID, const int rings);

template <std::size_t dimension>
auto getTangentH(const SurfaceMesh &mesh,
                 const tensor<double, dimension, 1> &origin,
                 SurfaceMesh::SimplexID<SurfaceMesh::topLevel> curr) {
  return (*curr).orientation;
}

template <std::size_t level, std::size_t dimension>
auto getTangentH(const SurfaceMesh &mesh,
                 const tensor<double, dimension, 1> &origin,
                 SurfaceMesh::SimplexID<level> curr) {
  tensor<double, dimension, SurfaceMesh::topLevel - level> rval;
  auto cover = mesh.get_cover(curr);
  for (auto alpha : cover) {
    auto edge = *mesh.get_edge_up(curr, alpha);
    const auto &v = (*mesh.get_simplex_up({alpha})).position;
    auto next = mesh.get_simplex_up(curr, alpha);
    rval += edge.orientation * (v - origin) * getTangentH(mesh, origin, next);
  }
  return rval / cover.size();
}

template <std::size_t dimension>
auto getTangentF(const SurfaceMesh &mesh,
                 const tensor<double, dimension, 1> &origin,
                 SurfaceMesh::SimplexID<SurfaceMesh::topLevel> curr,
                 std::set<SurfaceMesh::KeyType> &cover) {
  return (*curr).orientation;
}

template <std::size_t level, std::size_t dimension>
auto getTangentF(const SurfaceMesh &mesh,
                 const tensor<double, dimension, 1> &origin,
                 SurfaceMesh::SimplexID<level> curr,
                 std::set<SurfaceMesh::KeyType> &cover) {
  tensor<double, dimension, SurfaceMesh::topLevel - level> rval;
  for (auto alpha : cover) {
    auto edge = *mesh.get_edge_up(curr, alpha);
    const auto &v = (*mesh.get_simplex_up({alpha})).position;
    auto next = mesh.get_simplex_up(curr, alpha);
    auto coverup = cover;
    coverup.erase(alpha);
    rval += edge.orientation * (v - origin) *
            getTangentF(mesh, origin, next, coverup);
  }
  return rval / cover.size();
}

template <class K> struct initLocalOrientation {};

template <std::size_t k>
struct initLocalOrientation<std::integral_constant<std::size_t, k>> {
  template <typename Iterator>
  static void apply(SurfaceMesh &mesh, const std::set<int> &&names,
                    Iterator begin, Iterator end) {
    std::vector<SurfaceMesh::SimplexID<k + 1>> next;
    for (auto curr = begin; curr != end; ++curr) {
      auto currSimplexID = *curr;
      for (auto a : mesh.get_cover(currSimplexID)) {
        // Look for key a in names
        auto find = names.find(a);
        if (find != names.end()) {
          next.push_back(mesh.get_simplex_up(currSimplexID, a));

          int orient = 1;
          for (auto b : mesh.get_name(currSimplexID)) {
            if (a > b) {
              if (a > b) {
                orient *= -1;
              } else {
                break;
              }
            }
          }
          (*mesh.get_edge_up(currSimplexID, a)).orientation = orient;
        }
      }
    }
    initLocalOrientation<std::integral_constant<std::size_t, k + 1>>::apply(
        mesh, std::move(names), next.begin(), next.end());
  }
};

template <>
struct initLocalOrientation<
    std::integral_constant<std::size_t, SurfaceMesh::topLevel>> {
  template <typename Iterator>
  static void apply(SurfaceMesh &mesh, const std::set<int> &&names,
                    Iterator begin, Iterator end) {}
};

bool computeLocalOrientation(
    SurfaceMesh &mesh, const std::vector<SurfaceMesh::SimplexID<2>> &edgeList);

void weightedVertexSmooth(SurfaceMesh &mesh, SurfaceMesh::SimplexID<1> vertexID,
                          int rings);

Vector weightedVertexSmoothCache(SurfaceMesh &mesh,
                                 SurfaceMesh::SimplexID<1> vertexID,
                                 std::size_t rings);

void barycenterVertexSmooth(SurfaceMesh &mesh,
                            SurfaceMesh::SimplexID<1> vertexID);

void edgeFlip(SurfaceMesh &mesh, SurfaceMesh::SimplexID<2> edgeID);

void edgeFlipCache(SurfaceMesh &mesh, SurfaceMesh::SimplexID<2> edgeID);

template <class Inserter>
void selectFlipEdges(
    const SurfaceMesh &mesh, bool preserveRidges,
    std::function<bool(const SurfaceMesh &, const SurfaceMesh::SimplexID<2> &)>
        &&checkFlip,
    Inserter iter) {
  casc::NodeSet<SurfaceMesh::SimplexID<2>> ignoredEdges;

  for (auto edgeID : mesh.get_level_id<2>()) {
    if ((*edgeID).selected == true) {
      if (!ignoredEdges.count(edgeID)) {
        auto up = mesh.get_cover(edgeID);
        // The mesh is not a surface mesh...
        if (up.size() > 2) {
          // std::cerr << "This edge participates in more than 2
          // faces. "
          //           << "Returning..." << std::endl;
          gamer_runtime_error("SurfaceMesh is not pseudomanifold. Found "
                                   "an edge connected to more than 2 faces.");
        } else if (up.size() < 2) // Edge is a boundary
        {
          // std::cerr << "This edge participates in fewer than 2
          // faces. "
          //           << "Returning..." << std::endl;
          continue;
        }

        // Check if the edge is a part of a tetrahedron.
        if (mesh.exists<2>({up[0], up[1]})) {
          // std::cerr << "Found a tetrahedron cannot edge flip."
          //           << std::endl;
          continue;
        }

        // Check if we're on a ridge. This prevents folding also.
        if (preserveRidges) {
          auto a = getNormal(mesh, mesh.get_simplex_up(edgeID, up[0]));
          auto b = getNormal(mesh, mesh.get_simplex_up(edgeID, up[1]));
          auto val = angle(a, b);
          if (val > 60) {
            continue;
          }
        }

        // Check the flip using user function
        if (checkFlip(mesh, edgeID)) {
          *iter++ = edgeID; // Insert into edges to flip

          // The local topology will be changed by edge flip.
          // Don't flip edges which share a common face.
          std::set<SurfaceMesh::SimplexID<2>> tmpIgnored;
          kneighbors(mesh, edgeID, 3, tmpIgnored);
          ignoredEdges.insert(tmpIgnored.begin(), tmpIgnored.end());

          // Local neighborhood append. Larger neighborhood selected
          // above appears to work better...
          // neighbors(mesh, edgeID, std::inserter(ignoredEdges,
          // ignoredEdges.end()));
        }
      }
    }
  }
}

bool checkEdgeFlip(
    const SurfaceMesh &mesh, bool preserveRidges,
    SurfaceMesh::SimplexID<2> edgeID,
    std::function<bool(const SurfaceMesh &, const SurfaceMesh::SimplexID<2> &)>
        &&checkFlip);
bool checkFlipAngle(const SurfaceMesh &mesh,
                    const SurfaceMesh::SimplexID<2> &edgeID);

int checkFlipValenceExcess(const SurfaceMesh &mesh,
                           const SurfaceMesh::SimplexID<2> &edgeID);

void normalSmoothH(SurfaceMesh &mesh, const SurfaceMesh::SimplexID<1> vertexID,
                   const double k);

void findHoles(const SurfaceMesh &mesh,
               std::vector<std::vector<SurfaceMesh::SimplexID<2>>> &holeList);

bool orderBoundaryEdgeRing(
    const SurfaceMesh &mesh,
    std::set<SurfaceMesh::SimplexID<2>> &unvisitedBdryEdges,
    std::vector<SurfaceMesh::SimplexID<1>> &visitedVerts,
    std::vector<SurfaceMesh::SimplexID<2>> &bdryRing);

void edgeRingToVertices(
    const SurfaceMesh &mesh, std::vector<SurfaceMesh::SimplexID<2>> &edgeRing,
    std::back_insert_iterator<std::vector<SurfaceMesh::SimplexID<1>>> iter);

void triangulateHoleHelper(
    SurfaceMesh &mesh, std::vector<SurfaceMesh::SimplexID<1>> &boundary,
    const SMFace &fdata,
    std::back_insert_iterator<std::vector<SurfaceMesh::SimplexID<2>>> iter);

void triangulateHole(SurfaceMesh &mesh,
                     std::vector<SurfaceMesh::SimplexID<1>> &sortedVerts,
                     const SMFace &fdata,
                     std::vector<SurfaceMesh::SimplexID<2>> &edgeList);

template <typename Complex> struct CopyHelper {
  using SimplexSet = typename casc::SimplexSet<Complex>;
  using KeyType = typename Complex::KeyType;

  template <std::size_t k>
  static void apply(Complex &before, Complex &after, const SimplexSet &S) {
    for (auto sID : casc::get<k>(S)) {
      auto name = before.get_name(sID);
      auto data = *sID;
      after.insert(name, data);
    }
  }
};

template <typename Iterator>
void vertexGrabber(const SurfaceMesh &F, int need,
                   std::vector<SurfaceMesh::SimplexID<1>> &nbors,
                   Iterator begin, Iterator end) {
  if (need <= 0)
    return;
  std::set<SurfaceMesh::SimplexID<1>> next;
  for (; begin != end; ++begin) {
    for (auto a : F.get_cover(*begin)) {
      auto id = F.get_simplex_up(*begin, a);
      for (auto b : F.get_name(id)) {
        if (b != a) {
          auto nbor = F.get_simplex_down(id, b);
          if (std::find(nbors.begin(), nbors.end(), nbor) == nbors.end()) {
            // Haven't visited so push into next ring and
            nbors.push_back(nbor);
            next.insert(nbor);
            if (--need == 0)
              return;
          }
        }
      }
    }
  }
  vertexGrabber(F, need, nbors, next.begin(), next.end());
}

template <class Complex>
void vertexGrabber(const Complex &F, int need,
                   std::vector<typename Complex::template SimplexID<1>> &nbors,
                   typename Complex::template SimplexID<1> vid) {
  if (need <= 0)
    return;
  std::set<SurfaceMesh::SimplexID<1>> next;
  for (auto a : F.get_cover(vid)) {
    auto id = F.get_simplex_up(vid, a);
    for (auto b : F.get_name(id)) {
      if (b != a) {
        auto nbor = F.get_simplex_down(id, b);
        if (std::find(nbors.begin(), nbors.end(), nbor) == nbors.end()) {
          // Haven't visited so push into next ring and
          nbors.push_back(nbor);
          next.insert(nbor);
          if (--need == 0)
            return;
        }
      }
    }
  }
  vertexGrabber(F, need, nbors, next.begin(), next.end());
}
} // end namespace surfacemesh_detail

std::unique_ptr<SurfaceMesh> readOFF(const std::string &filename);

void writeOFF(const std::string &filename, const SurfaceMesh &mesh);

std::unique_ptr<SurfaceMesh> readOBJ(const std::string &filename);

void writeOBJ(const std::string &filename, const SurfaceMesh &mesh);

void writeComsol(const std::string &filename, const SurfaceMesh &mesh);

void writeComsol(const std::string &filename,
                 const std::vector<SurfaceMesh const *> &meshes);

void print(const SurfaceMesh &mesh);

void printQualityInfo(const std::string &filename, const SurfaceMesh &mesh);

void generateHistogram(const SurfaceMesh &mesh);

std::tuple<double, double, int, int> getMinMaxAngles(const SurfaceMesh &mesh,
                                                     double maxMinAngle,
                                                     double minMaxAngle);

double getArea(const SurfaceMesh &mesh);

double getArea(const SurfaceMesh &mesh, SurfaceMesh::SimplexID<3> faceID);

double getArea(Vertex a, Vertex b, Vertex c);

double getArea(std::array<Vertex, 3> t);

double getVolume(const SurfaceMesh &mesh);

bool hasHole(const SurfaceMesh &mesh);

inline std::size_t getValence(const SurfaceMesh &mesh,
                              const SurfaceMesh::SimplexID<1> vertexID) {
  return mesh.get_cover(vertexID).size();
}

std::tuple<
    REAL *, REAL *, REAL *, REAL *,
    std::map<typename SurfaceMesh::KeyType, typename SurfaceMesh::KeyType>>
curvatureViaMDSB(const SurfaceMesh &mesh);

std::tuple<
    REAL *, REAL *, REAL *, REAL *,
    std::map<typename SurfaceMesh::KeyType, typename SurfaceMesh::KeyType>>
curvatureViaJets(const SurfaceMesh &mesh, std::size_t dJet = 2,
                 std::size_t dPrime = 2);

// void osculatingJets(const SurfaceMesh&mesh, std::size_t dJet = 2, std::size_t
// dPrime = 2);

//
// @param      mesh  The mesh
// @param[in]  v     Displacement vector
//
void translate(SurfaceMesh &mesh, Vector v);
void translate(SurfaceMesh &mesh, double dx, double dy, double dz);
void scale(SurfaceMesh &mesh, Vector v);
void scale(SurfaceMesh &mesh, double sx, double sy, double sz);
void scale(SurfaceMesh &mesh, double s);

std::pair<Vector, double> getCenterRadius(SurfaceMesh &mesh);
void centeralize(SurfaceMesh &mesh);

void smoothMesh(SurfaceMesh &mesh, int maxIter, bool preserveRidges,
                std::size_t rings = 2, bool verbose = false);

void coarse(SurfaceMesh &mesh, double coarseRate, double flatRate,
            double denseWeight, std::size_t rings = 2, bool verbose = false);

void coarse_dense(SurfaceMesh &mesh, REAL threshold, REAL weight,
                  std::size_t rings = 2, bool verbose = false);

void coarse_flat(SurfaceMesh &mesh, REAL threshold, REAL weight,
                 std::size_t rings = 2, bool verbose = false);

void normalSmooth(SurfaceMesh &mesh, double k = 1.0);

void fillHoles(SurfaceMesh &mesh);

void flipNormals(SurfaceMesh &mesh);

std::unique_ptr<SurfaceMesh> refineMesh(const SurfaceMesh &mesh);

std::unique_ptr<SurfaceMesh> sphere(int order);

std::unique_ptr<SurfaceMesh> cube(int order);

std::vector<std::unique_ptr<SurfaceMesh>> splitSurfaces(SurfaceMesh &mesh);

void cacheNormals(SurfaceMesh &mesh);

std::tuple<bool, int, int, int> getBettiNumbers(SurfaceMesh &mesh);
} // end namespace gamer