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