Compiler Explorer

Source code

#include <algorithm>
#include <array>
#include <cstdio>
#include <functional>
#include <numeric>
#include <utility>

#include <strong_type/strong_type.hpp>

#define Debug

#ifdef Debug
#include <iomanip>
#include <iostream>
#endif

namespace tri {

template <class T, T... Args>
constexpr T get(std::integer_sequence<T, Args...>, std::size_t i) noexcept {
 constexpr T arr[] = {Args...};
 return arr[i];
}

template <class Seq>
constexpr auto get(std::size_t i) noexcept {
 return get(Seq{}, i);
}

template <typename F, class T, T... args>
constexpr bool all(F &&f, std::integer_sequence<T, args...>) noexcept {
 return (f(args) && ...);
}
template <size_t N, typename F>
constexpr bool all(F &&f) noexcept {
 return all(std::forward<F>(f), std::make_index_sequence<N>{});
}
template <typename F, class T, T... args>
constexpr bool any(F &&f, std::integer_sequence<T, args...>) noexcept {
 return (f(args) || ...);
}
template <size_t N, typename F>
constexpr bool any(F &&f) noexcept {
 return any(std::forward<F>(f), std::make_index_sequence<N>{});
}

template <typename U, typename T, size_t D>
constexpr auto strong_cast(const std::array<T, D> &arg) noexcept {
 return [&]<size_t... inds>(std::index_sequence<inds...>)->std::array<U, D> {
 return {U{value_of(arg[inds])}...};
 }
 (std::make_index_sequence<D>{});
}

template <typename T, typename U, size_t D>
constexpr auto add(const std::array<T, D> &a,
 const std::array<U, D> &b) noexcept {
 return [&]<size_t... inds>(std::index_sequence<inds...>)
 ->std::array<decltype(std::declval<T>() + std::declval()), D> {
 return {(a[inds] + b[inds])...};
 }
 (std::make_index_sequence<D>{});
}

template <typename T, size_t D>
constexpr auto sum(const std::array<T, D> &a) noexcept {
 return [&]<size_t... inds>(std::index_sequence<inds...>) {
 return (a[inds] + ...);
 }
 (std::make_index_sequence<D>{});
}

template <typename T, size_t D, typename F>
constexpr auto sum(const std::array<T, D> &a, F &&fun) noexcept {
 return [&]<size_t... inds>(std::index_sequence<inds...>) {
 return (fun(a[inds]) + ...);
 }
 (std::make_index_sequence<D>{});
}

template <typename T, size_t D>
constexpr auto product(const std::array<T, D> &a) noexcept {
 return [&]<size_t... inds>(std::index_sequence<inds...>) {
 return (a[inds] * ... * T{1});
 }
 (std::make_index_sequence<D>{});
}
template <typename T>
constexpr T power(T base, T pow) noexcept {
 T res = base;
 if (pow == 0) return T{0};
 for (size_t i = 1; i < pow; ++i) res *= base;
 return res;
};

template <typename T, typename U, size_t D>
constexpr auto categorise(const std::array<T, D> &pos,
 const std::array<U, D> &dim) noexcept {
 return [&]<size_t... inds>(std::index_sequence<inds...>)->std::array<int, D> {
 return {
 (pos[inds] == T{0} ? -1 : (pos[inds] == T{dim[inds] - 1} ? 1 : 0))...};
 }
 (std::make_index_sequence<D>{});
}

template <size_t SelectionSize, typename T, size_t InputSize>
constexpr auto permutations(const std::array<T, InputSize> &values) {
 std::array<std::array<T, SelectionSize>, power(InputSize, SelectionSize)>
 res{};
 std::array<size_t, SelectionSize> index{0};
 for (size_t n = 0; n < res.size(); ++n) {
 for (size_t s = 0; s < SelectionSize; ++s) {
 res[n][s] = values[index[s]];
 }

size_t i = 0;
 while (i < index.size() && ++index[i] >= InputSize) {
 index[i] = 0;
 ++i;
 }
 }
 return res;
}

template <typename PBC, typename GridInd, typename I, size_t D>
constexpr auto pbc(std::array<GridInd, D> pos, const std::array<I, D> &dim) {
 if constexpr (all([](auto x) { return !x; }, PBC{})) {
 return pos;
 } else {
 for (size_t d = 0; d < PBC::size(); ++d) {
 if (get<PBC>(d)) {
 pos[d] =
 value_of(pos[d]) < 0
 ? pos[d] + dim[d]
 : value_of(pos[d]) >= dim[d] ? pos[d] - dim[d] : pos[d];
 }
 }
 return pos;
 }
}

namespace tags {

struct linearCellId;
struct gridCellId;

struct linearVertexId;
struct gridVertexId;

struct linearEdgeId;
struct gridEdgeId;
struct setEdgeId;

struct linearTriangleId;
struct gridTriangleId;

struct linearTetraId;
struct gridTetraId;

struct boundaryTypeId;
struct flipId;

}  // namespace tags

template <typename T, typename Tag>
using index_type =
 strong::type<T, Tag, strong::regular, strong::bicrementable,
 strong::affine_point<T>, strong::ostreamable, strong::ordered>;

template <typename PBC, typename I, bool Flip = false>
struct EdgeSet {
 static constexpr auto D = PBC::size();
 using BID = index_type<I, tags::boundaryTypeId>;
 using FID = index_type<I, tags::flipId>;

static constexpr size_t size() noexcept { return edges.size(); }

static constexpr std::array<I, D> edge(I setId, I id, FID fid) noexcept {
 return flippedEdges[value_of(fid)][setId][id];
 }

static constexpr const auto &getStar(BID bid, FID fid) noexcept {
    return stars[value_of(fid)][value_of(bid)];
  }

static constexpr bool extendsInDir(I setId, I dim) noexcept {
    return edges[setId][0][dim] == edges[setId][1][dim];
  }

static constexpr size_t getBoundaryType(
 const std::array<int, D> &onBoundary) {
 size_t res = 0;
 size_t factor = 1;
 for (int i = 0; i < D; ++i) {
 res += (onBoundary[i] + 1) * factor;
 factor *= D;
 }
 return res;
 }

template <typename T, typename U>
 static constexpr BID getBoundaryType(const std::array<T, D> &pos,
 const std::array<U, D> &dim) noexcept {
 size_t res = 0;
 size_t factor = 1;
 for (int i = 0; i < D; ++i) {
 if (pos[i] == T{0}) {
 res += 0 * factor;
 } else if (pos[i] == T{dim[i] - 1}) {
 res += 2 * factor;
 } else {
 res += 1 * factor;
 }
 factor *= D;
 }
 return BID{res};
 }

static constexpr FID getFlipType(const std::array<I, D> &pos) noexcept {
 if constexpr (Flip) {
 size_t res = 0;
 size_t factor = 1;
 for (int i = 0; i < D; ++i) {
 res += (pos[i] % 2) * factor;
 factor *= 2;
 }
 return FID{res};
 } else {
 return FID{0};
 }
 }

private:
 static constexpr auto makeEdges() noexcept {
 if constexpr (D == 1) {
 return std::array<std::array<std::array<I, D>, 2>, 1>{{{{{0}, {1}}}}};
 } else if constexpr (D == 2) {
 return std::array<std::array<std::array<I, D>, 2>, 3>{
 {{{{0, 0}, {1, 0}}}, {{{0, 0}, {0, 1}}}, {{{0, 0}, {1, 1}}}}};
 } else {
 return std::array<std::array<std::array<I, D>, 2>, 7>{{
 {{{0, 0, 0}, {1, 0, 0}}}, // L
 {{{0, 0, 0}, {0, 1, 0}}}, // H
 {{{0, 0, 0}, {0, 0, 1}}}, // P
 {{{1, 0, 0}, {0, 1, 0}}}, // D1
 {{{0, 0, 0}, {0, 1, 1}}}, // D2
 {{{1, 0, 0}, {0, 0, 1}}}, // D3
 {{{1, 0, 0}, {0, 1, 1}}} // D4
 }};
 }
 }

static constexpr auto makeBoundaryTypes() noexcept {
 return permutations<D>(std::array<I, 3>{{-1, 0, 1}});
 }

static constexpr auto makeFlipTypes() noexcept {
 if constexpr (Flip) {
 return permutations<D>(std::array<bool, 2>{{false, true}});
 } else {
 return std::array<std::array<bool, D>, 1>{{false, false, false}};
 }
 }

static constexpr std::array<I, D> flip(size_t setId,
 const std::array<bool, D> &flipType,
 std::array<I, D> edge) noexcept {
 const auto mirror = [](auto &item, size_t i) {
 item[i] = item[i] == 0 ? 1 : 0;
 };

if constexpr (D == 2 && Flip) {
      if (setId == 2 && flipType[0]) mirror(edge, 0);
      if (setId == 2 && flipType[1]) mirror(edge, 1);
    } else if constexpr (D == 3 && Flip) {
      if (setId == 3 && flipType[0]) mirror(edge, 0);
      if (setId == 3 && flipType[1]) mirror(edge, 1);

if (setId == 4 && flipType[1]) mirror(edge, 1);
      if (setId == 4 && flipType[2]) mirror(edge, 2);

if (setId == 5 && flipType[0]) mirror(edge, 0);
      if (setId == 5 && flipType[2]) mirror(edge, 2);

if (setId == 6 && flipType[0]) mirror(edge, 0);
      if (setId == 6 && flipType[1]) mirror(edge, 1);
      if (setId == 6 && flipType[2]) mirror(edge, 2);
    }
    return edge;
  }

static constexpr std::array<I, D> flip(
 size_t setId, FID fid, const std::array<I, D> &edge) noexcept {
 return flip(setId, flipTypes[value_of(fid)], edge);
 }

static constexpr auto makeFlippedEdges() noexcept {
 std::array<std::array<std::array<std::array<I, D>, 2>, edges.size()>,
 flipTypes.size()>
 flippedEdges{};

for (size_t f = 0; f < flipTypes.size(); ++f) {
 for (size_t e = 0; e < edges.size(); ++e) {
 flippedEdges[f][e][0] = flip(e, FID{f}, edges[e][0]);
 flippedEdges[f][e][1] = flip(e, FID{f}, edges[e][1]);
 }
 }
 return flippedEdges;
 }

// Needs to consider Flipping
 static constexpr auto makeStars() noexcept {
 constexpr auto offsets = permutations<D>(std::array<I, 2>{{-1, 0}});
 constexpr auto starMaxSize = power(size_t{3}, D) - 1;

std::array<
 std::array<std::pair<size_t, std::array<std::array<I, D>, starMaxSize>>,
 boundaryTypes.size()>,
 flipTypes.size()>
 stars{};

auto getFlippedId = [&](size_t fid, const std::array<I, D> &offset) {
 const auto &ftype = flipTypes[fid];
 size_t res = 0;
 size_t factor = 1;
 for (size_t i = 0; i < D; ++i) {
 if (offset[i] != 0) {
 res += size_t{!ftype[i]} * factor;
 } else {
 res += size_t{ftype[i]} * factor;
 }
 factor *= 2;
 }
 return res;
 };

const auto inStar = [](const std::array<I, D> &offset,
 const std::array<I, D> &edge) -> bool {
 return all<D>([&](size_t d) { return offset[d] + edge[d] == 0; });
 };

const auto boundayFilter = [](size_t b,
 const std::array<I, D> &edge) -> bool {
 return any<D>([&](size_t d) {
 return (!get<PBC>(d) && boundaryTypes[b][d] < 0 && edge[d] < 0) ||
 (!get<PBC>(d) && boundaryTypes[b][d] > 0 && edge[d] > 0);
 });
 };

for (size_t f = 0; f < flipTypes.size(); ++f) {
 for (size_t b = 0; b < boundaryTypes.size(); ++b) {
 size_t count = 0;
 const auto addEdge = [&](const std::array<I, D> &edge) {
 if (count < stars[f][b].second.size()) {
 stars[f][b].second[count] = edge;
 ++count;
 }
 };

for (const auto &offset : offsets) {
          const auto &fedges = flippedEdges[getFlippedId(f, offset)];
          for (const auto &edge : fedges) {
            if (inStar(offset, edge[0])) {
              auto vert = add(offset, edge[1]);
              if (!boundayFilter(b, vert)) addEdge(vert);
            } else if (inStar(offset, edge[1])) {
              auto vert = add(offset, edge[0]);
              if (!boundayFilter(b, vert)) addEdge(vert);
            }
          }
        }
        stars[f][b].first = count;
      }
    }

return stars;
  }

public:
  static constexpr decltype(makeEdges()) edges = makeEdges();
  static constexpr decltype(makeBoundaryTypes()) boundaryTypes =
      makeBoundaryTypes();
  static constexpr decltype(makeFlipTypes()) flipTypes = makeFlipTypes();

static constexpr decltype(makeFlippedEdges()) flippedEdges =
      makeFlippedEdges();

static constexpr decltype(makeStars()) stars = makeStars();
};

template <size_t N, typename LinearInd, typename GridInd, typename PBC>
struct IndexMapper {};

template <typename LinearInd, typename GridInd, typename PBC>
struct IndexMapper<1, LinearInd, GridInd, PBC> {
 static_assert(std::is_same_v<strong::underlying_type_t<LinearInd>,
 strong::underlying_type_t<GridInd>>);

using I = strong::underlying_type_t<GridInd>;
 using Grid = std::array<GridInd, 1>;

constexpr IndexMapper(const std::array<I, 1> &dim, const I offset = 0)
 : dim{dim}, offset{offset} {};

constexpr LinearInd operator()(const Grid &pos) const noexcept {
 const auto p = pbc<PBC>(pos, dim);
 return offset + LinearInd{value_of(p[0])};
 }
 constexpr Grid operator()(const LinearInd index) const noexcept {
 index -= offset;
 return {GridInd{value_of(index)}};
 }
 const std::array<I, 1> dim;
 const I offset;
};

template <typename LinearInd, typename GridInd, typename PBC>
struct IndexMapper<2, LinearInd, GridInd, PBC> {
 static_assert(std::is_same_v<strong::underlying_type_t<LinearInd>,
 strong::underlying_type_t<GridInd>>);

using I = strong::underlying_type_t<GridInd>;
 using Grid = std::array<GridInd, 2>;

constexpr IndexMapper(const std::array<I, 2> &dim, I offset = 0)
 : dim{dim}, offset{offset} {};

constexpr LinearInd operator()(const Grid &pos) const noexcept {
 const auto p = pbc<PBC>(pos, dim);
 return offset +
 LinearInd{value_of(p[0]) + value_of(p[1]) * dim[0]};
 }
 constexpr Grid operator()(LinearInd index) const noexcept {
 index -= offset;
 return {GridInd{value_of(index) % dim[0]},
 GridInd{value_of(index) / dim[0]}};
 }
 const std::array<I, 2> dim;
 const I offset;
};

template <typename LinearInd, typename GridInd, typename PBC>
struct IndexMapper<3, LinearInd, GridInd, PBC> {
 static_assert(std::is_same_v<strong::underlying_type_t<LinearInd>,
 strong::underlying_type_t<GridInd>>);

using I = strong::underlying_type_t<GridInd>;
 using Grid = std::array<GridInd, 3>;

constexpr IndexMapper(const std::array<I, 3> &dim,
 const I offset = 0) noexcept
 : dim{dim}, offset{offset}, dimxy{dim[0] * dim[1]} {};

constexpr LinearInd operator()(const Grid &pos) const noexcept {
 const auto p = pbc<PBC>(pos, dim);
 return offset +
 LinearInd{value_of(p[0]) + value_of(p[1]) * dim[0] +
 value_of(p[2]) * dimxy};
 }
 constexpr Grid operator()(LinearInd index) const noexcept {
 index -= offset;
 return {GridInd{value_of(index) % dim[0]},
 GridInd{(value_of(index) % dimxy) / dim[0]},
 GridInd{value_of(index) / dimxy}};
 }

const std::array<I, 3> dim;
 const I offset;

private:
  const I dimxy;
};  // namespace tri

template <size_t N, size_t Size, typename LinearInd, typename GridInd,
 typename SetInd, typename PBC>
struct IndexMapperSet {
 using IM = IndexMapper<N, LinearInd, GridInd, PBC>;
 using I = typename IM::I;
 using Grid = typename IM::Grid;

constexpr IndexMapperSet(
 std::array<std::array<I, N>, Size> dimensions) noexcept
 : indexMappers{build(std::make_index_sequence<dimensions.size()>{},
 dimensions)},
 size{count(std::make_index_sequence<dimensions.size()>{}, dimensions)} {
 }

constexpr LinearInd operator()(const Grid &pos, SetInd k) const noexcept {
 return indexMappers[k](pos);
 }
 constexpr std::pair<Grid, SetInd> operator()(LinearInd index) const noexcept {
 size_t k = Size - 1;
 for (; k >= 0; --k) {
 if (value_of(index) >= indexMappers[k].offset) break;
 }
 return {indexMappers[k](index), SetInd{k}};
 }

constexpr const IM &operator[](SetInd i) const noexcept {
    return indexMappers[value_of(i)];
  }

const std::array<IM, Size> indexMappers;
 const I size;

private:
 template <size_t... Ind>
 constexpr auto build(std::index_sequence<Ind...>,
 std::array<std::array<I, N>, Size> dimensions)
 -> std::array<IM, Size> {
 auto create = [size = I{0}](auto dim) mutable {
 const auto offset = size;
 size += product(dim);
 return IM{dim, offset};
 };
 return {{create(dimensions[Ind])...}};
 }

template <size_t... Ind>
 constexpr auto count(std::index_sequence<Ind...>,
 std::array<std::array<I, N>, Size> dim) {
 return (product(dim[Ind]) + ...);
 }
};

template <typename T, typename I, typename PBC, bool Flip>
class ImplicitTriangulation final {
 public:
 static constexpr auto D = PBC::size();
 static_assert(D < 4);

using Edges = EdgeSet<PBC, I, Flip>;

using CID = index_type<I, tags::linearCellId>;
 using CGID = index_type<I, tags::gridCellId>;

using VID = index_type<I, tags::linearVertexId>;
 using VGID = index_type<I, tags::gridVertexId>;

using EID = index_type<I, tags::linearEdgeId>;
 using EGID = index_type<I, tags::gridEdgeId>;
 using ESID = index_type<I, tags::setEdgeId>;

using TID = index_type<I, tags::linearTriangleId>;
 using TGID = index_type<I, tags::gridTriangleId>;

using RID = index_type<I, tags::linearTetraId>;
 using RGID = index_type<I, tags::gridTetraId>;

using VertexMapper = IndexMapper<D, VID, VGID, PBC>;
 using EdgeMapper = IndexMapperSet<D, Edges::size(), EID, EGID, ESID, PBC>;
 using TriangleMapper = IndexMapper<D, TID, TGID, PBC>;

explicit constexpr ImplicitTriangulation(
 const std::array<T, D> &origin, const std::array<T, D> &spacing,
 const std::array<I, D> &dimensions) noexcept
 : origin{origin},
 spacing{spacing},
 dimensions{dimensions},
 vertexMapper{dimensions},
 edgeMapper{createEdgesSetDimensions(dimensions)} {}

constexpr I getDimensionality() const noexcept { return D; };

VID getNumberOfVertices() const noexcept {
 return VID{std::accumulate(dimensions.begin(), dimensions.end(), 1,
 std::multiplies<>{})};
 }

constexpr auto vertexToPosition(const VID vertexId) const noexcept {
    return vertexMapper(vertexId);
  }

constexpr std::array<T, D> getVertexPoint(const VID vertexId) const noexcept {
 const auto pos = vertexMapper(vertexId);
 std::array<T, D> res{};
 for (size_t d = 0; d < D; ++d) {
 res[d] = origin[d] + value_of(pos[d]) * spacing[d];
 }
 return res;
 }

constexpr bool isVertexOnBoundary(const VID vertexId) const noexcept {
 const auto id = vertexMapper(vertexId);
 for (I i = 0; i < D; ++i) {
 if (!get<PBC>(i) &&
 (id[i] == VGID{0} || id[i] == VGID{dimensions[i] - 1})) {
 return true;
 }
 }
 return false;
 }

constexpr VID getVertexNeighborNumber(const VID vertexId) const noexcept {
 const auto pos = vertexMapper(vertexId);
 const auto type = Edges::getBoundaryType(pos, dimensions);
 const auto flip = Edges::getFlipType(strong_cast(pos));
 return VID{Edges::getStar(type, flip).first};
 }

constexpr VID getVertexNeighbor(const VID vertexId,
 const I localNeighborId) const noexcept {
 const auto pos = vertexMapper(vertexId);
 const auto type = Edges::getBoundaryType(pos, dimensions);
 const auto flip = Edges::getFlipType(strong_cast(pos));
 const auto edge = Edges::getStar(type, flip).second[localNeighborId];
 return vertexMapper(add(pos, edge));
 }

std::vector<VID> getVertexNeighbors(const VID vertexId) const noexcept {
 const auto pos = vertexMapper(vertexId);
 const auto type = Edges::getBoundaryType(pos, dimensions);
 const auto flip = Edges::getFlipType(strong_cast(pos));
 const auto &star = Edges::getStar(type, flip);

std::vector<VID> res(star.first);
 std::transform(
 star.second.begin(), star.second.begin() + star.first, res.begin(),
 [&](const auto &edge) { return vertexMapper(add(pos, edge)); });
 return res;
 }

constexpr EID getVertexEdgeNumber(const VID vertexId) const noexcept {
 const auto pos = vertexMapper(vertexId);
 const auto type = Edges::getBoundaryType(pos, dimensions);
 const auto flip = Edges::getFlipType(strong_cast(pos));
 return EID{Edges::getStar(type, flip).first};
 }

constexpr EID getVertexEdge(const VID vertexId, const EID id) const noexcept {
 const auto pos = vertexMapper(vertexId);
 const auto type = Edges::getBoundaryType(pos, dimensions);
 const auto flip = Edges::getFlipType(strong_cast(pos));
 return edgeMapper(string_cast<EGID>(pos), value_of(id));
 }

std::vector<std::vector<EID>> getVertexEdges() const {
 std::vector<std::vector<EID>> res;
 for (VID v = 0; v < getNumberOfVertices(); ++v) {
 res.emplace_back();
 for (EID e = 0; e < getVertexEdgeNumber(v); ++e) {
 res.back().push_back(getVertexEdge(v, e));
 }
 }
 return res;
 }

constexpr I getCellVertexNumber(const CID cellId) const noexcept {
    return D - 1;
  }

EID getNumberOfEdges() const { return EID{edgeMapper.size}; }

constexpr auto edgeToPosition(const EID edgeId) const noexcept {
    return edgeMapper(edgeId);
  }

constexpr VID getEdgeVertex(const EID edgeId, const I id) const {
 const auto [epos, set] = edgeMapper(edgeId);
 const auto fid = Edges::getFlipType(strong_cast(epos));
 const auto vpos = strong_cast<VGID>(
 add(epos, Edges::edge(value_of(set), id, fid)));
 return vertexMapper(vpos);
 }

constexpr bool isEdgeOnBoundary(const EID edgeId) const noexcept {
 const auto [epos, set] = edgeMapper(edgeId);
 const auto onBoundary = categorise(epos, edgeMapper[set].dim);
 return any<onBoundary.size()>([&, s = value_of(set)](auto i) {
 return !get<PBC>(i) && onBoundary[i] != 0 && Edges::extendsInDir(s, i);
 });
 }

private:
 static constexpr auto createEdgesSetDimensions(
 const std::array<I, D> &dim) noexcept {
 std::array<std::array<I, D>, Edges::size()> res{};
 for (size_t e = 0; e < Edges::size(); ++e) {
 for (size_t d = 0; d < D; ++d) {
 if (!get<PBC>(d) && !Edges::extendsInDir(e, d)) {
 res[e][d] = dim[d] - 1;
 } else {
 res[e][d] = dim[d];
 }
 }
 }
 return res;
 }

std::array<T, D> origin;
 std::array<T, D> spacing;
 std::array<I, D> dimensions;
 VertexMapper vertexMapper;
 EdgeMapper edgeMapper;
};

}  // namespace tri

#ifdef Debug

template <typename T, size_t D>
std::ostream &print(std::ostream &os, const std::array<T, D> &data) {
 os << "{";
 bool first = true;
 for (const auto &d : data) {
 if (!first) os << ",";
 first = false;
 os << d;
 }
 os << "}";
 return os;
}

template <typename T>
struct printer {
 printer(const T &data) : data{data} {}
 const T &data;
};

template <typename T>
std::ostream &operator<<(std::ostream &os, const printer<T> &data) {
 print(os, data.data);
 return os;
}

template <typename Trig>
void debug(const Trig &trig) {
 using VID = typename Trig::VID;
 using EID = typename Trig::EID;

std::cout << "Nvert: " << trig.getNumberOfVertices() << "\n";
 std::cout << "Nedge: " << trig.getNumberOfEdges() << "\n";

{
 std::cout << "vertices = {";
 bool first = true;
 for (VID i{0}; i < trig.getNumberOfVertices(); ++i) {
 if (!first) std::cout << ", ";
 first = false;

std::cout << "Vertex[{" << printer{trig.getVertexPoint(i)} << "}"
 << ", VID->" <True" : "")
 << ", Neighbours->" << trig.getVertexNeighborNumber(i)
 << ", Edges->" << trig.getVertexEdgeNumber(i) << "]";
 }
 std::cout << "};";
 }

{
 std::cout << "edges = {";
 bool first = true;
 for (EID i{0}; i < trig.getNumberOfEdges(); ++i) {
 if (!first) std::cout << ", ";
 first = false;

auto vid1 = trig.getEdgeVertex(i, 0);
      auto vid2 = trig.getEdgeVertex(i, 1);

std::cout << "Edge[{" << printer{trig.getVertexPoint(vid1)} << "},"
 << " {" << printer{trig.getVertexPoint(vid2)} << "}"
 << (trig.isEdgeOnBoundary(i) ? ", Boundary->True" : "")
 << ", EID->" << i << "]";
 }
 std::cout << "};";
 }
}
#endif

#if 0
extern int a;
extern int b;
extern int c;
#else
constexpr int a = 3;
constexpr int b = 3;
constexpr int c = 3;
#endif

int main() {
 using pbc = std::integer_sequence<bool, true, false, true>;
 using Trig = tri::ImplicitTriangulation<double, int, pbc, true>;
 using VID = typename Trig::VID;
 using EID = typename Trig::EID;
 const auto trig = Trig{{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}, {a, b, c}};

for (VID i{0}; i < trig.getNumberOfVertices(); ++i) {
 std::putchar(trig.isVertexOnBoundary(i) ? 'b' : '*');
 }
 std::putchar('\n');
 for (EID i{0}; i < trig.getNumberOfEdges(); ++i) {
 std::putchar(trig.isEdgeOnBoundary(i) ? 'b' : '*');
 if (value_of(i + 1) % (7 * a) == 0) std::putchar('\n');
 }

#ifdef Debug

std::putchar('\n');
 constexpr auto star = Trig::Edges::stars[3][14];
 std::cout << "Star size: " << star.first << '\n';

std::cout << "star = {";
 bool first = true;
 for (const auto &edge : star.second) {
 if (!first) std::cout << ',';
 first = false;
 std::cout << '{' << edge[0] << ',' << edge[1] << ',' << edge[2] << '}';
 }
 std::cout << "}\n";

std::cout << "Start sizes:\n";
 for (const auto &starset : Trig::Edges::stars) {
 for (const auto &star : starset) {
 std::cout << std::setw(2) << star.first << ' ';
 }
 std::cout << '\n';
 }
 std::cout << '\n';

const auto i = VID{19};
 {
 const auto pos = trig.getVertexPoint(i);
 std::cout << "node = {" << pos[0] << ',' << pos[1] << ',' << pos[2]
 << "};\n";
 }
 std::cout << "star = {";
 first = true;
 for (const auto vid : trig.getVertexNeighbors(i)) {
 if (!first) std::cout << ',';
 first = false;
 const auto pos = trig.getVertexPoint(vid);
 std::cout << '{' << pos[0] << ',' << pos[1] << ',' << pos[2] << '}';
 }
 std::cout << "};\n";

debug(trig);
#endif

return 0;
}