230 std::pair<T, T>
operator()(
const viennahrle::Index<D> &indices,
232 auto &grid = levelSet->getGrid();
233 double gridDelta = grid.getGridDelta();
235 Vec3D<T> coordinate{0., 0., 0.};
236 for (
unsigned i = 0; i <
D; ++i) {
237 coordinate[i] = indices[i] * gridDelta;
241 neighborIterator.goToIndicesSequential(indices);
246 Vec3D<T> normalVector;
248 for (
unsigned i = 0; i <
D; i++) {
249 viennahrle::Index<D> neighborIndex(0);
250 neighborIndex[i] = 1;
252 T pos = neighborIterator.getNeighbor(neighborIndex).getValue() -
253 neighborIterator.getCenter().getValue();
254 T neg = neighborIterator.getCenter().getValue() -
255 neighborIterator.getNeighbor(-neighborIndex).getValue();
256 normalVector[i] = (pos + neg) * 0.5;
258 denominator += normalVector[i] * normalVector[i];
260 denominator = 1 / std::sqrt(denominator);
261 for (
unsigned i = 0; i <
D; ++i) {
262 normalVector[i] *= denominator;
265 double scalarVelocity = velocities->getScalarVelocity(
266 coordinate, material, normalVector,
267 neighborIterator.getCenter().getPointId());
268 auto vectorVelocity = velocities->getVectorVelocity(
269 coordinate, material, normalVector,
270 neighborIterator.getCenter().getPointId());
273 for (
unsigned i = 0; i <
D; ++i) {
274 scalarVelocity += vectorVelocity[i] * normalVector[i];
277 if (scalarVelocity ==
T(0)) {
282 Norm(calculateGradient(viennahrle::Index<D>(0))) * scalarVelocity;
288 std::array<VectorType<T, D>, numStencilPoints> alphas{};
290 viennahrle::Index<D> currentIndex(-order);
291 for (
size_t i = 0; i < numStencilPoints; ++i) {
292 VectorType<T, D> alpha;
293 Vec3D<T> normal = calculateNormal(currentIndex);
296 if ((std::abs(normal[0]) < 1e-6) && (std::abs(normal[1]) < 1e-6) &&
297 (std::abs(normal[2]) < 1e-6)) {
302 Vec3D<T> normal_p{normal[0], normal[1], normal[2]};
303 Vec3D<T> normal_n{normal[0], normal[1], normal[2]};
305 VectorType<T, D> velocityDelta;
306 velocityDelta.fill(0.);
309 Vec3D<T> localCoordArray = coordinate;
310 for (
unsigned dir = 0; dir <
D; ++dir)
311 localCoordArray[dir] += currentIndex[dir];
313 T localScalarVelocity = velocities->getScalarVelocity(
314 localCoordArray, material, normal_p,
315 neighborIterator.getCenter().getPointId());
316 Vec3D<T> localVectorVelocity = velocities->getVectorVelocity(
317 localCoordArray, material, normal_p,
318 neighborIterator.getCenter().getPointId());
320 for (
unsigned dir = 0; dir <
D; ++dir) {
321 localScalarVelocity += localVectorVelocity[dir] * normal[dir];
326#if ADAPTIVE_EPSILON_STRATEGY == 0
328 DN = std::abs(baseNormalEpsilon * scalarVelocity);
329#elif ADAPTIVE_EPSILON_STRATEGY == 1
332 levelSet->getGrid().getGridDelta());
333#elif ADAPTIVE_EPSILON_STRATEGY == 2
335 DN = calculateAdaptiveEpsilon(localScalarVelocity,
336 levelSet->getGrid().getGridDelta(),
337 localCoordArray, material, normal);
340 DN = std::abs(baseNormalEpsilon * scalarVelocity);
343 for (
int k = 0; k <
D; ++k) {
348 T vp = velocities->getScalarVelocity(
349 localCoordArray, material, normal_p,
350 neighborIterator.getCenter().getPointId());
351 T vn = velocities->getScalarVelocity(
352 localCoordArray, material, normal_n,
353 neighborIterator.getCenter().getPointId());
355 velocityDelta[k] = (vn - vp) / (2.0 * DN);
362 for (
int k = 0; k <
D; ++k) {
369 VectorType<T, D> gradient = calculateGradient(currentIndex);
371 for (
int j = 0; j <
D - 1; ++j) {
372 int idx = (k + 1 + j) %
D;
373 monti += gradient[idx] * gradient[idx];
374 toifl += gradient[idx] * velocityDelta[idx];
377 T denom = DotProduct(gradient, gradient);
378 monti *= velocityDelta[k] / denom;
379 toifl *= -gradient[k] / denom;
382 T osher = localScalarVelocity * normal[k];
384 alpha[k] = std::fabs(monti + toifl + osher);
386 if (alpha[k] > maxDissipation) {
387 alpha[k] = maxDissipation;
396 for (; dim <
D - 1; ++dim) {
397 if (currentIndex[dim] < order)
399 currentIndex[dim] = -order;
406 VectorType<T, D> gradientDiff = calculateGradientDiff();
407 for (
int d = 0; d <
D; ++d) {
410 for (
size_t i = 0; i < numStencilPoints; ++i) {
411 maxAlpha = std::max(maxAlpha, alphas[i][d]);
414 finalAlphas[d] = maxAlpha;
415 dissipation += maxAlpha * gradientDiff[d];
418 return {hamiltonian, dissipation};