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