ViennaLS
Loading...
Searching...
No Matches
lsLaxFriedrichs.hpp
Go to the documentation of this file.
1#pragma once
2
3#include <hrleSparseStarIterator.hpp>
4#include <hrleTypes.hpp>
5
6#include <lsDomain.hpp>
7#include <lsExpand.hpp>
8
9#include <vcVectorType.hpp>
10
11namespace lsInternal {
12
13using namespace viennacore;
14
19template <class T, int D, int order> class LaxFriedrichs {
20 SmartPointer<viennals::Domain<T, D>> levelSet;
21 SmartPointer<viennals::VelocityField<T>> velocities;
22 viennahrle::SparseStarIterator<viennahrle::Domain<T, D>, order>
23 neighborIterator;
24 const double alphaFactor = 1.0;
25 VectorType<T, 3> const finalAlphas;
26 const bool calculateNormalVectors = true;
27
28 static T pow2(const T &value) { return value * value; }
29
30public:
31 // static const int order_ = order;
32 static void prepareLS(SmartPointer<viennals::Domain<T, D>> passedlsDomain) {
33 assert(order == 1 || order == 2);
34 viennals::Expand<T, D>(passedlsDomain, 2 * order + 1).apply();
35 }
36
37 LaxFriedrichs(SmartPointer<viennals::Domain<T, D>> passedlsDomain,
38 SmartPointer<viennals::VelocityField<T>> vel, double alpha,
39 VectorType<T, 3> &alphas, bool calcNormal)
40 : levelSet(passedlsDomain), velocities(vel),
41 neighborIterator(levelSet->getDomain()), alphaFactor(alpha),
42 finalAlphas(alphas), calculateNormalVectors(calcNormal) {}
43
44 std::pair<T, T> operator()(const viennahrle::Index<D> &indices,
45 int material) {
46
47 auto &grid = levelSet->getGrid();
48 double gridDelta = grid.getGridDelta();
49
50 VectorType<T, 3> coordinate{0., 0., 0.};
51 for (unsigned i = 0; i < D; ++i) {
52 coordinate[i] = indices[i] * gridDelta;
53 }
54
55 // move neighborIterator to current position
56 neighborIterator.goToIndicesSequential(indices);
57
58 T gradPos[D];
59 T gradNeg[D];
60
61 T grad = 0.;
62 T dissipation = 0.;
63
64 Vec3D<T> normalVector = {};
65 T normalModulus = 0;
66
67 for (int i = 0; i < D; i++) { // iterate over dimensions
68
69 const T deltaPos = gridDelta;
70 const T deltaNeg = -gridDelta;
71
72 const T phi0 = neighborIterator.getCenter().getValue();
73 const T phiPos = neighborIterator.getNeighbor(i).getValue();
74 const T phiNeg = neighborIterator.getNeighbor(i + D).getValue();
75
76 T diffPos = (phiPos - phi0) / deltaPos;
77 T diffNeg = (phiNeg - phi0) / deltaNeg;
78
79 if (order == 2) { // if second order time integration scheme is used
80 const T deltaPosPos = 2 * gridDelta;
81 const T deltaNegNeg = -2 * gridDelta;
82
83 const T diff00 =
84 (((deltaNeg * phiPos - deltaPos * phiNeg) / (deltaPos - deltaNeg) +
85 phi0)) /
86 (deltaPos * deltaNeg);
87 const T phiPosPos =
88 neighborIterator.getNeighbor((D * order) + i).getValue();
89 const T phiNegNeg =
90 neighborIterator.getNeighbor((D * order) + D + i).getValue();
91
92 const T diffNegNeg = (((deltaNeg * phiNegNeg - deltaNegNeg * phiNeg) /
93 (deltaNegNeg - deltaNeg) +
94 phi0)) /
95 (deltaNegNeg * deltaNeg);
96 const T diffPosPos = (((deltaPos * phiPosPos - deltaPosPos * phiPos) /
97 (deltaPosPos - deltaPos) +
98 phi0)) /
99 (deltaPosPos * deltaPos);
100
101 if (std::signbit(diff00) == std::signbit(diffPosPos)) {
102 if (std::abs(diffPosPos * deltaPos) < std::abs(diff00 * deltaNeg)) {
103 diffPos -= deltaPos * diffPosPos;
104 } else {
105 diffPos += deltaNeg * diff00;
106 }
107 }
108
109 if (std::signbit(diff00) == std::signbit(diffNegNeg)) {
110 if (std::abs(diffNegNeg * deltaNeg) < std::abs(diff00 * deltaPos)) {
111 diffNeg -= deltaNeg * diffNegNeg;
112 } else {
113 diffNeg += deltaPos * diff00;
114 }
115 }
116 }
117
118 gradPos[i] = diffNeg;
119 gradNeg[i] = diffPos;
120
121 if (calculateNormalVectors) {
122 normalVector[i] = (diffNeg + diffPos) * 0.5;
123 normalModulus += normalVector[i] * normalVector[i];
124 }
125
126 grad += pow2((diffNeg + diffPos) * 0.5);
127 dissipation += alphaFactor * finalAlphas[i] * (diffPos - diffNeg) * 0.5;
128 }
129
130 if (calculateNormalVectors) {
131 normalModulus = std::sqrt(normalModulus);
132 for (unsigned i = 0; i < D; ++i) {
133 normalVector[i] /= normalModulus;
134 }
135 }
136
137 // convert coordinate to std array for interface
138 Vec3D<T> coordArray{coordinate[0], coordinate[1], coordinate[2]};
139
140 double scalarVelocity = velocities->getScalarVelocity(
141 coordArray, material, normalVector,
142 neighborIterator.getCenter().getPointId());
143 Vec3D<T> vectorVelocity = velocities->getVectorVelocity(
144 coordArray, material, normalVector,
145 neighborIterator.getCenter().getPointId());
146
147 T totalGrad = 0.;
148 if (scalarVelocity != 0.) {
149 totalGrad = scalarVelocity * std::sqrt(grad);
150 }
151
152 for (int w = 0; w < D; w++) {
153 if (vectorVelocity[w] > 0.) {
154 totalGrad += vectorVelocity[w] * gradPos[w];
155 } else {
156 totalGrad += vectorVelocity[w] * gradNeg[w];
157 }
158 }
159
160 return {totalGrad, ((totalGrad != 0.) ? dissipation : 0)};
161 }
162
163 void reduceTimeStepHamiltonJacobi(double &MaxTimeStep,
164 double gridDelta) const {
165 constexpr double alpha_maxCFL = 1.0;
166 // second time step test, based on alphas
167
168 double timeStep = 0;
169 for (int i = 0; i < D; ++i) {
170 timeStep += finalAlphas[i] / gridDelta;
171 }
172
173 timeStep = alpha_maxCFL / timeStep;
174 MaxTimeStep = std::min(timeStep, MaxTimeStep);
175 }
176};
177} // namespace lsInternal
static void prepareLS(SmartPointer< viennals::Domain< T, D > > passedlsDomain)
Definition lsLaxFriedrichs.hpp:32
void reduceTimeStepHamiltonJacobi(double &MaxTimeStep, double gridDelta) const
Definition lsLaxFriedrichs.hpp:163
std::pair< T, T > operator()(const viennahrle::Index< D > &indices, int material)
Definition lsLaxFriedrichs.hpp:44
LaxFriedrichs(SmartPointer< viennals::Domain< T, D > > passedlsDomain, SmartPointer< viennals::VelocityField< T > > vel, double alpha, VectorType< T, 3 > &alphas, bool calcNormal)
Definition lsLaxFriedrichs.hpp:37
Class containing all information about the level set, including the dimensions of the domain,...
Definition lsDomain.hpp:27
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
Abstract class defining the interface for the velocity field used during advection using lsAdvect.
Definition lsVelocityField.hpp:11
Definition lsCurvatureFormulas.hpp:9
constexpr int D
Definition pyWrap.cpp:70
double T
Definition pyWrap.cpp:68