ViennaLS
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
lsAdvect.hpp
Go to the documentation of this file.
1#pragma once
2
4
5#include <limits>
6#include <vector>
7
8#include <hrleSparseIterator.hpp>
9#include <hrleSparseStarIterator.hpp>
10
11#include <vcLogger.hpp>
12#include <vcSmartPointer.hpp>
13
15#include <lsDomain.hpp>
16#include <lsMarkVoidPoints.hpp>
17#include <lsReduce.hpp>
18
19// Integration schemes
20#include <lsEngquistOsher.hpp>
21#include <lsLaxFriedrichs.hpp>
26
27// Velocity accessor
28#include <lsVelocityField.hpp>
29
30// #define DEBUG_LS_ADVECT_HPP
31#ifdef DEBUG_LS_ADVECT_HPP
32#include <lsToMesh.hpp>
33#include <lsVTKWriter.hpp>
34#endif
35
36namespace viennals {
37
38using namespace viennacore;
39
54
64template <class T, int D> class Advect {
65 using ConstSparseIterator =
66 viennahrle::ConstSparseIterator<typename Domain<T, D>::DomainType>;
67 using hrleIndexType = viennahrle::IndexType;
68
69 std::vector<SmartPointer<Domain<T, D>>> levelSets;
70 SmartPointer<VelocityField<T>> velocities = nullptr;
71 IntegrationSchemeEnum integrationScheme =
73 double timeStepRatio = 0.4999;
74 double dissipationAlpha = 1.0;
75 bool calculateNormalVectors = true;
76 bool ignoreVoids = false;
77 double advectionTime = 0.;
78 bool performOnlySingleStep = false;
79 double advectedTime = 0.;
80 unsigned numberOfTimeSteps = 0;
81 bool saveAdvectionVelocities = false;
82 bool updatePointData = true;
83 bool checkDissipation = true;
84 static constexpr double wrappingLayerEpsilon = 1e-4;
85
86 // this vector will hold the maximum time step for each point and the
87 // corresponding velocity
88 std::vector<std::vector<std::pair<std::pair<T, T>, T>>> storedRates;
89 double currentTimeStep = -1.;
90
91 VectorType<T, 3> findGlobalAlphas() const {
92
93 auto &topDomain = levelSets.back()->getDomain();
94 auto &grid = levelSets.back()->getGrid();
95
96 const T gridDelta = grid.getGridDelta();
97 const T deltaPos = gridDelta;
98 const T deltaNeg = -gridDelta;
99
100 VectorType<T, 3> finalAlphas = {0., 0., 0.};
101
102#pragma omp parallel num_threads((levelSets.back())->getNumberOfSegments())
103 {
104 VectorType<T, 3> localAlphas = {0., 0., 0.};
105 int p = 0;
106#ifdef _OPENMP
107 p = omp_get_thread_num();
108#endif
109 viennahrle::Index<D> startVector =
110 (p == 0) ? grid.getMinGridPoint()
111 : topDomain.getSegmentation()[p - 1];
112 viennahrle::Index<D> endVector =
113 (p != static_cast<int>(topDomain.getNumberOfSegments() - 1))
114 ? topDomain.getSegmentation()[p]
115 : grid.incrementIndices(grid.getMaxGridPoint());
116
117 // an iterator for each level set
118 std::vector<ConstSparseIterator> iterators;
119 for (auto const &ls : levelSets) {
120 iterators.emplace_back(ls->getDomain());
121 }
122
123 // neighborIterator for the top level set
124 viennahrle::ConstSparseStarIterator<typename Domain<T, D>::DomainType, 1>
125 neighborIterator(topDomain);
126
127 for (ConstSparseIterator it(topDomain, startVector);
128 it.getStartIndices() < endVector; ++it) {
129
130 if (!it.isDefined() || std::abs(it.getValue()) > 0.5)
131 continue;
132
133 const T value = it.getValue();
134 const auto indices = it.getStartIndices();
135
136 // check if there is any other levelset at the same point:
137 // if yes, take the velocity of the lowest levelset
138 for (unsigned lowerLevelSetId = 0; lowerLevelSetId < levelSets.size();
139 ++lowerLevelSetId) {
140 // put iterator to same position as the top levelset
141 iterators[lowerLevelSetId].goToIndicesSequential(indices);
142
143 // if the lower surface is actually outside, i.e. its LS value
144 // is lower or equal
145 if (iterators[lowerLevelSetId].getValue() <=
146 value + wrappingLayerEpsilon) {
147
148 // move neighborIterator to current position
149 neighborIterator.goToIndicesSequential(indices);
150
151 Vec3D<T> coords;
152 for (unsigned i = 0; i < D; ++i) {
153 coords[i] = indices[i] * gridDelta;
154 }
155
156 Vec3D<T> normal = {};
157 T normalModulus = 0.;
158 for (unsigned i = 0; i < D; ++i) {
159 const T phiPos = neighborIterator.getNeighbor(i).getValue();
160 const T phiNeg = neighborIterator.getNeighbor(i + D).getValue();
161
162 T diffPos = (phiPos - value) / deltaPos;
163 T diffNeg = (phiNeg - value) / deltaNeg;
164
165 normal[i] = (diffNeg + diffPos) * 0.5;
166 normalModulus += normal[i] * normal[i];
167 }
168 normalModulus = std::sqrt(normalModulus);
169 for (unsigned i = 0; i < D; ++i)
170 normal[i] /= normalModulus;
171
172 T scaVel = velocities->getScalarVelocity(
173 coords, lowerLevelSetId, normal,
174 neighborIterator.getCenter().getPointId());
175 auto vecVel = velocities->getVectorVelocity(
176 coords, lowerLevelSetId, normal,
177 neighborIterator.getCenter().getPointId());
178
179 for (unsigned i = 0; i < D; ++i) {
180 T tempAlpha = std::abs((scaVel + vecVel[i]) * normal[i]);
181 localAlphas[i] = std::max(localAlphas[i], tempAlpha);
182 }
183
184 // exit material loop
185 break;
186 }
187 }
188 }
189
190#pragma omp critical
191 {
192 for (unsigned i = 0; i < D; ++i) {
193 finalAlphas[i] = std::max(finalAlphas[i], localAlphas[i]);
194 }
195 }
196 } // end of parallel section
197
198 return finalAlphas;
199 }
200
201 void rebuildLS() {
202 // TODO: this function uses Manhattan distances for renormalisation,
203 // since this is the quickest. For visualisation applications, better
204 // renormalisation is needed, so it might be good to implement
205 // Euler distance renormalisation as an option
206 auto &grid = levelSets.back()->getGrid();
207 auto newlsDomain = SmartPointer<Domain<T, D>>::New(grid);
208 auto &newDomain = newlsDomain->getDomain();
209 auto &domain = levelSets.back()->getDomain();
210
211 newDomain.initialize(domain.getNewSegmentation(),
212 domain.getAllocation() *
213 (2.0 / levelSets.back()->getLevelSetWidth()));
214
215 const bool updateData = updatePointData;
216 // save how data should be transferred to new level set
217 // list of indices into the old pointData vector
218 std::vector<std::vector<unsigned>> newDataSourceIds;
219 if (updateData)
220 newDataSourceIds.resize(newDomain.getNumberOfSegments());
221
222#ifdef DEBUG_LS_ADVECT_HPP
223 {
224 auto mesh = SmartPointer<Mesh<T>>::New();
225 ToMesh<T, D>(levelSets.back(), mesh).apply();
226 VTKWriter<T>(mesh, "Advect_beforeRebuild.vtk").apply();
227 }
228#endif
229
230#pragma omp parallel num_threads(newDomain.getNumberOfSegments())
231 {
232 int p = 0;
233#ifdef _OPENMP
234 p = omp_get_thread_num();
235#endif
236 auto &domainSegment = newDomain.getDomainSegment(p);
237
238 viennahrle::Index<D> startVector =
239 (p == 0) ? grid.getMinGridPoint()
240 : newDomain.getSegmentation()[p - 1];
241
242 viennahrle::Index<D> endVector =
243 (p != static_cast<int>(newDomain.getNumberOfSegments() - 1))
244 ? newDomain.getSegmentation()[p]
245 : grid.incrementIndices(grid.getMaxGridPoint());
246
247 // reserve a bit more to avoid reallocation
248 // would expect number of points to roughly double
249 if (updateData)
250 newDataSourceIds[p].reserve(2.5 * domainSegment.getNumberOfPoints());
251
252 for (viennahrle::SparseStarIterator<typename Domain<T, D>::DomainType, 1>
253 it(domain, startVector);
254 it.getIndices() < endVector; ++it) {
255
256 // if the center is an active grid point
257 // <1.0 since it could have been change by 0.5 max
258 if (std::abs(it.getCenter().getValue()) <= 1.0) {
259
260 int k = 0;
261 for (; k < 2 * D; k++)
262 if (std::signbit(it.getNeighbor(k).getValue() - 1e-7) !=
263 std::signbit(it.getCenter().getValue() + 1e-7))
264 break;
265
266 // if there is at least one neighbor of opposite sign
267 if (k != 2 * D) {
268 if (it.getCenter().getDefinedValue() > 0.5) {
269 int j = 0;
270 for (; j < 2 * D; j++) {
271 if (std::abs(it.getNeighbor(j).getValue()) <= 1.0)
272 if (it.getNeighbor(j).getDefinedValue() < -0.5)
273 break;
274 }
275 if (j == 2 * D) {
276 domainSegment.insertNextDefinedPoint(
277 it.getIndices(), it.getCenter().getDefinedValue());
278 if (updateData)
279 newDataSourceIds[p].push_back(it.getCenter().getPointId());
280 // if there is at least one active grid point, which is < -0.5
281 } else {
282 domainSegment.insertNextDefinedPoint(it.getIndices(), 0.5);
283 if (updateData)
284 newDataSourceIds[p].push_back(it.getNeighbor(j).getPointId());
285 }
286 } else if (it.getCenter().getDefinedValue() < -0.5) {
287 int j = 0;
288 for (; j < 2 * D; j++) {
289 if (std::abs(it.getNeighbor(j).getValue()) <= 1.0)
290 if (it.getNeighbor(j).getDefinedValue() > 0.5)
291 break;
292 }
293
294 if (j == 2 * D) {
295 domainSegment.insertNextDefinedPoint(
296 it.getIndices(), it.getCenter().getDefinedValue());
297 if (updateData)
298 newDataSourceIds[p].push_back(it.getCenter().getPointId());
299 // if there is at least one active grid point, which is > 0.5
300 } else {
301 domainSegment.insertNextDefinedPoint(it.getIndices(), -0.5);
302 if (updateData)
303 newDataSourceIds[p].push_back(it.getNeighbor(j).getPointId());
304 }
305 } else {
306 domainSegment.insertNextDefinedPoint(
307 it.getIndices(), it.getCenter().getDefinedValue());
308 if (updateData)
309 newDataSourceIds[p].push_back(it.getCenter().getPointId());
310 }
311 } else {
312 domainSegment.insertNextUndefinedPoint(
313 it.getIndices(), (it.getCenter().getDefinedValue() < 0)
316 }
317
318 } else { // if the center is not an active grid point
319 if (it.getCenter().getValue() >= 0) {
320 int usedNeighbor = -1;
321 T distance = Domain<T, D>::POS_VALUE;
322 for (int i = 0; i < 2 * D; i++) {
323 T value = it.getNeighbor(i).getValue();
324 if (std::abs(value) <= 1.0 && (value < 0.)) {
325 if (distance > value + 1.0) {
326 distance = value + 1.0;
327 usedNeighbor = i;
328 }
329 }
330 }
331
332 if (distance <= 1.) {
333 domainSegment.insertNextDefinedPoint(it.getIndices(), distance);
334 if (updateData)
335 newDataSourceIds[p].push_back(
336 it.getNeighbor(usedNeighbor).getPointId());
337 } else {
338 domainSegment.insertNextUndefinedPoint(it.getIndices(),
340 }
341
342 } else {
343 int usedNeighbor = -1;
344 T distance = Domain<T, D>::NEG_VALUE;
345 for (int i = 0; i < 2 * D; i++) {
346 T value = it.getNeighbor(i).getValue();
347 if (std::abs(value) <= 1.0 && (value > 0)) {
348 if (distance < value - 1.0) {
349 // distance = std::max(distance, value - T(1.0));
350 distance = value - 1.0;
351 usedNeighbor = i;
352 }
353 }
354 }
355
356 if (distance >= -1.) {
357 domainSegment.insertNextDefinedPoint(it.getIndices(), distance);
358 if (updateData)
359 newDataSourceIds[p].push_back(
360 it.getNeighbor(usedNeighbor).getPointId());
361 } else {
362 domainSegment.insertNextUndefinedPoint(it.getIndices(),
364 }
365 }
366 }
367 }
368 }
369
370 // now copy old data into new level set
371 if (updateData) {
372 auto &pointData = levelSets.back()->getPointData();
373 newlsDomain->getPointData().translateFromMultiData(pointData,
374 newDataSourceIds);
375 }
376
377 newDomain.finalize();
378 newDomain.segment();
379 levelSets.back()->deepCopy(newlsDomain);
380 levelSets.back()->finalize(2);
381 }
382
385 double advect(double maxTimeStep) {
386 // check whether a level set and velocities have been given
387 if (levelSets.empty()) {
388 Logger::getInstance()
389 .addWarning("No level sets passed to Advect. Not advecting.")
390 .print();
391 return std::numeric_limits<double>::max();
392 }
393 if (velocities == nullptr) {
394 Logger::getInstance()
395 .addWarning("No velocity field passed to Advect. Not advecting.")
396 .print();
397 return std::numeric_limits<double>::max();
398 }
399
400 if (currentTimeStep < 0. || storedRates.empty())
401 applyIntegration(maxTimeStep);
402
403 moveSurface();
404
405 rebuildLS();
406
407 // Adjust all level sets below the advected one
408 // This means, that when the top level-set and one below
409 // are etched, the lower one is moved with the top level-set
410 // TODO: Adjust lower layers also when they have grown,
411 // to allow for two different growth rates of materials
412 if (integrationScheme !=
414 for (unsigned i = 0; i < levelSets.size() - 1; ++i) {
415 BooleanOperation<T, D>(levelSets[i], levelSets.back(),
417 .apply();
418 }
419 }
420
421 return getCurrentTimeStep();
422 }
423
427 template <class IntegrationSchemeType>
428 double integrateTime(IntegrationSchemeType IntegrationScheme,
429 double maxTimeStep) {
430
431 auto &topDomain = levelSets.back()->getDomain();
432 auto &grid = levelSets.back()->getGrid();
433
434 typename PointData<T>::ScalarDataType *voidMarkerPointer;
435 if (ignoreVoids) {
436 MarkVoidPoints<T, D>(levelSets.back()).apply();
437 auto &pointData = levelSets.back()->getPointData();
438 voidMarkerPointer =
440 if (voidMarkerPointer == nullptr) {
441 Logger::getInstance()
442 .addWarning("Advect: Cannot find void point markers. Not "
443 "ignoring void points.")
444 .print();
445 ignoreVoids = false;
446 }
447 }
448 const bool ignoreVoidPoints = ignoreVoids;
449
450 if (!storedRates.empty()) {
451 Logger::getInstance()
452 .addWarning("Advect: Overwriting previously stored rates.")
453 .print();
454 }
455
456 storedRates.resize(topDomain.getNumberOfSegments());
457
458#pragma omp parallel num_threads(topDomain.getNumberOfSegments())
459 {
460 int p = 0;
461#ifdef _OPENMP
462 p = omp_get_thread_num();
463#endif
464 viennahrle::Index<D> startVector =
465 (p == 0) ? grid.getMinGridPoint()
466 : topDomain.getSegmentation()[p - 1];
467
468 viennahrle::Index<D> endVector =
469 (p != static_cast<int>(topDomain.getNumberOfSegments() - 1))
470 ? topDomain.getSegmentation()[p]
471 : grid.incrementIndices(grid.getMaxGridPoint());
472
473 double tempMaxTimeStep = maxTimeStep;
474 auto &tempRates = storedRates[p];
475 tempRates.reserve(
476 topDomain.getNumberOfPoints() /
477 static_cast<double>((levelSets.back())->getNumberOfSegments()) +
478 10);
479
480 // an iterator for each level set
481 std::vector<ConstSparseIterator> iterators;
482 for (auto const &ls : levelSets) {
483 iterators.emplace_back(ls->getDomain());
484 }
485
486 IntegrationSchemeType scheme(IntegrationScheme);
487
488 for (ConstSparseIterator it(topDomain, startVector);
489 it.getStartIndices() < endVector; ++it) {
490
491 if (!it.isDefined() || std::abs(it.getValue()) > 0.5)
492 continue;
493
494 T value = it.getValue();
495 double maxStepTime = 0;
496 double cfl = timeStepRatio;
497
498 for (int currentLevelSetId = levelSets.size() - 1;
499 currentLevelSetId >= 0; --currentLevelSetId) {
500
501 std::pair<T, T> gradNDissipation;
502
503 if (!(ignoreVoidPoints && (*voidMarkerPointer)[it.getPointId()])) {
504 // check if there is any other levelset at the same point:
505 // if yes, take the velocity of the lowest levelset
506 for (unsigned lowerLevelSetId = 0;
507 lowerLevelSetId < levelSets.size(); ++lowerLevelSetId) {
508 // put iterator to same position as the top levelset
509 iterators[lowerLevelSetId].goToIndicesSequential(
510 it.getStartIndices());
511
512 // if the lower surface is actually outside, i.e. its LS value
513 // is lower or equal
514 if (iterators[lowerLevelSetId].getValue() <=
515 value + wrappingLayerEpsilon) {
516 gradNDissipation =
517 scheme(it.getStartIndices(), lowerLevelSetId);
518 break;
519 }
520 }
521 }
522
523 T valueBelow;
524 // get value of material below (before in levelSets list)
525 if (currentLevelSetId > 0) {
526 iterators[currentLevelSetId - 1].goToIndicesSequential(
527 it.getStartIndices());
528 valueBelow = iterators[currentLevelSetId - 1].getValue();
529 } else {
530 valueBelow = std::numeric_limits<T>::max();
531 }
532
533 // if velocity is positive, set maximum time step possible without
534 // violating the cfl condition
535 T velocity = gradNDissipation.first - gradNDissipation.second;
536 if (velocity > 0.) {
537 maxStepTime += cfl / velocity;
538 tempRates.push_back(std::make_pair(gradNDissipation,
539 -std::numeric_limits<T>::max()));
540 break;
541 // if velocity is 0, maximum time step is infinite
542 } else if (velocity == 0.) {
543 maxStepTime = std::numeric_limits<T>::max();
544 tempRates.push_back(std::make_pair(gradNDissipation,
545 std::numeric_limits<T>::max()));
546 break;
547 // if the velocity is negative apply the velocity for as long as
548 // possible without infringing on material below
549 } else {
550 T difference = std::abs(valueBelow - value);
551
552 if (difference >= cfl) {
553 maxStepTime -= cfl / velocity;
554 tempRates.push_back(std::make_pair(
555 gradNDissipation, std::numeric_limits<T>::max()));
556 break;
557 } else {
558 maxStepTime -= difference / velocity;
559 // the second part of the pair indicates how far we can move
560 // in this time step until the end of the material is reached
561 tempRates.push_back(std::make_pair(gradNDissipation, valueBelow));
562 cfl -= difference;
563 // use new LS value for next calculations
564 value = valueBelow;
565 }
566 }
567 }
568
569 if (maxStepTime < tempMaxTimeStep)
570 tempMaxTimeStep = maxStepTime;
571 }
572
573#pragma omp critical
574 {
575 // If a Lax Friedrichs scheme is selected the time step is
576 // reduced depending on the dissipation coefficients
577 // For Engquist Osher scheme this function is empty.
578 scheme.reduceTimeStepHamiltonJacobi(
579 tempMaxTimeStep, levelSets.back()->getGrid().getGridDelta());
580
581 // set global timestep maximum
582 if (tempMaxTimeStep < maxTimeStep)
583 maxTimeStep = tempMaxTimeStep;
584 }
585 } // end of parallel section
586
587 // maxTimeStep is now the maximum time step possible for all points
588 // and rates are stored in a vector
589 return maxTimeStep;
590 }
591
592 // Level Sets below are also considered in order to adjust the advection
593 // depth accordingly if there would be a material change.
594 void moveSurface() {
595 if (timeStepRatio >= 0.5) {
596 Logger::getInstance()
597 .addWarning("Integration time step ratio should be smaller than 0.5. "
598 "Advection might fail!")
599 .print();
600 }
601
602 auto &topDomain = levelSets.back()->getDomain();
603 auto &grid = levelSets.back()->getGrid();
604
605 assert(currentTimeStep >= 0. && "No time step set!");
606 assert(storedRates.size() == topDomain.getNumberOfSegments());
607
608 // reduce to one layer thickness and apply new values directly to the
609 // domain segments --> DO NOT CHANGE SEGMENTATION HERE (true parameter)
610 Reduce<T, D>(levelSets.back(), 1, true).apply();
611
612 const bool saveVelocities = saveAdvectionVelocities;
613 std::vector<std::vector<double>> dissipationVectors(
614 levelSets.back()->getNumberOfSegments());
615 std::vector<std::vector<double>> velocityVectors(
616 levelSets.back()->getNumberOfSegments());
617
618 const bool checkDiss = checkDissipation;
619
620#pragma omp parallel num_threads(topDomain.getNumberOfSegments())
621 {
622 int p = 0;
623#ifdef _OPENMP
624 p = omp_get_thread_num();
625#endif
626 auto itRS = storedRates[p].cbegin();
627 auto &segment = topDomain.getDomainSegment(p);
628 const unsigned maxId = segment.getNumberOfPoints();
629
630 if (saveVelocities) {
631 velocityVectors[p].resize(maxId);
632 dissipationVectors[p].resize(maxId);
633 }
634
635 for (unsigned localId = 0; localId < maxId; ++localId) {
636 T &value = segment.definedValues[localId];
637 double time = currentTimeStep;
638
639 // if there is a change in materials during one time step, deduct
640 // the time taken to advect up to the end of the top material and
641 // set the LS value to the one below
642 auto const [gradient, dissipation] = itRS->first;
643 T velocity = gradient - dissipation;
644 // check if dissipation is too high
645 if (checkDiss && (gradient < 0 && velocity > 0) ||
646 (gradient > 0 && velocity < 0)) {
647 velocity = 0;
648 }
649
650 T rate = time * velocity;
651 while (std::abs(itRS->second - value) < std::abs(rate)) {
652 time -= std::abs((itRS->second - value) / velocity);
653 value = itRS->second;
654 ++itRS; // advance the TempStopRates iterator by one
655
656 // recalculate velocity and rate
657 velocity = itRS->first.first - itRS->first.second;
658 if (checkDiss && (itRS->first.first < 0 && velocity > 0) ||
659 (itRS->first.first > 0 && velocity < 0)) {
660 velocity = 0;
661 }
662 rate = time * velocity;
663 }
664
665 // now deduct the velocity times the time step we take
666 value -= rate;
667
668 if (saveVelocities) {
669 velocityVectors[p][localId] = rate;
670 dissipationVectors[p][localId] = itRS->first.second;
671 }
672
673 // this is run when two materials are close but the velocity is too slow
674 // to actually reach the second material, to get rid of the extra
675 // entry in the TempRatesStop
676 while (std::abs(itRS->second) != std::numeric_limits<T>::max())
677 ++itRS;
678
679 // advance the TempStopRates iterator by one
680 ++itRS;
681 }
682 } // end of parallel section
683
684 if (saveVelocities) {
685 auto &pointData = levelSets.back()->getPointData();
686
687 typename PointData<T>::ScalarDataType vels;
688 typename PointData<T>::ScalarDataType diss;
689
690 for (unsigned i = 0; i < velocityVectors.size(); ++i) {
691 vels.insert(vels.end(),
692 std::make_move_iterator(velocityVectors[i].begin()),
693 std::make_move_iterator(velocityVectors[i].end()));
694 diss.insert(diss.end(),
695 std::make_move_iterator(dissipationVectors[i].begin()),
696 std::make_move_iterator(dissipationVectors[i].end()));
697 }
698 pointData.insertReplaceScalarData(std::move(vels), velocityLabel);
699 pointData.insertReplaceScalarData(std::move(diss), dissipationLabel);
700 }
701
702 // clear the stored rates since surface has changed
703 storedRates.clear();
704 }
705
706public:
707 static constexpr char velocityLabel[] = "AdvectionVelocities";
708 static constexpr char dissipationLabel[] = "Dissipation";
709
710 Advect() = default;
711
712 explicit Advect(SmartPointer<Domain<T, D>> passedlsDomain) {
713 levelSets.push_back(passedlsDomain);
714 }
715
716 Advect(SmartPointer<Domain<T, D>> passedlsDomain,
717 SmartPointer<VelocityField<T>> passedVelocities) {
718 levelSets.push_back(passedlsDomain);
719 velocities = passedVelocities;
720 }
721
722 Advect(std::vector<SmartPointer<Domain<T, D>>> passedlsDomains,
723 SmartPointer<VelocityField<T>> passedVelocities)
724 : levelSets(passedlsDomains) {
725 velocities = passedVelocities;
726 }
727
730 void insertNextLevelSet(SmartPointer<Domain<T, D>> passedlsDomain) {
731 levelSets.push_back(passedlsDomain);
732 }
733
736 void setVelocityField(SmartPointer<VelocityField<T>> passedVelocities) {
737 velocities = passedVelocities;
738 }
739
745 void setAdvectionTime(double time) { advectionTime = time; }
746
751 void setSingleStep(bool singleStep) { performOnlySingleStep = singleStep; }
752
757 void setTimeStepRatio(const double &cfl) { timeStepRatio = cfl; }
758
763 void setCalculateNormalVectors(bool cnv) { calculateNormalVectors = cnv; }
764
771 void setIgnoreVoids(bool iV) { ignoreVoids = iV; }
772
775 void setSaveAdvectionVelocities(bool sAV) { saveAdvectionVelocities = sAV; }
776
779 double getAdvectedTime() const { return advectedTime; }
780
782 double getCurrentTimeStep() const { return currentTimeStep; }
783
785 unsigned getNumberOfTimeSteps() const { return numberOfTimeSteps; }
786
788 double getTimeStepRatio() const { return timeStepRatio; }
789
791 bool getCalculateNormalVectors() const { return calculateNormalVectors; }
792
796 integrationScheme = scheme;
797 }
798
803 void setDissipationAlpha(const double &a) { dissipationAlpha = a; }
804
805 // Sets the velocity to 0 if the dissipation is too high
806 void setCheckDissipation(bool check) { checkDissipation = check; }
807
810 void setUpdatePointData(bool update) { updatePointData = update; }
811
812 // Prepare the levelset for advection, based on the provided integration
813 // scheme.
814 void prepareLS() {
815 // check whether a level set and velocities have been given
816 if (levelSets.empty()) {
817 Logger::getInstance()
818 .addWarning("No level sets passed to Advect.")
819 .print();
820 return;
821 }
822
823 if (integrationScheme == IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER) {
825 } else if (integrationScheme ==
828 } else if (integrationScheme ==
831 } else if (integrationScheme ==
834 } else if (integrationScheme ==
838 levelSets.back());
839 } else if (integrationScheme ==
842 } else if (integrationScheme ==
845 } else if (integrationScheme ==
848 } else if (integrationScheme ==
851 } else if (integrationScheme ==
854 levelSets.back());
855 } else {
856 Logger::getInstance()
857 .addWarning("Advect: Integration scheme not found.")
858 .print();
859 }
860 }
861
863 void apply() {
864 if (advectionTime == 0.) {
865 advectedTime = advect(std::numeric_limits<double>::max());
866 numberOfTimeSteps = 1;
867 } else {
868 double currentTime = 0.0;
869 numberOfTimeSteps = 0;
870 while (currentTime < advectionTime) {
871 currentTime += advect(advectionTime - currentTime);
872 ++numberOfTimeSteps;
873 if (performOnlySingleStep)
874 break;
875 }
876 advectedTime = currentTime;
877 }
878 }
879
882 void
883 applyIntegration(double maxTimeStep = std::numeric_limits<double>::max()) {
884 prepareLS();
885 if (integrationScheme == IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER) {
886 auto is = lsInternal::EngquistOsher<T, D, 1>(levelSets.back(), velocities,
887 calculateNormalVectors);
888 currentTimeStep = integrateTime(is, maxTimeStep);
889 } else if (integrationScheme ==
891 auto is = lsInternal::EngquistOsher<T, D, 2>(levelSets.back(), velocities,
892 calculateNormalVectors);
893 currentTimeStep = integrateTime(is, maxTimeStep);
894 } else if (integrationScheme ==
896 auto alphas = findGlobalAlphas();
897 auto is = lsInternal::LaxFriedrichs<T, D, 1>(levelSets.back(), velocities,
898 dissipationAlpha, alphas,
899 calculateNormalVectors);
900 currentTimeStep = integrateTime(is, maxTimeStep);
901 } else if (integrationScheme ==
903 auto alphas = findGlobalAlphas();
904 auto is = lsInternal::LaxFriedrichs<T, D, 2>(levelSets.back(), velocities,
905 dissipationAlpha, alphas,
906 calculateNormalVectors);
907 currentTimeStep = integrateTime(is, maxTimeStep);
908 } else if (integrationScheme ==
912 levelSets.back(), velocities);
913 currentTimeStep = integrateTime(is, maxTimeStep);
914 } else if (integrationScheme ==
917 levelSets.back(), velocities, dissipationAlpha);
918 currentTimeStep = integrateTime(is, maxTimeStep);
919 } else if (integrationScheme ==
922 levelSets.back(), velocities, dissipationAlpha);
923 currentTimeStep = integrateTime(is, maxTimeStep);
924 } else if (integrationScheme ==
927 levelSets.back(), velocities, dissipationAlpha);
928 currentTimeStep = integrateTime(is, maxTimeStep);
929 } else if (integrationScheme ==
932 levelSets.back(), velocities, dissipationAlpha);
933 currentTimeStep = integrateTime(is, maxTimeStep);
934 } else if (integrationScheme ==
937 levelSets.back(), velocities, dissipationAlpha);
938 currentTimeStep = integrateTime(is, maxTimeStep);
939 } else {
940 Logger::getInstance()
941 .addWarning("Advect: Integration scheme not found. Not integrating.")
942 .print();
943 currentTimeStep = -1.;
944 }
945 }
946};
947
948// add all template specializations for this class
950
951} // namespace viennals
Engquist-Osher integration scheme based on the upwind integration scheme. Offers high performance but...
Definition lsEngquistOsher.hpp:18
static void prepareLS(SmartPointer< viennals::Domain< T, D > > passedlsDomain)
Definition lsEngquistOsher.hpp:28
Lax Friedrichs integration scheme with constant alpha value for dissipation. This alpha value should ...
Definition lsLaxFriedrichs.hpp:19
static void prepareLS(SmartPointer< viennals::Domain< T, D > > passedlsDomain)
Definition lsLaxFriedrichs.hpp:32
Lax Friedrichs integration scheme, which uses alpha values provided by the user in getDissipationAlph...
Definition lsLocalLaxFriedrichsAnalytical.hpp:20
static void prepareLS(SmartPointer< viennals::Domain< T, D > > passedlsDomain)
Definition lsLocalLaxFriedrichsAnalytical.hpp:48
Lax Friedrichs integration scheme, which uses a first neighbour stencil to calculate the alpha values...
Definition lsLocalLaxFriedrichs.hpp:20
static void prepareLS(SmartPointer< viennals::Domain< T, D > > passedlsDomain)
Definition lsLocalLaxFriedrichs.hpp:49
Lax Friedrichs integration scheme, which considers only the current point for alpha calculation....
Definition lsLocalLocalLaxFriedrichs.hpp:18
static void prepareLS(SmartPointer< viennals::Domain< T, D > > passedlsDomain)
Definition lsLocalLocalLaxFriedrichs.hpp:29
Stencil Local Lax Friedrichs Integration Scheme. It uses a stencil of order around active points,...
Definition lsStencilLocalLaxFriedrichsScalar.hpp:25
static void prepareLS(LevelSetType passedlsDomain)
Definition lsStencilLocalLaxFriedrichsScalar.hpp:116
void setVelocityField(SmartPointer< VelocityField< T > > passedVelocities)
Set the velocity field used for advection. This should be a concrete implementation of lsVelocityFiel...
Definition lsAdvect.hpp:736
static constexpr char velocityLabel[]
Definition lsAdvect.hpp:707
void insertNextLevelSet(SmartPointer< Domain< T, D > > passedlsDomain)
Pushes the passed level set to the back of the list of level sets used for advection.
Definition lsAdvect.hpp:730
void setUpdatePointData(bool update)
Set whether the point data in the old LS should be translated to the advected LS. Defaults to true.
Definition lsAdvect.hpp:810
bool getCalculateNormalVectors() const
Get whether normal vectors were calculated.
Definition lsAdvect.hpp:791
void setIgnoreVoids(bool iV)
Set whether level set values, which are not part of the "top" geometrically connected part of values,...
Definition lsAdvect.hpp:771
Advect()=default
void setCalculateNormalVectors(bool cnv)
Set whether normal vectors should be calculated at each level set point. Defaults to true....
Definition lsAdvect.hpp:763
void apply()
Perform the advection.
Definition lsAdvect.hpp:863
void setSingleStep(bool singleStep)
If set to true, only a single advection step will be performed, even if the advection time set with s...
Definition lsAdvect.hpp:751
void prepareLS()
Definition lsAdvect.hpp:814
double getAdvectedTime() const
Get by how much the physical time was advanced during the last apply() call.
Definition lsAdvect.hpp:779
void setDissipationAlpha(const double &a)
Set the alpha dissipation coefficient. For lsLaxFriedrichs, this is used as the alpha value....
Definition lsAdvect.hpp:803
void setCheckDissipation(bool check)
Definition lsAdvect.hpp:806
Advect(std::vector< SmartPointer< Domain< T, D > > > passedlsDomains, SmartPointer< VelocityField< T > > passedVelocities)
Definition lsAdvect.hpp:722
void setIntegrationScheme(IntegrationSchemeEnum scheme)
Set which integration scheme should be used out of the ones specified in IntegrationSchemeEnum.
Definition lsAdvect.hpp:795
double getCurrentTimeStep() const
Return the last applied time step.
Definition lsAdvect.hpp:782
void setAdvectionTime(double time)
Set the time until when the level set should be advected. If this takes more than one advection step,...
Definition lsAdvect.hpp:745
unsigned getNumberOfTimeSteps() const
Get how many advection steps were performed during the last apply() call.
Definition lsAdvect.hpp:785
void setTimeStepRatio(const double &cfl)
Set the CFL condition to use during advection. The CFL condition sets the maximum distance a surface ...
Definition lsAdvect.hpp:757
void applyIntegration(double maxTimeStep=std::numeric_limits< double >::max())
This function applies the integration scheme and calculates the rates and the maximum time step,...
Definition lsAdvect.hpp:883
Advect(SmartPointer< Domain< T, D > > passedlsDomain)
Definition lsAdvect.hpp:712
Advect(SmartPointer< Domain< T, D > > passedlsDomain, SmartPointer< VelocityField< T > > passedVelocities)
Definition lsAdvect.hpp:716
double getTimeStepRatio() const
Get the value of the CFL number.
Definition lsAdvect.hpp:788
void setSaveAdvectionVelocities(bool sAV)
Set whether the velocities applied to each point should be saved in the level set for debug purposes.
Definition lsAdvect.hpp:775
static constexpr char dissipationLabel[]
Definition lsAdvect.hpp:708
This class is used to perform boolean operations on two level sets and write the resulting level set ...
Definition lsBooleanOperation.hpp:43
void apply()
Perform operation.
Definition lsBooleanOperation.hpp:316
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
static constexpr T NEG_VALUE
Definition lsDomain.hpp:52
static constexpr T POS_VALUE
Definition lsDomain.hpp:51
This class is used to mark points of the level set which are enclosed in a void.
Definition lsMarkVoidPoints.hpp:28
void apply()
Definition lsMarkVoidPoints.hpp:153
static constexpr char voidPointLabel[]
Definition lsMarkVoidPoints.hpp:85
std::vector< T > ScalarDataType
Definition lsPointData.hpp:22
ScalarDataType * getScalarData(int index)
Definition lsPointData.hpp:130
Reduce the level set size to the specified width. This means all level set points with value <= 0....
Definition lsReduce.hpp:14
void apply()
Reduces the leveleSet to the specified number of layers. The largest value in the levelset is thus wi...
Definition lsReduce.hpp:55
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
Abstract class defining the interface for the velocity field used during advection using lsAdvect.
Definition lsVelocityField.hpp:11
#define PRECOMPILE_PRECISION_DIMENSION(className)
Definition lsPreCompileMacros.hpp:24
Definition lsAdvect.hpp:36
IntegrationSchemeEnum
Enumeration for the different Integration schemes used by the advection kernel.
Definition lsAdvect.hpp:42
@ LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER
Definition lsAdvect.hpp:49
@ STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER
Definition lsAdvect.hpp:52
@ LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER
Definition lsAdvect.hpp:48
@ LAX_FRIEDRICHS_2ND_ORDER
Definition lsAdvect.hpp:46
@ LOCAL_LAX_FRIEDRICHS_1ST_ORDER
Definition lsAdvect.hpp:50
@ ENGQUIST_OSHER_2ND_ORDER
Definition lsAdvect.hpp:44
@ LAX_FRIEDRICHS_1ST_ORDER
Definition lsAdvect.hpp:45
@ LOCAL_LAX_FRIEDRICHS_2ND_ORDER
Definition lsAdvect.hpp:51
@ ENGQUIST_OSHER_1ST_ORDER
Definition lsAdvect.hpp:43
@ LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER
Definition lsAdvect.hpp:47
@ INTERSECT
Definition lsBooleanOperation.hpp:26
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68