From 420ab854e7cc40de95abeff3db6ddb3d84c5b49a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Sun, 15 Mar 2020 15:55:30 +0100
Subject: [PATCH] add coupling terms to mmdg
---
dune/mmdg/mmdg.hh | 214 +++++++++++++++++++++++++++++++++++++++++-----
1 file changed, 192 insertions(+), 22 deletions(-)
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 24a4c4b..044da98 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -105,6 +105,8 @@ private:
//basis functions (0, phi_iElem,i)
const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;
+ // === interface entries ===
+
//fill load vector b with interface entries using a quadrature rule
for (const auto& ip : rule_qInterface)
{
@@ -112,17 +114,17 @@ private:
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
- const Scalar weight(ip.weight());
- const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+ const Scalar quadratureFactor =
+ ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal));
//quadrature for int_elem q*phi_elem,0 dV
- Base::b[iElemIdxSLE] += weight * integrationElem * qInterfaceEvaluation;
+ Base::b[iElemIdxSLE] += quadratureFactor * qInterfaceEvaluation;
//quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim - 1; i++)
{
- Base::b[iElemIdxSLE + i + 1] += weight * integrationElem *
+ Base::b[iElemIdxSLE + i + 1] += quadratureFactor *
(qpGlobal * iFrame[i]) * qInterfaceEvaluation;
}
}
@@ -138,8 +140,8 @@ private:
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
- const Scalar weight(ip.weight());
- const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+ const Scalar quadratureFactor =
+ ip.weight() * iGeo.integrationElement(qpLocal);
for (int i = 0; i < dim-1 ; i++)
{
@@ -149,7 +151,7 @@ private:
for (int j = 0; j <= i; j++) //TODO: symmetrize
{
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
- weight * integrationElem * Base::problem_.aperture(qpGlobal)
+ quadratureFactor * Base::problem_.aperture(qpGlobal)
* ( K_parDotTau_i * iFrame[j] );
}
}
@@ -157,52 +159,220 @@ private:
//TODO: check interface bilinear form
+ // === coupling entries ===
+
//evaluate the coupling term
// int_iElem beta * phi_iElem,i * phi_iElem,j ds
- //using a quadrature rule for i = 0 and j = 0,...,dim-1 //TODO or vice versa (symmetrize)
+ //using a quadrature rule for i = 0 and j = 0,...,dim-1 or vice versa //TODO symmetrize more efficiently
for (const auto& ip : rule_Kperp)
{
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
- const Scalar weight(ip.weight());
- const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+ const Scalar quadratureFactor =
+ ip.weight() * iGeo.integrationElement(qpLocal);
+ const Scalar betaEvaluation = beta(qpGlobal);
Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
- weight * integrationElem * beta(qpGloabal);
+ quadratureFactor * betaEvaluation;
for (int i = 0; i < dim - 1; i++)
{
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
- weight * integrationElem * beta(qpGlobal) * (iFrame[i] * qpGlobal)
- * (iFrame[j] * qpGlobal);
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
+ quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
+ quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
}
}
//evaluate the coupling term
// int_iElem beta * phi_iElem,i * phi_iElem,j ds
- //using a quadrature rule for i,j = 1,...,dim-1 //TODO symmetrize
+ //using a quadrature rule for i,j = 1,...,dim-1 //TODO symmetrize more efficiently
for (const auto& ip : rule_KperpPlus1)
{
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
- const Scalar weight(ip.weight());
- const Scalar integrationElem(iGeo.integrationElement(qpLocal));
+ const Scalar quadratureFactor =
+ ip.weight() * iGeo.integrationElement(qpLocal);
+ const Scalar betaEvaluation = beta(qpGlobal);
for (int i = 0; i < dim - 1; i++)
{
for (int j = 0; j <= i; j++)
{
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
- weight * integrationElem * beta(qpGlobal) * (iFrame[i] * qpGlobal) ;
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
+ * (iFrame[j] * qpGlobal);
+
+ if (i != j)
+ {
+ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
+ * (iFrame[j] * qpGlobal);
+ }
+ }
+ }
+ }
+
+ //get neighboring bulk grid elements of iELem
+ const auto elemIn = intersct.inside();
+ const auto elemOut = intersct.outside();
+
+ const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn);
+ const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut);
+
+ //evalute the coupling terms
+ // int_iElem alpha * jump(phi_elem1,0) * jump(phi_elem2,i) ds
+ //and
+ // int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds
+ //for elem1, elem2 in {elemIn, elemOut} and i = 0,...,dim
+ //using a quadrature rule
+ for (const auto& ip : rule_Kperp)
+ {
+ const auto& qpLocal = ip.position();
+ const auto& qpGlobal = iGeo.global(qpLocal);
+
+ const Scalar quadratureFactor =
+ ip.weight() * iGeo.integrationElement(qpLocal);
+ const Scalar alphaEvaluation = alpha(qpGlobal);
+ const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
+
+ const Scalar bulkUpdate1 =
+ quadratureFactor * (alphaEvaluation + betaEvaluation4);
+ const Scalar bulkUpdate2 =
+ quadratureFactor * (-alphaEvaluation + betaEvaluation4);
+
+ Base::A->entry(elemInIdxSLE, elemInIdxSLE) += bulkUpdate1;
+ Base::A->entry(elemOutIdxSLE, elemOutIdxSLE) += bulkUpdate1;
+
+ Base::A->entry(elemInIdxSLE, elemOutIdxSLE) += bulkUpdate2;
+ Base::A->entry(elemOutIdxSLE, elemInIdxSLE) += bulkUpdate2;
+
+ for (int i = 0; i < dim; i++)
+ {
+ const Scalar bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
+ const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
+
+ Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE) += bulkUpdate3;
+ Base::A->entry(elemInIdxSLE, elemInIdxSLE + i + 1) += bulkUpdate3;
+ Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE) += bulkUpdate3;
+ Base::A->entry(elemOutIdxSLE, elemOutIdxSLE + i + 1) += bulkUpdate3;
+
+ Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE) += bulkUpdate4;
+ Base::A->entry(elemOutIdxSLE, elemInIdxSLE + i + 1) += bulkUpdate4;
+ Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE) += bulkUpdate4;
+ Base::A->entry(elemInIdxSLE, elemOutIdxSLE + i + 1) += bulkUpdate4;
+ }
+ }
+
+ //evalute the coupling terms
+ // int_iElem alpha * jump(phi_elem1,i) * jump(phi_elem2,j) ds
+ //and
+ // int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,j) ds
+ //for elem1, elem2 in {elemIn, elemOut} and i,j = 1,...,dim
+ //using a quadrature rule
+ for (const auto& ip : rule_KperpPlus1)
+ {
+ const auto& qpLocal = ip.position();
+ const auto& qpGlobal = iGeo.global(qpLocal);
+
+ const Scalar quadratureFactor =
+ ip.weight() * iGeo.integrationElement(qpLocal);
+ const Scalar alphaEvaluation = alpha(qpGlobal);
+ const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
+
+ const Scalar bulkUpdate1 =
+ quadratureFactor * (alphaEvaluation + betaEvaluation4);
+ const Scalar bulkUpdate2 =
+ quadratureFactor * (-alphaEvaluation + betaEvaluation4);
+
+ for (int i = 0; i < dim; i++)
+ {
+ const Scalar bulkUpdate3 = qpGlobal[i] * qpGlobal[i] * bulkUpdate1;
+ const Scalar bulkUpdate4 = qpGlobal[i] * qpGlobal[i] * bulkUpdate2;
+
+ Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
+ bulkUpdate3;
+ Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
+ bulkUpdate3;
+
+ Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
+ bulkUpdate4;
+ Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
+ bulkUpdate4;
+
+ for (int j = 0; j < i; j++)
+ {
+ const Scalar bulkUpdate5 = qpGlobal[i] * qpGlobal[j] * bulkUpdate1;
+ const Scalar bulkUpdate6 = qpGlobal[i] * qpGlobal[j] * bulkUpdate2;
+
+ Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
+ bulkUpdate5;
+ Base::A->entry(elemInIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
+ bulkUpdate5;
+ Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
+ bulkUpdate5;
+ Base::A->entry(elemOutIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
+ bulkUpdate5;
+
+ Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
+ bulkUpdate6;
+ Base::A->entry(elemOutIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
+ bulkUpdate6;
+ Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
+ bulkUpdate6;
+ Base::A->entry(elemInIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
+ bulkUpdate6;
}
}
}
- const auto elem1 = intersct.inside();
- const auto elem2 = intersct.ouside();
- //TODO coupling entries
+ //evalute the coupling terms
+ // -int_iElem beta * phi_iElem,i * avg(phi_elem,j) ds
+ //for (i = 0 and j = 0,...,dim-1) or (i = 0,...,dim and j = 0)
+ //and for elem in {elemIn, elemOut} using a quadrature rule
+ for (const auto& ip : rule_Kperp) //TODO: symmetrize more efficiently
+ //TODO: apply quadrature rule only once
+ {
+ 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);
+
+ const Scalar couplingUpdate1 = quadratureFactor * betaEvaluation2;
+
+ Base::A->entry(elemInIdxSLE, iElemIdx) -= couplingUpdate1;
+ Base::A->entry(iElemIdx, elemInIdxSLE) -= couplingUpdate1;
+ Base::A->entry(elemOutIdxSLE, iElemIdx) -= couplingUpdate1;
+ Base::A->entry(iElemIdx, elemOutIdxSLE) -= couplingUpdate1;
+
+ for (int i = 0; i < dim - 1; i++)
+ {
+ const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
+ const Scalar couplingUpdate3 =
+ (iFrame[i] * qpGlobal) * couplingUpdate1;
+
+ Base::A->entry(elemInIdxSLE + i + 1, iElemIdx) -= couplingUpdate2;
+ Base::A->entry(iElemIdx, elemInIdxSLE + i + 1) -= couplingUpdate2;
+ Base::A->entry(elemOutIdxSLE + i + 1, iElemIdx) -= couplingUpdate2;
+ Base::A->entry(iElemIdx, elemOutIdxSLE + i + 1) -= couplingUpdate2;
+
+ Base::A->entry(elemInIdxSLE, iElemIdx + i + 1) -= couplingUpdate3;
+ Base::A->entry(iElemIdx + i + 1, elemInIdxSLE) -= couplingUpdate3;
+ Base::A->entry(elemOutIdxSLE, iElemIdx + i + 1) -= couplingUpdate3;
+ Base::A->entry(iElemIdx + i + 1, elemOutIdxSLE) -= couplingUpdate3;
+ }
+
+ //i = dim - 1
+ const Scalar couplingUpdate4 = qpGlobal[dim] * couplingUpdate1;
+ Base::A->entry(elemInIdxSLE + dim + 1, iElemIdx) -= couplingUpdate4;
+ Base::A->entry(iElemIdx, elemInIdxSLE + dim + 1) -= couplingUpdate4;
+ Base::A->entry(elemOutIdxSLE + dim + 1, iElemIdx) -= couplingUpdate4;
+ Base::A->entry(iElemIdx, elemOutIdxSLE + dim + 1) -= couplingUpdate4;
+ }
}
}
--
GitLab