ViennaLS
Loading...
Searching...
No Matches
lsToSurfaceMesh.hpp
Go to the documentation of this file.
1#pragma once
2
4
5#include <iostream>
6#include <map>
7
8#include <hrleSparseCellIterator.hpp>
9#include <lsDomain.hpp>
10#include <lsExpand.hpp>
11#include <lsMarchingCubes.hpp>
12#include <lsMesh.hpp>
13
14namespace viennals {
15
16using namespace viennacore;
17
21template <class T, int D> class ToSurfaceMesh {
22 typedef typename Domain<T, D>::DomainType hrleDomainType;
23
24 SmartPointer<Domain<T, D>> levelSet = nullptr;
25 SmartPointer<Mesh<T>> mesh = nullptr;
26 // std::vector<hrleIndexType> meshNodeToPointIdMapping;
27 const T epsilon;
28 bool updatePointData = true;
29
30public:
31 explicit ToSurfaceMesh(double eps = 1e-12) : epsilon(eps) {}
32
33 ToSurfaceMesh(const SmartPointer<Domain<T, D>> passedLevelSet,
34 SmartPointer<Mesh<T>> passedMesh, double eps = 1e-12)
35 : levelSet(passedLevelSet), mesh(passedMesh), epsilon(eps) {}
36
37 void setLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
38 levelSet = passedlsDomain;
39 }
40
41 void setMesh(SmartPointer<Mesh<T>> passedMesh) { mesh = passedMesh; }
42
43 void setUpdatePointData(bool update) { updatePointData = update; }
44
45 void apply() {
46 if (levelSet == nullptr) {
47 Logger::getInstance()
48 .addWarning("No level set was passed to ToSurfaceMesh.")
49 .print();
50 return;
51 }
52 if (mesh == nullptr) {
53 Logger::getInstance()
54 .addWarning("No mesh was passed to ToSurfaceMesh.")
55 .print();
56 return;
57 }
58
59 if (levelSet->getNumberOfPoints() == 0) {
60 return;
61 }
62
63 mesh->clear();
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};
67
68 // test if level set function consists of at least 2 layers of
69 // defined grid points
70 if (levelSet->getLevelSetWidth() < 2) {
71 Logger::getInstance()
72 .addWarning("Levelset is less than 2 layers wide. Expanding levelset "
73 "to 2 layers.")
74 .print();
75 Expand<T, D>(levelSet, 2).apply();
76 }
77
78 typedef std::map<viennahrle::Index<D>, unsigned> nodeContainerType;
79
80 nodeContainerType nodes[D];
81 typename nodeContainerType::iterator nodeIt;
82 const bool updateData = updatePointData;
83
84 // save how data should be transferred to new level set
85 // list of indices into the old pointData vector
86 std::vector<std::vector<unsigned>> newDataSourceIds;
87 // there is no multithreading here, so just use 1
88 if (updateData)
89 newDataSourceIds.resize(1);
90
91 // iterate over all active points
92 for (viennahrle::ConstSparseCellIterator<hrleDomainType> cellIt(
93 levelSet->getDomain());
94 !cellIt.isFinished(); cellIt.next()) {
95
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());
101 }
102
103 unsigned signs = 0;
104 for (int i = 0; i < (1 << D); i++) {
105 if (cellIt.getCorner(i).getValue() >= T(0))
106 signs |= (1 << i);
107 }
108
109 // all corners have the same sign, so no surface here
110 if (signs == 0)
111 continue;
112 if (signs == (1 << (1 << D)) - 1)
113 continue;
114
115 // for each element
116 const int *Triangles =
119
120 for (; Triangles[0] != -1; Triangles += D) {
121 std::array<unsigned, D> nod_numbers;
122
123 // for each node
124 for (int n = 0; n < D; n++) {
125 const int edge = Triangles[n];
126
127 unsigned p0 = corner0[edge];
128 unsigned p1 = corner1[edge];
129
130 // determine direction of edge
131 auto dir = direction[edge];
132
133 // look for existing surface node
134 viennahrle::Index<D> d(cellIt.getIndices());
135 d += viennahrle::BitMaskToIndex<D>(p0);
136
137 nodeIt = nodes[dir].find(d);
138 if (nodeIt != nodes[dir].end()) {
139 nod_numbers[n] = nodeIt->second;
140 } else { // if node does not exist yet
141
142 // calculate coordinate of new node
143 Vec3D<T> cc = {0., 0., 0.}; // initialise with zeros
144 std::size_t currentPointId = 0;
145 for (int z = 0; z < D; z++) {
146 if (z != dir) {
147 // TODO might not need BitMaskToVector here, just check if z bit
148 // is set
149 cc[z] =
150 static_cast<double>(cellIt.getIndices(z) +
151 viennahrle::BitMaskToIndex<D>(p0)[z]);
152 } else {
153 T d0, d1;
154
155 d0 = cellIt.getCorner(p0).getValue();
156 d1 = cellIt.getCorner(p1).getValue();
157
158 // calculate the surface-grid intersection point
159 if (d0 == -d1) { // includes case where d0=d1=0
160 currentPointId = cellIt.getCorner(p0).getPointId();
161 cc[z] = static_cast<T>(cellIt.getIndices(z)) + 0.5;
162 } else {
163 if (std::abs(d0) <= std::abs(d1)) {
164 currentPointId = cellIt.getCorner(p0).getPointId();
165 cc[z] =
166 static_cast<T>(cellIt.getIndices(z)) + (d0 / (d0 - d1));
167 } else {
168 currentPointId = cellIt.getCorner(p1).getPointId();
169 cc[z] = static_cast<T>(cellIt.getIndices(z) + 1) -
170 (d1 / (d1 - d0));
171 }
172 }
173 cc[z] = std::max(cc[z], cellIt.getIndices(z) + epsilon);
174 cc[z] = std::min(cc[z], (cellIt.getIndices(z) + 1) - epsilon);
175 }
176 cc[z] = levelSet->getGrid().getGridDelta() * cc[z];
177 }
178
179 // insert new node
180 nod_numbers[n] =
181 mesh->insertNextNode(cc); // insert new surface node
182 nodes[dir][d] = nod_numbers[n];
183
184 if (updateData)
185 newDataSourceIds[0].push_back(currentPointId);
186 }
187 }
188
189 if (!triangleMisformed(nod_numbers))
190 mesh->insertNextElement(nod_numbers); // insert new surface element
191 }
192 }
193
194 // now copy old data into new level set
195 if (updateData) {
196 mesh->getPointData().translateFromMultiData(levelSet->getPointData(),
197 newDataSourceIds);
198 }
199 }
200
201private:
202 static bool triangleMisformed(const std::array<unsigned, D> &nodeNumbers) {
203 if constexpr (D == 3) {
204 return nodeNumbers[0] == nodeNumbers[1] ||
205 nodeNumbers[0] == nodeNumbers[2] ||
206 nodeNumbers[1] == nodeNumbers[2];
207 } else {
208 return nodeNumbers[0] == nodeNumbers[1];
209 }
210 }
211};
212
213// add all template specialisations for this class
215
216} // namespace viennals
static const int * polygonize2d(unsigned int signs)
signs = signs of the corners in lexicographic order (1bit per corner)
Definition lsMarchingCubes.hpp:312
static const int * polygonize3d(unsigned int signs)
signs = signs of the corners in lexicographic order (1bit per corner)
Definition lsMarchingCubes.hpp:324
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
viennahrle::Domain< T, D > DomainType
Definition lsDomain.hpp:32
Expands the levelSet to the specified number of layers. The largest value in the levelset is thus wid...
Definition lsExpand.hpp:17
void apply()
Apply the expansion to the specified width.
Definition lsExpand.hpp:44
This class holds an explicit mesh, which is always given in 3 dimensions. If it describes a 2D mesh,...
Definition lsMesh.hpp:22
Extract an explicit Mesh<> instance from an lsDomain. The interface is then described by explicit sur...
Definition lsToSurfaceMesh.hpp:21
void setMesh(SmartPointer< Mesh< T > > passedMesh)
Definition lsToSurfaceMesh.hpp:41
ToSurfaceMesh(const SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< Mesh< T > > passedMesh, double eps=1e-12)
Definition lsToSurfaceMesh.hpp:33
void apply()
Definition lsToSurfaceMesh.hpp:45
void setUpdatePointData(bool update)
Definition lsToSurfaceMesh.hpp:43
void setLevelSet(SmartPointer< Domain< T, D > > passedlsDomain)
Definition lsToSurfaceMesh.hpp:37
ToSurfaceMesh(double eps=1e-12)
Definition lsToSurfaceMesh.hpp:31
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
Definition lsAdvect.hpp:36
constexpr int D
Definition pyWrap.cpp:71
double T
Definition pyWrap.cpp:69