From 121bc4feaa487f6ac3cbeee02edc340321818cdb Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Wed, 7 Aug 2024 09:44:24 +0200
Subject: [PATCH] Fix memory access errors when density does not depend on
function value
If InterpolationDerivatives is told not to compute the derivatives
for the interpolation value, then it sometimes leaves the corresponding
matrix entries uninitialized. That's okay: they shall not be used
anyway.
However, the code still did use them: It multiplied them by zero.
This is okay if the matrix entries are actual values, but they may
also be NaNs. In that case, the NaNs propagate, leading to wrong
results. Fix this by modifying the matrix products to really omit
the uninitialized values.
---
dune/gfe/assemblers/localintegralstiffness.hh | 13 +++++++------
1 file changed, 7 insertions(+), 6 deletions(-)
diff --git a/dune/gfe/assemblers/localintegralstiffness.hh b/dune/gfe/assemblers/localintegralstiffness.hh
index e7ed2cfa..8d7cfd0c 100644
--- a/dune/gfe/assemblers/localintegralstiffness.hh
+++ b/dune/gfe/assemblers/localintegralstiffness.hh
@@ -177,6 +177,11 @@ namespace Dune::GFE
const auto& element = localView.element();
+ // The range of input variables that the density depends on
+ const size_t begin = (localDensity_->dependsOnValue()) ? 0 : TargetSpace::CoordinateType::dimension;
+ const size_t end = (localDensity_->dependsOnDerivative()) ? m : TargetSpace::CoordinateType::dimension;
+
+ // The quadrature rule
int quadOrder = (localFiniteElement.type().isSimplex())
? (localFiniteElement.localBasis().order()-1) * 2
: (localFiniteElement.localBasis().order() * gridDim - 1) * 2;
@@ -237,21 +242,17 @@ namespace Dune::GFE
// Store the Euclidean gradient for the conversion from the Euclidean to the Riemannian Hesse matrix.
for (size_t i = 0; i < nDofs; i++)
for (size_t ii=0; ii<embeddedBlocksize; ++ii)
- for (size_t j = 0; j < m; j++)
+ for (size_t j = begin; j < end; j++)
localEmbeddedGradient[i][ii] += densityGradient[j] * interpolationGradient[j][i*embeddedBlocksize+ii];
for (size_t i = 0; i < nDofs*blocksize; i++)
- for (size_t j = 0; j < m; j++)
+ for (size_t j = begin; j < end; j++)
localGradient[i] += densityGradient[j] * interpolationGradientShort[j][i];
///////////////////////////////////////////////////////////////////////
// Chain rule to construct Hessian
///////////////////////////////////////////////////////////////////////
- // The range of input variables that the density depends on
- const size_t begin = (localDensity_->dependsOnValue()) ? 0 : TargetSpace::CoordinateType::dimension;
- const size_t end = (localDensity_->dependsOnDerivative()) ? m : TargetSpace::CoordinateType::dimension;
-
// tmp = hessianDensity * interpolationGradientShort
Matrix<double> tmp(m,nDofs*blocksize);
tmp = 0;
--
GitLab