flatsurf · saraedum · Jun 16, 2022 · Jun 17, 2022 · videlec · Aug 9, 2022
@@ -1,7 +1,7 @@
 /**********************************************************************
  *  This file is part of flatsurf.
  *
- *        Copyright (C) 2019-2021 Julian Rüth
+ *        Copyright (C) 2019-2022 Julian Rüth
  *
  *  Flatsurf is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU General Public License as published by
@@ -36,7 +36,7 @@
 
 namespace flatsurf {
 
-// A triangulated translation surface. For most purposes this is the central
+// A triangulated half-dilation surface. For most purposes this is the central
 // object of the flatsurf library.
 template <class T>
 class FlatTriangulation : public FlatTriangulationCombinatorics<FlatTriangulation<T>>,
@@ -48,9 +48,24 @@ class FlatTriangulation : public FlatTriangulationCombinatorics<FlatTriangulatio
   using Coordinate = T;
 
   FlatTriangulation() noexcept;
+
+  // Create a surface with the given combinatorics and assign the vectors to the half edges,
+  // one vector per half edge in the order 1, 2, ..., n, -1, -2, ..., -n.
+  // Optionally, only the vectors for the positive half edges 1, ..., n can be
+  // given. The vectors for the negative half edges are then taken to be their
+  // negatives, i.e., a translation surface is produced from the data.
   FlatTriangulation(FlatTriangulationCombinatorial &&, const std::vector<Vector<T>> &vectors);
+
   FlatTriangulation(FlatTriangulationCombinatorial &&, const std::function<Vector<T>(HalfEdge)> &vectors);
 
+  bool translation_surface() const;
+
+  bool half_translation_surface() const;
+
+  bool dilation_surface() const;
+
+  bool half_dilation_surface() const;
+
   // Create an independent clone of this triangulation that is built from the
   // same data. There is no copy-constructor since it is too likely that this
   // is would not update the associated HalfEdgeMaps in the way that the caller

@@ -33,7 +33,6 @@
 #include "external/rx-ranges/include/rx/ranges.hpp"
 #include "impl/chain.impl.hpp"
 #include "impl/chain_iterator.impl.hpp"
-#include "util/assert.ipp"
 #include "util/hash.ipp"
 
 namespace flatsurf {

@@ -73,6 +73,9 @@ using std::end;
 
 template <typename T>
 Deformation<FlatTriangulation<T>> FlatTriangulation<T>::operator+(const OddHalfEdgeMap<Vector<T>> &shift) const {
+  if (!translation_surface())
+    throw std::logic_error("operator+ only implemented for translation surfaces");
+
   // Half edges that collapse at the end of the shift.
   EdgeSet collapsing;
 
@@ -312,6 +315,9 @@ Deformation<FlatTriangulation<T>> FlatTriangulation<T>::operator+(const OddHalfE
 
 template <typename T>
 Deformation<FlatTriangulation<T>> FlatTriangulation<T>::eliminateMarkedPoints() const {
+  if (!translation_surface())
+    throw std::logic_error("eliminating marked points only implemented for translation surfaces");
+
   std::optional<HalfEdge> collapse;
 
   for (const auto &vertex : this->vertices()) {
@@ -409,6 +415,9 @@ Deformation<FlatTriangulation<T>> FlatTriangulation<T>::eliminateMarkedPoints()
 
 template <typename T>
 Vector<T> FlatTriangulation<T>::shortest() const {
+  if (!translation_surface())
+    throw std::logic_error("shortest() only implemented for translation surfaces");
+
   const auto edges = this->edges();
   Edge shortest = *std::min_element(begin(edges), end(edges), [&](const auto &a, const auto &b) {
     const Vector x = fromHalfEdge(a.positive());
@@ -440,11 +449,17 @@ Vector<T> FlatTriangulation<T>::shortest(const Vector<T> &direction) const {
 
 template <typename T>
 const Vector<T> &FlatTriangulation<T>::fromHalfEdge(const HalfEdge e) const {
+  if (!translation_surface())
+    throw std::logic_error("fromHalfEdge() only implemented for translation surfaces");
+
   return self->vectors->get(e);
 }
 
 template <typename T>
 const flatsurf::Vector<exactreal::Arb> &FlatTriangulation<T>::fromHalfEdgeApproximate(HalfEdge e) const {
+  if (!translation_surface())
+    throw std::logic_error("fromHalfEdgeApproximate() only implemented for translation surfaces");
+
   return self->approximations->get(e);
 }
 
@@ -477,6 +492,9 @@ FlatTriangulation<T> FlatTriangulation<T>::clone() const {
 
 template <typename T>
 Deformation<FlatTriangulation<T>> FlatTriangulation<T>::slit(HalfEdge slit) const {
+  if (!translation_surface())
+    throw std::logic_error("slit() only implemented for translation surfaces");
+
   const auto codomain = FlatTriangulation(
       static_cast<const FlatTriangulationCombinatorial &>(*this).slit(slit),
       [&](HalfEdge e) {
@@ -491,6 +509,9 @@ Deformation<FlatTriangulation<T>> FlatTriangulation<T>::slit(HalfEdge slit) cons
 
 template <typename T>
 Deformation<FlatTriangulation<T>> FlatTriangulation<T>::insertAt(HalfEdge &nextTo, const Vector<T> &vector) const {
+  if (!translation_surface())
+    throw std::logic_error("insertAt() only implemented for translation surfaces");
+
   LIBFLATSURF_CHECK_ARGUMENT(inSector(nextTo, vector), "vector must be contained in the sector next to the half edge");
 
   Deformation<FlatTriangulation<T>> deformation{*this};
@@ -628,6 +649,9 @@ Deformation<FlatTriangulation<T>> FlatTriangulation<T>::insertAt(HalfEdge &nextT
 
 template <typename T>
 void FlatTriangulation<T>::delaunay() {
+  if (!translation_surface())
+    throw std::logic_error("delaunay() only implemented for translation surfaces");
+
   bool isDelaunay;
   do {
     isDelaunay = true;
@@ -642,6 +666,9 @@ void FlatTriangulation<T>::delaunay() {
 
 template <typename T>
 DELAUNAY FlatTriangulation<T>::delaunay(const Edge edge) const {
+  if (!translation_surface())
+    throw std::logic_error("delaunay() only implemented for translation surfaces");
+
   // We could eventually use Vector::insideCircumcircle() so vectors can
   // provide optimized implementations of this. However, at the moment, it
   // does not seem worth it.
@@ -675,6 +702,9 @@ DELAUNAY FlatTriangulation<T>::delaunay(const Edge edge) const {
 
 template <typename T>
 T FlatTriangulation<T>::area() const {
+  if (!translation_surface())
+    throw std::logic_error("area() only implemented for translation surfaces");
+
   T area = T();
   for (auto e : this->halfEdges()) {
     if (this->boundary(e)) continue;
@@ -690,13 +720,19 @@ T FlatTriangulation<T>::area() const {
 
 template <typename T>
 FlatTriangulation<T> FlatTriangulation<T>::scale(const mpz_class &scalar) const {
+  if (!translation_surface())
+    throw std::logic_error("scale() only implemented for translation surfaces");
+
   return FlatTriangulation(static_cast<const FlatTriangulationCombinatorial &>(*this).clone(), [&](HalfEdge e) {
     return scalar * fromHalfEdge(e);
   });
 }
 
 template <typename T>
 bool FlatTriangulation<T>::convex(HalfEdge e, bool strict) const {
+  if (!translation_surface())
+    throw std::logic_error("convex() only implemented for translation surfaces");
+
   if (strict)
     return fromHalfEdge(this->previousAtVertex(e)).ccw(fromHalfEdge(this->nextAtVertex(e))) == CCW::COUNTERCLOCKWISE &&
            fromHalfEdge(this->previousAtVertex(-e)).ccw(fromHalfEdge(this->nextAtVertex(-e))) == CCW::COUNTERCLOCKWISE;
@@ -707,21 +743,33 @@ bool FlatTriangulation<T>::convex(HalfEdge e, bool strict) const {
 
 template <typename T>
 bool FlatTriangulation<T>::inSector(const HalfEdge sector, const Vector<T> &vector) const {
+  if (!translation_surface())
+    throw std::logic_error("inSector() only implemented for translation surfaces");
+
   return fromHalfEdge(sector).ccw(vector) != CCW::CLOCKWISE && fromHalfEdge(this->nextAtVertex(sector)).ccw(vector) == CCW::CLOCKWISE;
 }
 
 template <typename T>
 bool FlatTriangulation<T>::inSector(const HalfEdge sector, const Vertical<FlatTriangulation<T>> &vector) const {
+  if (!translation_surface())
+    throw std::logic_error("inSector() only implemented for translation surfaces");
+
   return inSector(sector, vector.vertical());
 }
 
 template <typename T>
 SaddleConnections<FlatTriangulation<T>> FlatTriangulation<T>::connections() const {
+  if (!translation_surface())
+    throw std::logic_error("connections() only implemented for translation surfaces");
+
   return SaddleConnections(*this);
 }
 
 template <typename T>
 bool FlatTriangulation<T>::operator==(const FlatTriangulation<T> &rhs) const {
+  if (!translation_surface())
+    throw std::logic_error("equality testing only implemented for translation surfaces");
+
   if (this == &rhs)
     return true;
 
@@ -739,6 +787,9 @@ bool FlatTriangulation<T>::operator==(const FlatTriangulation<T> &rhs) const {
 
 template <typename T>
 int FlatTriangulation<T>::angle(const Vertex &vertex) const {
+  if (!translation_surface())
+    throw std::logic_error("angle() only implemented for translation surfaces");
+
   int angle = 0;
 
   const HalfEdge first = *begin(this->atVertex(vertex));
@@ -759,6 +810,9 @@ int FlatTriangulation<T>::angle(const Vertex &vertex) const {
 
 template <typename T>
 std::optional<Deformation<FlatTriangulation<T>>> FlatTriangulation<T>::isomorphism(const FlatTriangulation<T> &other, ISOMORPHISM kind, std::function<bool(const T &, const T &, const T &, const T &)> filterMatrix, std::function<bool(HalfEdge, HalfEdge)> filterHalfEdgeMap) const {
+  if (!translation_surface())
+    throw std::logic_error("isomorphism() only implemented for translation surfaces");
+
   if (this->hasBoundary() != other.hasBoundary())
     return std::nullopt;
 
@@ -955,11 +1009,21 @@ ImplementationOf<FlatTriangulation<T>>::ImplementationOf(FlatTriangulationCombin
 
 template <typename T>
 void ImplementationOf<FlatTriangulation<T>>::updateAfterFlip(OddHalfEdgeMap<Vector<T>> &vectors, const FlatTriangulationCombinatorial &parent, HalfEdge flip) {
+  const auto &surface = reinterpret_cast<const FlatTriangulation<T> &>(parent);
+
+  if (!surface.translation_surface())
+    throw std::logic_error("flipping only implemented for translation surfaces");
+
   vectors.set(flip, vectors.get(-parent.nextInFace(flip)) + vectors.get(-parent.previousInFace(flip)));
 }
 
 template <typename T>
 std::string ImplementationOf<FlatTriangulation<T>>::yaml() const {
+  const auto self = const_cast<ImplementationOf<FlatTriangulation<T>>*>(this)->from_this();
+
+  if (!self.translation_surface())
+    throw std::logic_error("yaml export only implemented for translation surfaces");
+
   std::stringstream out;
   out << std::setprecision(std::numeric_limits<double>::digits10);
   out << "vertices:" << std::endl;
@@ -987,13 +1051,20 @@ std::string ImplementationOf<FlatTriangulation<T>>::yaml() const {
 template <typename T>
 void ImplementationOf<FlatTriangulation<T>>::updateApproximationAfterFlip(OddHalfEdgeMap<flatsurf::Vector<exactreal::Arb>> &vectors, const FlatTriangulationCombinatorial &combinatorial, HalfEdge flip) {
   const auto &surface = reinterpret_cast<const FlatTriangulation<T> &>(combinatorial);
+
+  if (!surface.translation_surface())
+    throw std::logic_error("flipping only implemented for translation surfaces");
+
   vectors.set(flip, static_cast<flatsurf::Vector<exactreal::Arb>>(surface.fromHalfEdge(-surface.nextInFace(flip)) + surface.fromHalfEdge(-surface.previousInFace(flip))));
 }
 
 template <typename T>
 void ImplementationOf<FlatTriangulation<T>>::flip(HalfEdge e) {
   const auto self = from_this();
 
+  if (!self.translation_surface())
+    throw std::logic_error("flip only implemented for translation surfaces");
+
   LIBFLATSURF_CHECK_ARGUMENT(self.convex(e, true), "cannot flip this edge as a resulting face would not be strictly convex");
 
   ImplementationOf<FlatTriangulationCombinatorial>::flip(e);
@@ -1003,6 +1074,9 @@ template <typename T>
 void ImplementationOf<FlatTriangulation<T>>::check() {
   const auto self = from_this();
 
+  if (!self.translation_surface())
+    throw std::logic_error("checking only implemented for translation surfaces");
+
   // check that faces are closed
   for (auto edge : self.halfEdges()) {
     if (self.boundary(edge)) continue;
@@ -1025,6 +1099,9 @@ void ImplementationOf<FlatTriangulation<T>>::check() {
 
 template <typename T>
 std::ostream &operator<<(std::ostream &os, const FlatTriangulation<T> &self) {
+  if (!self.translation_surface())
+    throw std::logic_error("printing only implemented for translation surfaces");
+
   return os << static_cast<const FlatTriangulationCombinatorial &>(self)
             << " with vectors " << self.self->vectors;
 }

@@ -52,7 +52,6 @@
 #include "external/rx-ranges/include/rx/ranges.hpp"
 #include "impl/assert_connection.hpp"
 #include "impl/collapsed_half_edge.hpp"
-#include "impl/flat_triangulation.impl.hpp"
 #include "impl/flat_triangulation_collapsed.impl.hpp"
 #include "impl/flat_triangulation_combinatorial.impl.hpp"
 #include "impl/saddle_connection.impl.hpp"

@@ -29,7 +29,6 @@
 #include "../../flatsurf/edge_map.hpp"
 #include "../../flatsurf/vector.hpp"
 #include "chain_vector.hpp"
-#include "flat_triangulation.impl.hpp"
 #include "read_only.hpp"
 
 namespace flatsurf {

@@ -26,8 +26,6 @@
 #include <unordered_set>
 
 #include "../../flatsurf/interval_exchange_transformation.hpp"
-#include "flat_triangulation.impl.hpp"
-#include "flat_triangulation_collapsed.impl.hpp"
 #include "lengths.hpp"
 #include "read_only.hpp"
 

@@ -24,7 +24,6 @@
 #include "../../flatsurf/flat_triangulation.hpp"
 #include "../../flatsurf/flat_triangulation_collapsed.hpp"
 #include "../../flatsurf/saddle_connection.hpp"
-#include "flat_triangulation.impl.hpp"
 #include "read_only.hpp"
 
 namespace flatsurf {

@@ -25,7 +25,6 @@
 #include "../../flatsurf/bound.hpp"
 #include "../../flatsurf/saddle_connections.hpp"
 #include "../../flatsurf/vector.hpp"
-#include "flat_triangulation.impl.hpp"
 #include "read_only.hpp"
 
 namespace flatsurf {

@@ -24,8 +24,7 @@
 #include <unordered_set>
 
 #include "../../flatsurf/vertical.hpp"
-#include "flat_triangulation.impl.hpp"
-#include "flat_triangulation_collapsed.impl.hpp"
+#include "../../flatsurf/tracked.hpp"
 #include "managed_movable.impl.hpp"
 #include "read_only.hpp"
 

@@ -25,6 +25,7 @@
 #include "../flatsurf/saddle_connections_sample.hpp"
 #include "impl/saddle_connections_sample.impl.hpp"
 #include "impl/saddle_connections_sample_iterator.impl.hpp"
+#include "util/assert.ipp"
 
 namespace flatsurf {
 

@@ -58,7 +58,10 @@ class SurfaceGenerator : public Catch::Generators::IGenerator<std::tuple<std::st
     GOLDEN_L,
     _1072221,
     _12345034015233,
-    LAST = _12345034015233,
+    HALF_TRANSLATION_SPHERE,
+    DILATION_TORUS,
+    HALF_DILATION_SPHERE,
+    LAST = HALF_DILATION_SPHERE,
   };
 
   Surface current;
@@ -207,6 +210,18 @@ class SurfaceGenerator : public Catch::Generators::IGenerator<std::tuple<std::st
         } else
           value = nullptr;
         return;
+      case Surface::HALF_TRANSLATION_SPHERE:
+        name = "HalfTranslationSphere";
+        value = makeHalfTranslationSphere<R2>();
+        return;
+      case Surface::DILATION_TORUS:
+        name = "DilationTorus";
+        value = makeDilationTorus<R2>();
+        return;
+      case Surface::HALF_DILATION_SPHERE:
+        name = "HalfDilationSphere";
+        value = makeHalfDilationSphere<R2>();
+        return;
       default:
         throw std::logic_error("unknown Surface type");
     }