132 std::pair<T, T>
operator()(
const viennahrle::Index<D> &indices,
134 auto &grid = levelSet->getGrid();
135 double gridDelta = grid.getGridDelta();
137 Vec3D<T> coordinate{0., 0., 0.};
138 for (
unsigned i = 0; i <
D; ++i) {
139 coordinate[i] = indices[i] * gridDelta;
143 neighborIterator.goToIndicesSequential(indices);
146 Vec3D<T> coordArray{coordinate[0], coordinate[1], coordinate[2]};
151 Vec3D<T> normalVector;
153 for (
unsigned i = 0; i <
D; i++) {
154 viennahrle::Index<D> neighborIndex(0);
155 neighborIndex[i] = 1;
157 T pos = neighborIterator.getNeighbor(neighborIndex).getValue() -
158 neighborIterator.getCenter().getValue();
159 T neg = neighborIterator.getCenter().getValue() -
160 neighborIterator.getNeighbor(-neighborIndex).getValue();
161 normalVector[i] = (pos + neg) * 0.5;
163 denominator += normalVector[i] * normalVector[i];
165 denominator = std::sqrt(denominator);
167 for (
unsigned i = 0; i <
D; ++i) {
168 normalVector[i] /= denominator;
171 double scalarVelocity = velocities->getScalarVelocity(
172 coordArray, material, normalVector,
173 neighborIterator.getCenter().getPointId());
174 auto vectorVelocity = velocities->getVectorVelocity(
175 coordArray, material, normalVector,
176 neighborIterator.getCenter().getPointId());
179 for (
unsigned i = 0; i <
D; ++i) {
180 scalarVelocity += vectorVelocity[i] * normalVector[i];
183 if (scalarVelocity ==
T(0)) {
188 Norm(calculateGradient(viennahrle::Index<D>(0))) * scalarVelocity;
194 std::vector<VectorType<T, D>> alphas;
195 alphas.reserve(numStencilPoints);
197 viennahrle::Index<D> currentIndex(-order);
198 for (
size_t i = 0; i < numStencilPoints; ++i) {
199 VectorType<T, D> alpha;
200 Vec3D<T> normal = calculateNormal(currentIndex);
203 if ((std::abs(normal[0]) < 1e-6) && (std::abs(normal[1]) < 1e-6) &&
204 (std::abs(normal[2]) < 1e-6)) {
208 Vec3D<T> normal_p{normal[0], normal[1], normal[2]};
209 Vec3D<T> normal_n{normal[0], normal[1], normal[2]};
211 VectorType<T, D> velocityDelta;
212 std::fill(velocityDelta.begin(), velocityDelta.end(), 0.);
215 Vec3D<T> localCoordArray = coordArray;
216 for (
unsigned dir = 0; dir <
D; ++dir)
217 localCoordArray[dir] += currentIndex[dir];
219 T localScalarVelocity = velocities->getScalarVelocity(
220 localCoordArray, material, normal_p,
221 neighborIterator.getCenter().getPointId());
222 Vec3D<T> localVectorVelocity = velocities->getVectorVelocity(
223 localCoordArray, material, normal_p,
224 neighborIterator.getCenter().getPointId());
226 for (
unsigned dir = 0; dir <
D; ++dir) {
227 localScalarVelocity += localVectorVelocity[dir] * normal[dir];
230 const T DN = std::abs(normalEpsilon * scalarVelocity);
232 for (
int k = 0; k <
D; ++k) {
237 T vp = velocities->getScalarVelocity(
238 localCoordArray, material, normal_p,
239 neighborIterator.getCenter().getPointId());
240 T vn = velocities->getScalarVelocity(
241 localCoordArray, material, normal_n,
242 neighborIterator.getCenter().getPointId());
244 velocityDelta[k] = (vn - vp) / (2.0 * DN);
251 for (
int k = 0; k <
D; ++k) {
258 VectorType<T, D> gradient = calculateGradient(currentIndex);
260 for (
int j = 0; j <
D - 1; ++j) {
261 int idx = (k + 1 + j) %
D;
262 monti += gradient[idx] * gradient[idx];
263 toifl += gradient[idx] * velocityDelta[idx];
266 T denom = DotProduct(gradient, gradient);
267 monti *= velocityDelta[k] / denom;
268 toifl *= -gradient[k] / denom;
271 T osher = localScalarVelocity * normal[k];
273 alpha[k] = std::fabs(monti + toifl + osher);
276 alphas.push_back(alpha);
281 for (; dim <
D - 1; ++dim) {
282 if (currentIndex[dim] < order)
284 currentIndex[dim] = -order;
291 VectorType<T, D> gradientDiff = calculateGradientDiff();
292 for (
int d = 0; d <
D; ++d) {
295 for (
size_t i = 0; i < numStencilPoints; ++i) {
296 maxAlpha = std::max(maxAlpha, alphas[i][d]);
299 finalAlphas[d] = maxAlpha;
300 dissipation += maxAlpha * gradientDiff[d];
303 return {hamiltonian, dissipation};