From 24b01f7e6eead8ead4af2f9cf6ee01ed043e98fb 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 19:53:10 +0100
Subject: [PATCH] [bugfix] add missing coupling entries, [bugfix] fix index
errors for coupling entries
---
dune/mmdg/mmdg.hh | 43 ++++++++++++++++++++++++++++++++++++++-----
1 file changed, 38 insertions(+), 5 deletions(-)
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 03077fd..102d0e2 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -384,11 +384,44 @@ private:
}
//i = dim - 1
- const Scalar couplingUpdate4 = qpGlobal[dim] * couplingUpdate1;
- (*Base::A)[elemInIdxSLE + dim + 1][iElemIdxSLE] -= couplingUpdate4;
- (*Base::A)[iElemIdxSLE][elemInIdxSLE + dim + 1] -= couplingUpdate4;
- (*Base::A)[elemOutIdxSLE + dim + 1][iElemIdxSLE] -= couplingUpdate4;
- (*Base::A)[iElemIdxSLE][elemOutIdxSLE + dim + 1] -= couplingUpdate4;
+ const Scalar couplingUpdate4 = qpGlobal[dim - 1] * couplingUpdate1;
+
+ (*Base::A)[elemInIdxSLE + dim][iElemIdxSLE] -= couplingUpdate4;
+ (*Base::A)[iElemIdxSLE][elemInIdxSLE + dim] -= couplingUpdate4;
+ (*Base::A)[elemOutIdxSLE + dim][iElemIdxSLE] -= couplingUpdate4;
+ (*Base::A)[iElemIdxSLE][elemOutIdxSLE + dim] -= couplingUpdate4;
+ }
+
+ //evalute the coupling terms
+ // -int_iElem beta * phi_iElem,i * avg(phi_elem,j) ds
+ //for i = 1,...,dim-1 and j = 1,...,dim
+ //and for elem in {elemIn, elemOut} using a quadrature rule
+ for (const auto& ip : rule_KperpPlus1) //TODO: symmetrize more efficiently
+ //TODO: apply quadrature rule only once
+ { //TODO TODO TODO
+ const auto& qpLocal = ip.position();
+ const auto& qpGlobal = iGeo.global(qpLocal);
+
+ const Scalar quadratureFactor =
+ ip.weight() * iGeo.integrationElement(qpLocal);
+ const Scalar betaEvaluation2 = 0.5 * beta(qpGlobal);
+
+ for (int i = 0; i < dim; i++)
+ {
+ for (int j = 0; j < dim - 1; j++)
+ {
+ const Scalar couplingUpdate = quadratureFactor * betaEvaluation2
+ * (iFrame[j] * qpGlobal) * qpGlobal[i];
+ (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
+ couplingUpdate;
+ (*Base::A)[iElemIdxSLE + j + 1][elemInIdxSLE + i + 1] -=
+ couplingUpdate;
+ (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
+ couplingUpdate;
+ (*Base::A)[iElemIdxSLE + j + 1][elemOutIdxSLE + i + 1] -=
+ couplingUpdate;
+ }
+ }
}
// === interface entries ===
--
GitLab