ViennaLS
Loading...
Searching...
No Matches
lsGeometricAdvect.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <hrleSparseIterator.hpp>
4
6#include <lsConcepts.hpp>
7#include <lsDomain.hpp>
8#include <lsExpand.hpp>
9#include <lsFromMesh.hpp>
12#include <lsToDiskMesh.hpp>
13
14#include <vcLogger.hpp>
15#include <vcSmartPointer.hpp>
16#include <vcVectorType.hpp>
17
18#ifndef NDEBUG // if in debug build
19#include <lsCheck.hpp>
20#include <lsToMesh.hpp>
21#include <lsVTKWriter.hpp>
22#endif
23
24namespace viennals {
25
26using namespace viennacore;
27
35template <class T, int D> class GeometricAdvect {
36 using hrleIndexType = viennahrle::IndexType;
37 using hrleCoordType = viennahrle::CoordType;
38
39 SmartPointer<Domain<T, D>> levelSet = nullptr;
40 SmartPointer<Domain<T, D>> maskLevelSet = nullptr;
41 SmartPointer<const GeometricAdvectDistribution<hrleCoordType, D>> dist =
42 nullptr;
43 static constexpr T cutoffValue =
44 T(1.) + std::numeric_limits<T>::epsilon() * T(100);
45
46 static void incrementIndices(viennahrle::Index<D> &indices,
47 const viennahrle::Index<D> &min,
48 const viennahrle::Index<D> &max) {
49 int dim = 0;
50 for (; dim < D - 1; ++dim) {
51 if (indices[dim] < max[dim])
52 break;
53 indices[dim] = min[dim];
54 }
55 ++indices[dim];
56 }
57
58 template <class K, class V, template <class...> class MapType, class... Ts>
59 MapType<V, K> inverseTranslator(MapType<K, V, Ts...> &map) {
60 MapType<V, K> inv;
61 std::for_each(map.begin(), map.end(), [&inv](const std::pair<K, V> &p) {
62 inv.insert(std::make_pair(p.second, p.first));
63 });
64 return inv;
65 }
66
67public:
68 GeometricAdvect() = default;
69
70 template <class DistType,
72 DistType> = lsConcepts::assignable>
73 GeometricAdvect(SmartPointer<Domain<T, D>> passedLevelSet,
74 SmartPointer<DistType> passedDist,
75 SmartPointer<Domain<T, D>> passedMaskLevelSet = nullptr)
76 : levelSet(passedLevelSet), maskLevelSet(passedMaskLevelSet) {
77 dist = std::dynamic_pointer_cast<
79 }
80
82 void setLevelSet(SmartPointer<Domain<T, D>> passedLevelSet) {
83 levelSet = passedLevelSet;
84 }
85
88 template <class DistType,
90 DistType> = lsConcepts::assignable>
91 void setAdvectionDistribution(SmartPointer<DistType> passedDist) {
92 dist = std::dynamic_pointer_cast<
94 }
95
99 void setMaskLevelSet(SmartPointer<Domain<T, D>> passedMaskLevelSet) {
100 maskLevelSet = passedMaskLevelSet;
101 }
102
104 void apply() {
105 if (levelSet == nullptr) {
106 Logger::getInstance()
107 .addWarning("No level set passed to GeometricAdvect. Not Advecting.")
108 .print();
109 return;
110 }
111 if (dist == nullptr) {
112 Logger::getInstance()
113 .addWarning("No GeometricAdvectDistribution passed to "
114 "GeometricAdvect. Not "
115 "Advecting.")
116 .print();
117 return;
118 }
119
120 // levelSet must have at least a width of 3
121 Expand<T, D>(levelSet, 3).apply();
122
123 if (maskLevelSet != nullptr) {
124 Expand<T, D>(maskLevelSet, 3).apply();
125 }
126
127 typedef typename Domain<T, D>::DomainType DomainType;
128
129 auto &domain = levelSet->getDomain();
130
131 auto &grid = levelSet->getGrid();
132 auto gridDelta = grid.getGridDelta();
133
134 // Extract the original surface as a point cloud of grid
135 // points shifted to the surface (disk mesh)
136 auto surfaceMesh = SmartPointer<Mesh<hrleCoordType>>::New();
137 auto pointIdTranslator =
138 SmartPointer<typename ToDiskMesh<T, D>::TranslatorType>::New();
139 ToDiskMesh<T, D, hrleCoordType>(levelSet, surfaceMesh, pointIdTranslator)
140 .apply();
141 *pointIdTranslator = inverseTranslator(*pointIdTranslator);
142
143 // find bounds of distribution
144 auto distBounds = dist->getBounds();
145
146 // TODO: need to add support for periodic boundary conditions!
147 viennahrle::Index<D> distMin, distMax;
148
149 bool minPointNegative = domain.getDomainSegment(0).definedValues[0] < 0.;
150 bool maxPointNegative =
151 domain.getDomainSegment(domain.getNumberOfSegments() - 1)
152 .definedValues.back() < 0.;
153 bool distIsPositive = true;
154
155 // find bounding box of old domain
156 hrleIndexType bounds[6];
157 domain.getDomainBounds(bounds);
158 viennahrle::Index<D> min, max;
159 for (unsigned i = 0; i < D; ++i) {
160 // translate from coords to indices
161 distMin[i] =
162 distBounds[2 * i] / gridDelta + ((distBounds[2 * i] < 0) ? -2 : 2);
163 distMax[i] = distBounds[2 * i + 1] / gridDelta +
164 ((distBounds[2 * i + 1] < 0) ? -2 : 2);
165 if (distBounds[2 * i] >= 0) {
166 distIsPositive = false;
167 }
168
169 // use the extent of the diskMesh to identify bounding box of new
170 // level set
171 // TODO: respect periodic boundary condition
172 min[i] = surfaceMesh->minimumExtent[i] / gridDelta;
173 // TODO also do the same thing for positive point and etching
174 if (grid.isNegBoundaryInfinite(i) && minPointNegative && distMin[i] < 0) {
175 min[i] -= 2;
176 } else {
177 if (distIsPositive) {
178 min[i] += distMin[i];
179 } else {
180 min[i] -= distMin[i];
181 }
182 }
183 // if calculated index is out of bounds, set the extent
184 // TODO: need to add periodic BNC handling here
185 if (min[i] < grid.getMinGridPoint(i)) {
186 min[i] = grid.getMinGridPoint(i);
187 }
188
189 max[i] = surfaceMesh->maximumExtent[i] / gridDelta;
190 if (grid.isPosBoundaryInfinite(i) && maxPointNegative && distMax[i] > 0) {
191 max[i] += 2;
192 } else {
193 if (distIsPositive) {
194 max[i] += distMax[i];
195 } else {
196 max[i] -= distMax[i];
197 }
198 }
199 if (max[i] > grid.getMaxGridPoint(i)) {
200 max[i] = grid.getMaxGridPoint(i);
201 }
202 }
203
204 // Remove contribute points if they are part of the mask
205 // If a mask is supplied, remove all contribute points which
206 // lie on (or inside) the mask
207 if (maskLevelSet != nullptr) {
208 // Go over all contribute points and see if they are on the mask surface
209 auto &maskDomain = maskLevelSet->getDomain();
210 auto values = surfaceMesh->cellData.getScalarData("LSValues");
211 auto valueIt = values->begin();
212
213 auto newSurfaceMesh = SmartPointer<Mesh<hrleCoordType>>::New();
215 viennahrle::ConstSparseIterator<DomainType> maskIt(maskDomain);
216 for (auto &node : surfaceMesh->getNodes()) {
217 viennahrle::Index<D> index;
218 for (unsigned i = 0; i < D; ++i) {
219 index[i] = std::round(node[i] / gridDelta);
220 }
221 // can do sequential, because surfaceNodes are lexicographically sorted
222 // from lsToDiskMesh
223 maskIt.goToIndicesSequential(index);
224 // if it is a mask point, mark it to maybe use it in new level set
225 if (!maskIt.isDefined() || !(maskIt.getValue() < *valueIt + 1e-5)) {
226 newSurfaceMesh->insertNextNode(node);
227 newValues.push_back(*valueIt);
228 // insert vertex
229 std::array<unsigned, 1> vertex{};
230 vertex[0] = newSurfaceMesh->nodes.size();
231 newSurfaceMesh->insertNextVertex(vertex);
232 }
233 ++valueIt;
234 }
235 newSurfaceMesh->cellData.insertNextScalarData(newValues, "LSValues");
236 // use new mesh as surfaceMesh
237 newSurfaceMesh->minimumExtent = surfaceMesh->minimumExtent;
238 newSurfaceMesh->maximumExtent = surfaceMesh->maximumExtent;
239 surfaceMesh = newSurfaceMesh;
240 }
241
242#ifndef NDEBUG // if in debug build
243 {
244 Logger::getInstance()
245 .addDebug("GeomAdvect: Writing debug meshes")
246 .print();
248 "DEBUG_lsGeomAdvectMesh_contributewoMask.vtp")
249 .apply();
250 auto mesh = SmartPointer<Mesh<T>>::New();
251 if (maskLevelSet != nullptr) {
252 ToMesh<T, D>(maskLevelSet, mesh).apply();
254 "DEBUG_lsGeomAdvectMesh_mask.vtp")
255 .apply();
256 }
257 ToMesh<T, D>(levelSet, mesh).apply();
259 "DEBUG_lsGeomAdvectMesh_initial.vtp")
260 .apply();
261 }
262
263#endif
264
265 const auto &surfaceNodes = surfaceMesh->getNodes();
266
267 // initialize with segmentation for whole range
268 typename viennahrle::Domain<T, D>::IndexPoints segmentation;
269
270 {
271 unsigned long long numPoints = 1;
272 unsigned long long pointsPerDimension[D];
273 for (unsigned i = 0; i < D; ++i) {
274 pointsPerDimension[i] = numPoints;
275 numPoints *= max[i] - min[i];
276 }
277 unsigned long numberOfSegments = domain.getNumberOfSegments();
278 unsigned long long pointsPerSegment = numPoints / numberOfSegments;
279 unsigned long long pointId = 0;
280 for (unsigned i = 0; i < numberOfSegments - 1; ++i) {
281 pointId = pointsPerSegment * (i + 1);
282 viennahrle::Index<D> segmentPoint;
283 for (int j = D - 1; j >= 0; --j) {
284 segmentPoint[j] = pointId / (pointsPerDimension[j]) + min[j];
285 pointId %= pointsPerDimension[j];
286 }
287 segmentation.push_back(segmentPoint);
288 }
289 }
290
291 typedef std::vector<std::pair<viennahrle::Index<D>, T>> PointValueVector;
292 std::vector<PointValueVector> newPoints;
293 newPoints.resize(domain.getNumberOfSegments());
294
295 const T initialDistance = (distIsPositive)
296 ? std::numeric_limits<double>::max()
297 : std::numeric_limits<double>::lowest();
298
299#ifndef NDEBUG
300 {
301 std::ostringstream oss;
302 oss << "GeomAdvect: Min: " << min << ", Max: " << max << std::endl;
303 Logger::getInstance().addDebug(oss.str()).print();
304 }
305#endif
306// set up multithreading
307#pragma omp parallel num_threads(domain.getNumberOfSegments())
308 {
309 int p = 0;
310#ifdef _OPENMP
311 p = omp_get_thread_num();
312#endif
313
314 viennahrle::Index<D> startVector;
315 if (p == 0) {
316 startVector = min;
317 } else {
318 startVector = segmentation[p - 1];
319 incrementIndices(startVector, min, max);
320 }
321
322 viennahrle::Index<D> endVector =
323 (p != static_cast<int>(domain.getNumberOfSegments() - 1))
324 ? segmentation[p]
325 : grid.incrementIndices(max);
326
327 viennahrle::ConstSparseIterator<DomainType> checkIt(levelSet->getDomain(),
328 startVector);
329
330 // Mask iterator for checking whether inside mask or not
331 SmartPointer<viennahrle::ConstSparseIterator<DomainType>> maskIt =
332 nullptr;
333 if (maskLevelSet != nullptr) {
334 maskIt = SmartPointer<viennahrle::ConstSparseIterator<DomainType>>::New(
335 maskLevelSet->getDomain(), startVector);
336 }
337
338 // Iterate through the bounds of new lsDomain lexicographically
339 for (viennahrle::Index<D> currentIndex = startVector;
340 currentIndex <= endVector;
341 incrementIndices(currentIndex, min, max)) {
342 // if point is already full in old level set, skip it
343 checkIt.goToIndicesSequential(currentIndex);
344 T oldValue = checkIt.getValue();
345 // if run is already negative undefined, just ignore the point
346 if (distIsPositive) {
347 if (oldValue < -cutoffValue) {
348 continue;
349 }
350 } else if (oldValue > cutoffValue) {
351 continue;
352 }
353
354 VectorType<hrleCoordType, 3> currentCoords{};
355 VectorType<hrleCoordType, 3> currentDistMin{};
356 VectorType<hrleCoordType, 3> currentDistMax{};
357
358 for (unsigned i = 0; i < D; ++i) {
359 currentCoords[i] = currentIndex[i] * gridDelta;
360
361 currentDistMin[i] = currentIndex[i] - std::abs(distMin[i]);
362 if (currentDistMin[i] < grid.getMinGridPoint(i)) {
363 currentDistMin[i] = grid.getMinGridPoint(i);
364 }
365 currentDistMin[i] *= gridDelta;
366
367 currentDistMax[i] = currentIndex[i] + std::abs(distMax[i]);
368 if (currentDistMin[i] > grid.getMaxGridPoint(i)) {
369 currentDistMin[i] = grid.getMaxGridPoint(i);
370 }
371 currentDistMax[i] *= gridDelta;
372 }
373
374 T distance = initialDistance;
375
376 unsigned long currentPointId = 0;
377 // now check which surface points contribute to currentIndex
378 for (auto surfIt = surfaceNodes.begin(); surfIt != surfaceNodes.end();
379 ++surfIt, ++currentPointId) {
380
381 auto &currentNode = *surfIt;
382
383 // if we are outside min/max go to next index inside
384 {
385 bool outside = false;
386 for (unsigned i = 0; i < D; ++i) {
387 if ((currentNode[i] < currentDistMin[i]) ||
388 (currentNode[i] > currentDistMax[i])) {
389 outside = true;
390 break;
391 }
392 }
393 if (outside) {
394 continue;
395 }
396 }
397
398 // TODO: does this really save time? Try without it.
399 if (!dist->isInside(currentNode, currentCoords, 2 * gridDelta)) {
400 continue;
401 }
402
403 // get filling fraction from distance to dist surface
404 T tmpDistance = dist->getSignedDistance(
405 currentNode, currentCoords,
406 pointIdTranslator->find(currentPointId)->second) /
407 gridDelta;
408
409 // if cell is far within a distribution, set it filled
410 if (distIsPositive) {
411 if (tmpDistance <= -cutoffValue) {
412 distance = std::numeric_limits<T>::lowest();
413 break;
414 }
415
416 if (tmpDistance < distance) {
417 distance = tmpDistance;
418 }
419 } else {
420 if (tmpDistance >= cutoffValue) {
421 distance = std::numeric_limits<T>::max();
422 break;
423 }
424
425 if (tmpDistance > distance) {
426 distance = tmpDistance;
427 }
428 }
429 }
430
431 // TODO: There are still issues with positive box distributions
432 // if there is a mask used!
433 // if point is part of the mask, keep smaller value
434 if (maskLevelSet != nullptr) {
435 maskIt->goToIndicesSequential(currentIndex);
436
437 // if dist is positive, flip logic of comparison
438 if (distIsPositive ^
439 (std::abs(oldValue - maskIt->getValue()) < 1e-6)) {
440 if (!distIsPositive && std::abs(oldValue) <= cutoffValue) {
441 newPoints[p].push_back(std::make_pair(currentIndex, oldValue));
442 continue;
443 }
444 } else {
445 if (distance != initialDistance) {
446 distance = std::min(maskIt->getValue(), distance);
447 } else if (distIsPositive || oldValue >= 0.) {
448 newPoints[p].push_back(std::make_pair(currentIndex, oldValue));
449 continue;
450 }
451 }
452 }
453
454 if (std::abs(distance) <= cutoffValue) {
455 // avoid using distribution in wrong direction
456 if (distIsPositive && oldValue >= 0.) {
457 newPoints[p].push_back(std::make_pair(currentIndex, distance));
458 } else if (!distIsPositive && oldValue <= 0.) {
459 // if we are etching, need to make sure, we are not inside mask
460 if (maskIt == nullptr || maskIt->getValue() > -cutoffValue) {
461 newPoints[p].push_back(std::make_pair(currentIndex, distance));
462 }
463 } else {
464 // this only happens if distribution is very small, < 2 * gridDelta
465 newPoints[p].push_back(std::make_pair(currentIndex, oldValue));
466 }
467 }
468 }
469 }
470
471 // copy all points into the first vector
472 {
473 unsigned long long numberOfPoints = newPoints[0].size();
474 for (unsigned i = 1; i < domain.getNumberOfSegments(); ++i) {
475 numberOfPoints += newPoints[i].size();
476 }
477 newPoints[0].reserve(numberOfPoints);
478 for (unsigned i = 1; i < domain.getNumberOfSegments(); ++i) {
479 std::move(std::begin(newPoints[i]), std::end(newPoints[i]),
480 std::back_inserter(newPoints[0]));
481 }
482 }
483
484 auto mesh = SmartPointer<Mesh<T>>::New();
485 // output all points directly to mesh
486 {
487 std::vector<T> scalarData;
488 for (auto it = newPoints[0].begin(); it != newPoints[0].end(); ++it) {
489 Vec3D<T> node{0., 0., 0.};
490 for (unsigned i = 0; i < D; ++i) {
491 node[i] = T((it->first)[i]) * gridDelta;
492 }
493
494 mesh->insertNextNode(node);
495 std::array<unsigned, 1> vertex{};
496 vertex[0] = mesh->vertices.size();
497 mesh->insertNextVertex(vertex);
498 scalarData.push_back(it->second);
499 }
500 mesh->cellData.insertNextScalarData(scalarData, "LSValues");
501 }
502
503#ifndef NDEBUG // if in debug build
504 Logger::getInstance().addDebug("GeomAdvect: Writing final mesh...").print();
505 VTKWriter<T>(mesh, FileFormatEnum::VTP, "DEBUG_lsGeomAdvectMesh_final.vtp")
506 .apply();
507#endif
508
509 FromMesh<T, D>(levelSet, mesh).apply();
510
511#ifndef NDEBUG // if in debug build
512 Logger::getInstance().addDebug("GeomAdvect: Writing final LS...").print();
513 ToMesh<T, D>(levelSet, mesh).apply();
514 VTKWriter<T>(mesh, FileFormatEnum::VTP, "DEBUG_lsGeomAdvectLS_final.vtp")
515 .apply();
516#endif
517
518 Prune<T, D>(levelSet).apply();
519
520 levelSet->getDomain().segment();
521 levelSet->finalize(1);
522
523 Expand<T, D>(levelSet, 2).apply();
524 }
525};
526
527// add all template specialisations for this class
528PRECOMPILE_PRECISION_DIMENSION(GeometricAdvect)
529
530} // namespace viennals
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
DomainType & getDomain()
get const reference to the underlying hrleDomain data structure
Definition lsDomain.hpp:142
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
Import the regular grid, on which the level set values are defined, from an explicit Mesh<>....
Definition lsFromMesh.hpp:16
void apply()
Definition lsFromMesh.hpp:40
Base class for distributions used by lsGeometricAdvect. All functions are pure virtual and must be im...
Definition lsGeometricAdvectDistributions.hpp:15
void setLevelSet(SmartPointer< Domain< T, D > > passedLevelSet)
Set the levelset which should be advected.
Definition lsGeometricAdvect.hpp:82
void setMaskLevelSet(SmartPointer< Domain< T, D > > passedMaskLevelSet)
Set the levelset, which should be used as a mask. This level set has to be wrapped by the levelset se...
Definition lsGeometricAdvect.hpp:99
void setAdvectionDistribution(SmartPointer< DistType > passedDist)
Set which advection distribution to use. Must be derived from GeometricAdvectDistribution.
Definition lsGeometricAdvect.hpp:91
GeometricAdvect(SmartPointer< Domain< T, D > > passedLevelSet, SmartPointer< DistType > passedDist, SmartPointer< Domain< T, D > > passedMaskLevelSet=nullptr)
Definition lsGeometricAdvect.hpp:73
void apply()
Perform geometrical advection.
Definition lsGeometricAdvect.hpp:104
std::vector< T > ScalarDataType
Definition lsPointData.hpp:22
Removes all level set points, which do not have at least one oppositely signed neighbour (Meaning the...
Definition lsPrune.hpp:17
void apply()
removes all grid points, which do not have at least one opposite signed neighbour returns the number ...
Definition lsPrune.hpp:74
This class creates a mesh from the level set with all grid points with a level set value <= 0....
Definition lsToDiskMesh.hpp:24
void apply()
Definition lsToDiskMesh.hpp:83
Extract the regular grid, on which the level set values are defined, to an explicit Mesh<>....
Definition lsToMesh.hpp:19
void apply()
Definition lsToMesh.hpp:49
Class handling the output of an Mesh<> to VTK file types.
Definition lsVTKWriter.hpp:32
void apply()
Definition lsVTKWriter.hpp:92
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
constexpr AssignType assignable
Definition lsConcepts.hpp:10
std::enable_if_t< std::is_base_of_v< Base, Derived >, AssignType > IsBaseOf
Definition lsConcepts.hpp:13
Definition lsAdvect.hpp:36
@ VTP
Definition lsFileFormats.hpp:6
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68