From 32f043fcc3644642f5b3f8d9d52b0176e7637f5c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Sun, 22 Mar 2020 14:38:53 +0100
Subject: [PATCH] [bugfix] fix index error in an interface entry of the load
vector
---
dune/mmdg/mmdg.hh | 37 +++++++++++++++++++++----------------
1 file changed, 21 insertions(+), 16 deletions(-)
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 7db35c5..824846f 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -12,7 +12,7 @@ public:
static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1);
using Base = DG<GridView, Mapper, Problem>;
- using Scalar = typename Base::Scalar; //NOTE using Base::Scalar
+ using Scalar = typename Base::Scalar;
using Matrix = typename Base::Matrix;
using Vector = typename Base::Vector;
using Coordinate = Dune::FieldVector<Scalar, dim>;
@@ -203,22 +203,25 @@ private:
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
- const Scalar betaEvaluation = beta(qpGlobal);
+ const Scalar interfaceUpdate1 = quadratureFactor * beta(qpGlobal);
for (int i = 0; i < dim - 1; i++)
{
- for (int j = 0; j <= i; j++)
+ const Scalar qpI = iFrame[i] * qpGlobal;
+ const Scalar interfaceUpdate2 = interfaceUpdate1 * qpI;
+
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+ interfaceUpdate2 * qpI;
+
+ for (int j = 0; j < i; j++)
{
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
- quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
- * (iFrame[j] * qpGlobal);
+ const Scalar interfaceUpdate3 =
+ interfaceUpdate2 * (iFrame[j] * qpGlobal);
- if (i != j)
- {
- (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
- quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
- * (iFrame[j] * qpGlobal);
- }
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+ interfaceUpdate3;
+ (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ interfaceUpdate3;
}
}
}
@@ -571,9 +574,10 @@ private:
// int_intersct d * g * (mu * phi_iElem,i +
// (Kparallel * grad(phi_i,iElem)) * normal) dr
//for i = 0,...,dim-1
- const Scalar dgGamma = Base::problem_.aperture(center)
- * Base::problem_.boundary(center);
- Base::b[iElemIdxSLE] += mu * dgGamma;
+ const Scalar dEvaluation = Base::problem_.aperture(center);
+ const Scalar dgGamma = dEvaluation * Base::problem_.boundary(center);
+
+ Base::b[iElemIdxSLE] += mu * dgGamma * dEvaluation;
for (int i = 0; i < dim - 1; i++)
{
@@ -584,7 +588,8 @@ private:
const Scalar interfaceUpdate3 =
interfaceUpdate1 - interfaceUpdate2;
- Base::b[iElemIdxSLE] += dgGamma * (mu * centerI - interfaceUpdate2);
+ Base::b[iElemIdxSLE + i + 1] +=
+ dgGamma * dEvaluation * (mu * centerI - interfaceUpdate2);
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
interfaceUpdate3;
--
GitLab