105 if (levelSet ==
nullptr) {
106 Logger::getInstance()
107 .addWarning(
"No level set passed to GeometricAdvect. Not Advecting.")
111 if (dist ==
nullptr) {
112 Logger::getInstance()
113 .addWarning(
"No GeometricAdvectDistribution passed to "
114 "GeometricAdvect. Not "
123 if (maskLevelSet !=
nullptr) {
131 auto &grid = levelSet->getGrid();
132 auto gridDelta = grid.getGridDelta();
136 auto surfaceMesh = SmartPointer<Mesh<hrleCoordType>>::New();
137 auto pointIdTranslator =
138 SmartPointer<typename ToDiskMesh<T, D>::TranslatorType>::New();
141 *pointIdTranslator = inverseTranslator(*pointIdTranslator);
144 auto distBounds = dist->getBounds();
147 viennahrle::Index<D> distMin, distMax;
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;
156 hrleIndexType bounds[6];
157 domain.getDomainBounds(bounds);
158 viennahrle::Index<D> min, max;
159 for (
unsigned i = 0; i <
D; ++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;
172 min[i] = surfaceMesh->minimumExtent[i] / gridDelta;
174 if (grid.isNegBoundaryInfinite(i) && minPointNegative && distMin[i] < 0) {
177 if (distIsPositive) {
178 min[i] += distMin[i];
180 min[i] -= distMin[i];
185 if (min[i] < grid.getMinGridPoint(i)) {
186 min[i] = grid.getMinGridPoint(i);
189 max[i] = surfaceMesh->maximumExtent[i] / gridDelta;
190 if (grid.isPosBoundaryInfinite(i) && maxPointNegative && distMax[i] > 0) {
193 if (distIsPositive) {
194 max[i] += distMax[i];
196 max[i] -= distMax[i];
199 if (max[i] > grid.getMaxGridPoint(i)) {
200 max[i] = grid.getMaxGridPoint(i);
207 if (maskLevelSet !=
nullptr) {
209 auto &maskDomain = maskLevelSet->getDomain();
210 auto values = surfaceMesh->cellData.getScalarData(
"LSValues");
211 auto valueIt = values->begin();
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);
223 maskIt.goToIndicesSequential(index);
225 if (!maskIt.isDefined() || !(maskIt.getValue() < *valueIt + 1e-5)) {
226 newSurfaceMesh->insertNextNode(node);
227 newValues.push_back(*valueIt);
229 std::array<unsigned, 1> vertex{};
230 vertex[0] = newSurfaceMesh->nodes.size();
231 newSurfaceMesh->insertNextVertex(vertex);
235 newSurfaceMesh->cellData.insertNextScalarData(newValues,
"LSValues");
237 newSurfaceMesh->minimumExtent = surfaceMesh->minimumExtent;
238 newSurfaceMesh->maximumExtent = surfaceMesh->maximumExtent;
239 surfaceMesh = newSurfaceMesh;
244 Logger::getInstance()
245 .addDebug(
"GeomAdvect: Writing debug meshes")
248 "DEBUG_lsGeomAdvectMesh_contributewoMask.vtp")
250 auto mesh = SmartPointer<Mesh<T>>::New();
251 if (maskLevelSet !=
nullptr) {
254 "DEBUG_lsGeomAdvectMesh_mask.vtp")
259 "DEBUG_lsGeomAdvectMesh_initial.vtp")
265 const auto &surfaceNodes = surfaceMesh->getNodes();
268 typename viennahrle::Domain<T, D>::IndexPoints segmentation;
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];
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];
287 segmentation.push_back(segmentPoint);
291 typedef std::vector<std::pair<viennahrle::Index<D>,
T>> PointValueVector;
292 std::vector<PointValueVector> newPoints;
293 newPoints.resize(domain.getNumberOfSegments());
295 const T initialDistance = (distIsPositive)
296 ? std::numeric_limits<double>::max()
297 : std::numeric_limits<double>::lowest();
301 std::ostringstream oss;
302 oss <<
"GeomAdvect: Min: " << min <<
", Max: " << max << std::endl;
303 Logger::getInstance().addDebug(oss.str()).print();
307#pragma omp parallel num_threads(domain.getNumberOfSegments())
311 p = omp_get_thread_num();
314 viennahrle::Index<D> startVector;
318 startVector = segmentation[p - 1];
319 incrementIndices(startVector, min, max);
322 viennahrle::Index<D> endVector =
323 (p !=
static_cast<int>(domain.getNumberOfSegments() - 1))
325 : grid.incrementIndices(max);
327 viennahrle::ConstSparseIterator<DomainType> checkIt(levelSet->getDomain(),
331 SmartPointer<viennahrle::ConstSparseIterator<DomainType>> maskIt =
333 if (maskLevelSet !=
nullptr) {
334 maskIt = SmartPointer<viennahrle::ConstSparseIterator<DomainType>>::New(
335 maskLevelSet->getDomain(), startVector);
339 for (viennahrle::Index<D> currentIndex = startVector;
340 currentIndex <= endVector;
341 incrementIndices(currentIndex, min, max)) {
343 checkIt.goToIndicesSequential(currentIndex);
344 T oldValue = checkIt.getValue();
346 if (distIsPositive) {
347 if (oldValue < -cutoffValue) {
350 }
else if (oldValue > cutoffValue) {
354 VectorType<hrleCoordType, 3> currentCoords{};
355 VectorType<hrleCoordType, 3> currentDistMin{};
356 VectorType<hrleCoordType, 3> currentDistMax{};
358 for (
unsigned i = 0; i <
D; ++i) {
359 currentCoords[i] = currentIndex[i] * gridDelta;
361 currentDistMin[i] = currentIndex[i] - std::abs(distMin[i]);
362 if (currentDistMin[i] < grid.getMinGridPoint(i)) {
363 currentDistMin[i] = grid.getMinGridPoint(i);
365 currentDistMin[i] *= gridDelta;
367 currentDistMax[i] = currentIndex[i] + std::abs(distMax[i]);
368 if (currentDistMin[i] > grid.getMaxGridPoint(i)) {
369 currentDistMin[i] = grid.getMaxGridPoint(i);
371 currentDistMax[i] *= gridDelta;
374 T distance = initialDistance;
376 unsigned long currentPointId = 0;
378 for (
auto surfIt = surfaceNodes.begin(); surfIt != surfaceNodes.end();
379 ++surfIt, ++currentPointId) {
381 auto ¤tNode = *surfIt;
385 bool outside =
false;
386 for (
unsigned i = 0; i <
D; ++i) {
387 if ((currentNode[i] < currentDistMin[i]) ||
388 (currentNode[i] > currentDistMax[i])) {
399 if (!dist->isInside(currentNode, currentCoords, 2 * gridDelta)) {
404 T tmpDistance = dist->getSignedDistance(
405 currentNode, currentCoords,
406 pointIdTranslator->find(currentPointId)->second) /
410 if (distIsPositive) {
411 if (tmpDistance <= -cutoffValue) {
412 distance = std::numeric_limits<T>::lowest();
416 if (tmpDistance < distance) {
417 distance = tmpDistance;
420 if (tmpDistance >= cutoffValue) {
421 distance = std::numeric_limits<T>::max();
425 if (tmpDistance > distance) {
426 distance = tmpDistance;
434 if (maskLevelSet !=
nullptr) {
435 maskIt->goToIndicesSequential(currentIndex);
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));
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));
454 if (std::abs(distance) <= cutoffValue) {
456 if (distIsPositive && oldValue >= 0.) {
457 newPoints[p].push_back(std::make_pair(currentIndex, distance));
458 }
else if (!distIsPositive && oldValue <= 0.) {
460 if (maskIt ==
nullptr || maskIt->getValue() > -cutoffValue) {
461 newPoints[p].push_back(std::make_pair(currentIndex, distance));
465 newPoints[p].push_back(std::make_pair(currentIndex, oldValue));
473 unsigned long long numberOfPoints = newPoints[0].size();
474 for (
unsigned i = 1; i < domain.getNumberOfSegments(); ++i) {
475 numberOfPoints += newPoints[i].size();
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]));
484 auto mesh = SmartPointer<Mesh<T>>::New();
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;
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);
500 mesh->cellData.insertNextScalarData(scalarData,
"LSValues");
504 Logger::getInstance().addDebug(
"GeomAdvect: Writing final mesh...").print();
512 Logger::getInstance().addDebug(
"GeomAdvect: Writing final LS...").print();
520 levelSet->getDomain().segment();
521 levelSet->finalize(1);