44 std::pair<T, T>
operator()(
const viennahrle::Index<D> &indices,
47 auto &grid = levelSet->getGrid();
48 double gridDelta = grid.getGridDelta();
50 VectorType<T, 3> coordinate{0., 0., 0.};
51 for (
unsigned i = 0; i <
D; ++i) {
52 coordinate[i] = indices[i] * gridDelta;
56 neighborIterator.goToIndicesSequential(indices);
64 Vec3D<T> normalVector = {};
67 for (
int i = 0; i <
D; i++) {
69 const T deltaPos = gridDelta;
70 const T deltaNeg = -gridDelta;
72 const T phi0 = neighborIterator.getCenter().getValue();
73 const T phiPos = neighborIterator.getNeighbor(i).getValue();
74 const T phiNeg = neighborIterator.getNeighbor(i +
D).getValue();
76 T diffPos = (phiPos - phi0) / deltaPos;
77 T diffNeg = (phiNeg - phi0) / deltaNeg;
80 const T deltaPosPos = 2 * gridDelta;
81 const T deltaNegNeg = -2 * gridDelta;
84 (((deltaNeg * phiPos - deltaPos * phiNeg) / (deltaPos - deltaNeg) +
86 (deltaPos * deltaNeg);
88 neighborIterator.getNeighbor((
D * order) + i).getValue();
90 neighborIterator.getNeighbor((
D * order) +
D + i).getValue();
92 const T diffNegNeg = (((deltaNeg * phiNegNeg - deltaNegNeg * phiNeg) /
93 (deltaNegNeg - deltaNeg) +
95 (deltaNegNeg * deltaNeg);
96 const T diffPosPos = (((deltaPos * phiPosPos - deltaPosPos * phiPos) /
97 (deltaPosPos - deltaPos) +
99 (deltaPosPos * deltaPos);
101 if (std::signbit(diff00) == std::signbit(diffPosPos)) {
102 if (std::abs(diffPosPos * deltaPos) < std::abs(diff00 * deltaNeg)) {
103 diffPos -= deltaPos * diffPosPos;
105 diffPos += deltaNeg * diff00;
109 if (std::signbit(diff00) == std::signbit(diffNegNeg)) {
110 if (std::abs(diffNegNeg * deltaNeg) < std::abs(diff00 * deltaPos)) {
111 diffNeg -= deltaNeg * diffNegNeg;
113 diffNeg += deltaPos * diff00;
118 gradPos[i] = diffNeg;
119 gradNeg[i] = diffPos;
121 if (calculateNormalVectors) {
122 normalVector[i] = (diffNeg + diffPos) * 0.5;
123 normalModulus += normalVector[i] * normalVector[i];
126 grad += pow2((diffNeg + diffPos) * 0.5);
127 dissipation += alphaFactor * finalAlphas[i] * (diffPos - diffNeg) * 0.5;
130 if (calculateNormalVectors) {
131 normalModulus = std::sqrt(normalModulus);
132 for (
unsigned i = 0; i <
D; ++i) {
133 normalVector[i] /= normalModulus;
138 Vec3D<T> coordArray{coordinate[0], coordinate[1], coordinate[2]};
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());
148 if (scalarVelocity != 0.) {
149 totalGrad = scalarVelocity * std::sqrt(grad);
152 for (
int w = 0; w <
D; w++) {
153 if (vectorVelocity[w] > 0.) {
154 totalGrad += vectorVelocity[w] * gradPos[w];
156 totalGrad += vectorVelocity[w] * gradNeg[w];
160 return {totalGrad, ((totalGrad != 0.) ? dissipation : 0)};