46 if (levelSet ==
nullptr) {
48 .addWarning(
"No level set was passed to ToSurfaceMesh.")
52 if (mesh ==
nullptr) {
54 .addWarning(
"No mesh was passed to ToSurfaceMesh.")
59 if (levelSet->getNumberOfPoints() == 0) {
64 const unsigned int corner0[12] = {0, 1, 2, 0, 4, 5, 6, 4, 0, 1, 3, 2};
65 const unsigned int corner1[12] = {1, 3, 3, 2, 5, 7, 7, 6, 4, 5, 7, 6};
66 const unsigned int direction[12] = {0, 1, 0, 1, 0, 1, 0, 1, 2, 2, 2, 2};
70 if (levelSet->getLevelSetWidth() < 2) {
72 .addWarning(
"Levelset is less than 2 layers wide. Expanding levelset "
78 typedef std::map<viennahrle::Index<D>,
unsigned> nodeContainerType;
80 nodeContainerType nodes[
D];
81 typename nodeContainerType::iterator nodeIt;
82 const bool updateData = updatePointData;
86 std::vector<std::vector<unsigned>> newDataSourceIds;
89 newDataSourceIds.resize(1);
92 for (viennahrle::ConstSparseCellIterator<hrleDomainType> cellIt(
93 levelSet->getDomain());
94 !cellIt.isFinished(); cellIt.next()) {
96 for (
int u = 0; u <
D; u++) {
97 while (!nodes[u].empty() &&
98 nodes[u].begin()->first <
99 viennahrle::Index<D>(cellIt.getIndices()))
100 nodes[u].erase(nodes[u].begin());
104 for (
int i = 0; i < (1 <<
D); i++) {
105 if (cellIt.getCorner(i).getValue() >=
T(0))
112 if (signs == (1 << (1 <<
D)) - 1)
116 const int *Triangles =
120 for (; Triangles[0] != -1; Triangles +=
D) {
121 std::array<unsigned, D> nod_numbers;
124 for (
int n = 0; n <
D; n++) {
125 const int edge = Triangles[n];
127 unsigned p0 = corner0[edge];
128 unsigned p1 = corner1[edge];
131 auto dir = direction[edge];
134 viennahrle::Index<D> d(cellIt.getIndices());
135 d += viennahrle::BitMaskToIndex<D>(p0);
137 nodeIt = nodes[dir].find(d);
138 if (nodeIt != nodes[dir].end()) {
139 nod_numbers[n] = nodeIt->second;
143 Vec3D<T> cc = {0., 0., 0.};
144 std::size_t currentPointId = 0;
145 for (
int z = 0; z <
D; z++) {
150 static_cast<double>(cellIt.getIndices(z) +
151 viennahrle::BitMaskToIndex<D>(p0)[z]);
155 d0 = cellIt.getCorner(p0).getValue();
156 d1 = cellIt.getCorner(p1).getValue();
160 currentPointId = cellIt.getCorner(p0).getPointId();
161 cc[z] =
static_cast<T>(cellIt.getIndices(z)) + 0.5;
163 if (std::abs(d0) <= std::abs(d1)) {
164 currentPointId = cellIt.getCorner(p0).getPointId();
166 static_cast<T>(cellIt.getIndices(z)) + (d0 / (d0 - d1));
168 currentPointId = cellIt.getCorner(p1).getPointId();
169 cc[z] =
static_cast<T>(cellIt.getIndices(z) + 1) -
173 cc[z] = std::max(cc[z], cellIt.getIndices(z) + epsilon);
174 cc[z] = std::min(cc[z], (cellIt.getIndices(z) + 1) - epsilon);
176 cc[z] = levelSet->getGrid().getGridDelta() * cc[z];
181 mesh->insertNextNode(cc);
182 nodes[dir][d] = nod_numbers[n];
185 newDataSourceIds[0].push_back(currentPointId);
189 if (!triangleMisformed(nod_numbers))
190 mesh->insertNextElement(nod_numbers);
196 mesh->getPointData().translateFromMultiData(levelSet->getPointData(),