135 std::pair<T, T>
operator()(
const viennahrle::Index<D> &indices,
137 auto &grid = levelSet->getGrid();
138 double gridDelta = grid.getGridDelta();
140 Vec3D<T> coordinate{0., 0., 0.};
141 for (
unsigned i = 0; i <
D; ++i) {
142 coordinate[i] = indices[i] * gridDelta;
146 neighborIterator.goToIndicesSequential(indices);
149 Vec3D<T> coordArray{coordinate[0], coordinate[1], coordinate[2]};
154 Vec3D<T> normalVector;
156 for (
unsigned i = 0; i <
D; i++) {
157 viennahrle::Index<D> neighborIndex(0);
158 neighborIndex[i] = 1;
160 T pos = neighborIterator.getNeighbor(neighborIndex).getValue() -
161 neighborIterator.getCenter().getValue();
162 T neg = neighborIterator.getCenter().getValue() -
163 neighborIterator.getNeighbor(-neighborIndex).getValue();
164 normalVector[i] = (pos + neg) * 0.5;
166 denominator += normalVector[i] * normalVector[i];
168 denominator = std::sqrt(denominator);
170 for (
unsigned i = 0; i <
D; ++i) {
171 normalVector[i] /= denominator;
174 double scalarVelocity = velocities->getScalarVelocity(
175 coordArray, material, normalVector,
176 neighborIterator.getCenter().getPointId());
177 auto vectorVelocity = velocities->getVectorVelocity(
178 coordArray, material, normalVector,
179 neighborIterator.getCenter().getPointId());
182 for (
unsigned i = 0; i <
D; ++i) {
183 scalarVelocity += vectorVelocity[i] * normalVector[i];
186 if (scalarVelocity ==
T(0)) {
191 Norm(calculateGradient(viennahrle::Index<D>(0))) * scalarVelocity;
197 std::vector<VectorType<T, D>> alphas;
198 alphas.reserve(numStencilPoints);
200 viennahrle::Index<D> currentIndex(-order);
201 for (
size_t i = 0; i < numStencilPoints; ++i) {
202 VectorType<T, D> alpha;
203 Vec3D<T> normal = calculateNormal(currentIndex);
206 if ((std::abs(normal[0]) < 1e-6) && (std::abs(normal[1]) < 1e-6) &&
207 (std::abs(normal[2]) < 1e-6)) {
211 Vec3D<T> normal_p{normal[0], normal[1], normal[2]};
212 Vec3D<T> normal_n{normal[0], normal[1], normal[2]};
214 VectorType<T, D> velocityDelta;
215 std::fill(velocityDelta.begin(), velocityDelta.end(), 0.);
218 Vec3D<T> localCoordArray = coordArray;
219 for (
unsigned dir = 0; dir <
D; ++dir)
220 localCoordArray[dir] += currentIndex[dir];
222 T localScalarVelocity = velocities->getScalarVelocity(
223 localCoordArray, material, normal_p,
224 neighborIterator.getCenter().getPointId());
225 Vec3D<T> localVectorVelocity = velocities->getVectorVelocity(
226 localCoordArray, material, normal_p,
227 neighborIterator.getCenter().getPointId());
229 for (
unsigned dir = 0; dir <
D; ++dir) {
230 localScalarVelocity += localVectorVelocity[dir] * normal[dir];
233 const T DN = std::abs(normalEpsilon * scalarVelocity);
235 for (
int k = 0; k <
D; ++k) {
240 T vp = velocities->getScalarVelocity(
241 localCoordArray, material, normal_p,
242 neighborIterator.getCenter().getPointId());
243 T vn = velocities->getScalarVelocity(
244 localCoordArray, material, normal_n,
245 neighborIterator.getCenter().getPointId());
247 velocityDelta[k] = (vn - vp) / (2.0 * DN);
254 for (
int k = 0; k <
D; ++k) {
261 VectorType<T, D> gradient = calculateGradient(currentIndex);
263 for (
int j = 0; j <
D - 1; ++j) {
264 int idx = (k + 1 + j) %
D;
265 monti += gradient[idx] * gradient[idx];
266 toifl += gradient[idx] * velocityDelta[idx];
269 T denom = DotProduct(gradient, gradient);
270 monti *= velocityDelta[k] / denom;
271 toifl *= -gradient[k] / denom;
274 T osher = localScalarVelocity * normal[k];
276 alpha[k] = std::fabs(monti + toifl + osher);
278 if (alpha[k] > maxDissipation) {
283 alphas.push_back(alpha);
288 for (; dim <
D - 1; ++dim) {
289 if (currentIndex[dim] < order)
291 currentIndex[dim] = -order;
298 VectorType<T, D> gradientDiff = calculateGradientDiff();
299 for (
int d = 0; d <
D; ++d) {
302 for (
size_t i = 0; i < numStencilPoints; ++i) {
303 maxAlpha = std::max(maxAlpha, alphas[i][d]);
306 finalAlphas[d] = maxAlpha;
307 dissipation += maxAlpha * gradientDiff[d];
310 return {hamiltonian, dissipation};