96 if (levelSet ==
nullptr) {
98 .addError(
"No level set passed to GeometricAdvect. Not Advecting.")
102 if (dist ==
nullptr) {
103 Logger::getInstance()
104 .addError(
"No GeometricAdvectDistribution passed to "
105 "GeometricAdvect. Not "
114 if (maskLevelSet !=
nullptr) {
118 dist->prepare(levelSet);
124 auto &grid = levelSet->getGrid();
125 const auto gridDelta = grid.getGridDelta();
126 const bool useSurfacePointId = dist->useSurfacePointId();
130 auto surfaceMesh = SmartPointer<Mesh<hrleCoordType>>::New();
131 auto pointIdTranslator =
132 SmartPointer<typename ToDiskMesh<T, D>::TranslatorType>::New();
135 if (!useSurfacePointId)
136 *pointIdTranslator = inverseTranslator(*pointIdTranslator);
139 auto distBounds = dist->getBounds();
142 viennahrle::Index<D> distMin, distMax;
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;
151 hrleIndexType bounds[6];
152 domain.getDomainBounds(bounds);
153 viennahrle::Index<D> min, max;
154 for (
unsigned i = 0; i <
D; ++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;
167 min[i] = surfaceMesh->minimumExtent[i] / gridDelta;
169 if (grid.isNegBoundaryInfinite(i) && minPointNegative && distMin[i] < 0) {
172 if (distIsPositive) {
173 min[i] += distMin[i];
175 min[i] -= distMin[i];
180 if (min[i] < grid.getMinGridPoint(i)) {
181 min[i] = grid.getMinGridPoint(i);
184 max[i] = surfaceMesh->maximumExtent[i] / gridDelta;
185 if (grid.isPosBoundaryInfinite(i) && maxPointNegative && distMax[i] > 0) {
188 if (distIsPositive) {
189 max[i] += distMax[i];
191 max[i] -= distMax[i];
194 if (max[i] > grid.getMaxGridPoint(i)) {
195 max[i] = grid.getMaxGridPoint(i);
202 if (maskLevelSet !=
nullptr) {
204 auto &maskDomain = maskLevelSet->getDomain();
205 auto values = surfaceMesh->cellData.getScalarData(
"LSValues");
206 auto valueIt = values->begin();
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);
218 maskIt.goToIndicesSequential(index);
220 if (!maskIt.isDefined() || !(maskIt.getValue() < *valueIt + 1e-5)) {
221 newSurfaceMesh->insertNextNode(node);
222 newValues.push_back(*valueIt);
224 std::array<unsigned, 1> vertex{};
225 vertex[0] = newSurfaceMesh->nodes.size();
226 newSurfaceMesh->insertNextVertex(vertex);
230 newSurfaceMesh->cellData.insertNextScalarData(newValues,
"LSValues");
232 newSurfaceMesh->minimumExtent = surfaceMesh->minimumExtent;
233 newSurfaceMesh->maximumExtent = surfaceMesh->maximumExtent;
234 surfaceMesh = newSurfaceMesh;
239 Logger::getInstance()
240 .addDebug(
"GeomAdvect: Writing debug meshes")
243 "DEBUG_lsGeomAdvectMesh_contributewoMask.vtp")
245 auto mesh = SmartPointer<Mesh<T>>::New();
246 if (maskLevelSet !=
nullptr) {
249 "DEBUG_lsGeomAdvectMesh_mask.vtp")
254 "DEBUG_lsGeomAdvectMesh_initial.vtp")
260 const auto &surfaceNodes = surfaceMesh->getNodes();
263 typename viennahrle::Domain<T, D>::IndexPoints segmentation;
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];
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];
282 segmentation.push_back(segmentPoint);
286 typedef std::vector<std::pair<viennahrle::Index<D>,
T>> PointValueVector;
287 std::vector<PointValueVector> newPoints;
288 newPoints.resize(domain.getNumberOfSegments());
290 const T initialDistance = (distIsPositive)
291 ? std::numeric_limits<double>::max()
292 : std::numeric_limits<double>::lowest();
296 std::ostringstream oss;
297 oss <<
"GeomAdvect: Min: " << min <<
", Max: " << max << std::endl;
298 Logger::getInstance().addDebug(oss.str()).print();
302#pragma omp parallel num_threads(domain.getNumberOfSegments())
306 p = omp_get_thread_num();
309 viennahrle::Index<D> startVector;
313 startVector = segmentation[p - 1];
314 incrementIndices(startVector, min, max);
317 viennahrle::Index<D> endVector =
318 (p !=
static_cast<int>(domain.getNumberOfSegments() - 1))
320 : grid.incrementIndices(max);
322 viennahrle::ConstSparseIterator<DomainType> checkIt(levelSet->getDomain(),
326 SmartPointer<viennahrle::ConstSparseIterator<DomainType>> maskIt =
328 if (maskLevelSet !=
nullptr) {
329 maskIt = SmartPointer<viennahrle::ConstSparseIterator<DomainType>>::New(
330 maskLevelSet->getDomain(), startVector);
334 for (viennahrle::Index<D> currentIndex = startVector;
335 currentIndex <= endVector;
336 incrementIndices(currentIndex, min, max)) {
338 checkIt.goToIndicesSequential(currentIndex);
339 T oldValue = checkIt.getValue();
341 if (distIsPositive) {
342 if (oldValue < -cutoffValue) {
345 }
else if (oldValue > cutoffValue) {
349 VectorType<hrleCoordType, 3> currentCoords{};
350 VectorType<hrleCoordType, 3> currentDistMin{};
351 VectorType<hrleCoordType, 3> currentDistMax{};
353 for (
unsigned i = 0; i <
D; ++i) {
354 currentCoords[i] = currentIndex[i] * gridDelta;
356 currentDistMin[i] = currentIndex[i] - std::abs(distMin[i]);
357 if (currentDistMin[i] < grid.getMinGridPoint(i)) {
358 currentDistMin[i] = grid.getMinGridPoint(i);
360 currentDistMin[i] *= gridDelta;
362 currentDistMax[i] = currentIndex[i] + std::abs(distMax[i]);
363 if (currentDistMin[i] > grid.getMaxGridPoint(i)) {
364 currentDistMin[i] = grid.getMaxGridPoint(i);
366 currentDistMax[i] *= gridDelta;
369 T distance = initialDistance;
371 unsigned long currentPointId = 0;
373 for (
auto surfIt = surfaceNodes.begin(); surfIt != surfaceNodes.end();
374 ++surfIt, ++currentPointId) {
376 auto ¤tNode = *surfIt;
380 bool outside =
false;
381 for (
unsigned i = 0; i <
D; ++i) {
382 if ((currentNode[i] < currentDistMin[i]) ||
383 (currentNode[i] > currentDistMax[i])) {
394 if (!dist->isInside(currentNode, currentCoords, 2 * gridDelta)) {
399 auto pointId = currentPointId;
400 if (!useSurfacePointId) {
401 pointId = pointIdTranslator->find(currentPointId)->second;
404 dist->getSignedDistance(currentNode, currentCoords, pointId) /
408 if (distIsPositive) {
409 if (tmpDistance <= -cutoffValue) {
410 distance = std::numeric_limits<T>::lowest();
414 if (tmpDistance < distance) {
415 distance = tmpDistance;
418 if (tmpDistance >= cutoffValue) {
419 distance = std::numeric_limits<T>::max();
423 if (tmpDistance > distance) {
424 distance = tmpDistance;
432 if (maskLevelSet !=
nullptr) {
433 maskIt->goToIndicesSequential(currentIndex);
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));
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));
452 if (std::abs(distance) <= cutoffValue) {
454 if (distIsPositive && oldValue >= 0.) {
455 newPoints[p].push_back(std::make_pair(currentIndex, distance));
456 }
else if (!distIsPositive && oldValue <= 0.) {
458 if (maskIt ==
nullptr || maskIt->getValue() > -cutoffValue) {
459 newPoints[p].push_back(std::make_pair(currentIndex, distance));
463 newPoints[p].push_back(std::make_pair(currentIndex, oldValue));
471 unsigned long long numberOfPoints = newPoints[0].size();
472 for (
unsigned i = 1; i < domain.getNumberOfSegments(); ++i) {
473 numberOfPoints += newPoints[i].size();
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]));
482 auto mesh = SmartPointer<Mesh<T>>::New();
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;
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);
498 mesh->cellData.insertNextScalarData(scalarData,
"LSValues");
502 Logger::getInstance().addDebug(
"GeomAdvect: Writing final mesh...").print();
510 Logger::getInstance().addDebug(
"GeomAdvect: Writing final LS...").print();
518 levelSet->getDomain().segment();
519 levelSet->finalize(1);