diff --git a/dune/mmdg/dg.hh b/dune/mmdg/dg.hh
index 03387bd35f0d6335eb673d0a6304619676be3625..b3372b7d897b46f5449017539fcd379390e5da5a 100644
--- a/dune/mmdg/dg.hh
+++ b/dune/mmdg/dg.hh
@@ -80,7 +80,7 @@ public:
}
//compute L2 error using a quadrature rule
- Scalar computeL2error(const int quadratureOrder = 5) //const
+ Scalar computeL2error(const int quadratureOrder = 5)
{
if (!problem_.hasExactSolution())
{
@@ -97,7 +97,7 @@ public:
const QRuleElem& qrule = QRulesElem::rule(simplex, quadratureOrder);
//determine L2-error on element elem using a quadrature rule
- const auto& qlambda =
+ const auto qlambda =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
Scalar numericalSolutionAtQP = d[elemIdxSLE];
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 489b352db701aada8fe0078a92bffd3c4bd918f8..8b05f05bd75ef8a3c12a45eef422be15fc3dada9 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -48,7 +48,7 @@ public:
//compute L2 error using quadrature rules
// norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) )
- Scalar computeL2error(const int quadratureOrder = 5) const
+ Scalar computeL2error(const int quadratureOrder = 5)
{
const Scalar bulkError = Base::computeL2error(quadratureOrder);
Scalar interfaceError(0.0);
@@ -65,27 +65,24 @@ public:
//determine L2-error on the interface element iElem using a
//quadrature rule
- for (const auto& ip : qrule)
- {
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = iGeo.global(qpLocal);
+ const auto qlambda =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
+ {
+ Scalar numericalSolutionAtQP = Base::d[iElemIdxSLE];
- const Scalar quadratureFactor =
- ip.weight() * iGeo.integrationElement(qpLocal);
+ for (int i = 0; i < dim - 1; i++)
+ {
+ numericalSolutionAtQP +=
+ Base::d[iElemIdxSLE + i + 1] * (qpGlobal * iFrame[i]);
+ }
- Scalar numericalSolutionAtQP = Base::d[iElemIdxSLE];
+ const Scalar errorEvaluation = numericalSolutionAtQP
+ - Base::problem_.exactInterfaceSolution(qpGlobal);
- for (int i = 0; i < dim - 1; i++)
- {
- numericalSolutionAtQP +=
- Base::d[iElemIdxSLE + i + 1] * (qpGlobal * iFrame[i]);
- }
+ interfaceError += weight * errorEvaluation * errorEvaluation;
+ };
- const Scalar errorEvaluation = numericalSolutionAtQP
- - Base::problem_.exactInterfaceSolution(qpGlobal);
-
- interfaceError += quadratureFactor * errorEvaluation * errorEvaluation;
- }
+ Base::applyQuadratureRule(iGeo, qrule, qlambda);
}
//print bulk and interface error
@@ -159,68 +156,9 @@ private:
//in the total system of linear equations (SLE) Ad = b,
//the index iElemIdxSLE refers to the basis function (0, phi_iElem,0)
//and the indices iElemIdxSLE + i + 1, i = 1,...,dim-1, refer to the
- //basis functions (0, phi_iElem,i)interfaceinterface
+ //basis functions (0, phi_iElem,i)
const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;
- // === 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 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 quadratureFactor =
- ip.weight() * iGeo.integrationElement(qpLocal);
- const Scalar betaEvaluation = beta(qpGlobal);
-
- (*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
- quadratureFactor * betaEvaluation;
-
- for (int i = 0; i < dim - 1; i++)
- {
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
- quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
- (*Base::A)[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 more efficiently
- 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 interfaceUpdate1 = quadratureFactor * beta(qpGlobal);
-
- for (int i = 0; i < dim - 1; i++)
- {
- 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++)
- {
- const Scalar interfaceUpdate3 =
- interfaceUpdate2 * (iFrame[j] * qpGlobal);
-
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
- interfaceUpdate3;
- (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
- interfaceUpdate3;
- }
- }
- }
-
//get neighboring bulk grid elements of iELem
const auto elemIn = intersct.inside();
const auto elemOut = intersct.outside();
@@ -228,249 +166,301 @@ private:
const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn);
const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut);
- //evalute the coupling terms
+ // === coupling entries ===
+
+ //lambda to evaluate the coupling terms
+ // int_iElem beta * phi_iElem,0 * phi_iElem,j ds
+ //and
// int_iElem K_perp * 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 KEvaluation = Base::problem_.Kperp(qpGlobal);
- const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
-
- const Scalar bulkUpdate1 =
- quadratureFactor * (KEvaluation + betaEvaluation4);
- const Scalar bulkUpdate2 =
- quadratureFactor * (-KEvaluation + betaEvaluation4);
-
- (*Base::A)[elemInIdxSLE][elemInIdxSLE] += bulkUpdate1;
- (*Base::A)[elemOutIdxSLE][elemOutIdxSLE] += bulkUpdate1;
-
- (*Base::A)[elemInIdxSLE][elemOutIdxSLE] += bulkUpdate2;
- (*Base::A)[elemOutIdxSLE][elemInIdxSLE] += bulkUpdate2;
-
- for (int i = 0; i < dim; i++)
- {
- const Scalar bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
- const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
-
- (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate3;
- (*Base::A)[elemInIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate3;
- (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate3;
- (*Base::A)[elemOutIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate3;
-
- (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate4;
- (*Base::A)[elemOutIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate4;
- (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate4;
- (*Base::A)[elemInIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate4;
- }
- }
-
- //evalute the coupling terms
- // int_iElem K_perp * 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 KEvaluation = Base::problem_.Kperp(qpGlobal);
- const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
-
- const Scalar bulkUpdate1 =
- quadratureFactor * (KEvaluation + betaEvaluation4);
- const Scalar bulkUpdate2 =
- quadratureFactor * (-KEvaluation + betaEvaluation4);
-
- for (int i = 0; i < dim; i++)
+ // -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
+ //and
+ // -int_iElem beta * phi_iElem,j * avg(phi_elem,0) ds
+ //for i = 0,...,dim and j = 0,...,dim-1
+ //and for elem in {elemIn, elemOut} using a quadrature rule
+ //TODO symmetrize more efficiently
+ const auto lambda_Kperp =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
{
- const Scalar bulkUpdate3 = qpGlobal[i] * qpGlobal[i] * bulkUpdate1;
- const Scalar bulkUpdate4 = qpGlobal[i] * qpGlobal[i] * bulkUpdate2;
+ const Scalar betaEvaluation = beta(qpGlobal);
+ const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
- (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
- bulkUpdate3;
- (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
- bulkUpdate3;
+ const Scalar bulkUpdate1 =
+ weight * (KEvaluation + 0.25 * betaEvaluation);
+ const Scalar bulkUpdate2 =
+ weight * (-KEvaluation + 0.25 * betaEvaluation);
+ const Scalar couplingUpdate1 = 0.5 * weight * betaEvaluation;
- (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
- bulkUpdate4;
- (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
- bulkUpdate4;
+ //quadrature for
+ // int_iElem beta * phi_iElem,0 * phi_iElem,0 ds
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE] += weight * betaEvaluation;
- for (int j = 0; j < i; j++)
+ //quadrature for
+ // int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds
+ //and
+ // int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,0) ds
+ //for elem1, elem2 in {elemIn, elemOut}
+ (*Base::A)[elemInIdxSLE][elemInIdxSLE] += bulkUpdate1;
+ (*Base::A)[elemOutIdxSLE][elemOutIdxSLE] += bulkUpdate1;
+ (*Base::A)[elemInIdxSLE][elemOutIdxSLE] += bulkUpdate2;
+ (*Base::A)[elemOutIdxSLE][elemInIdxSLE] += bulkUpdate2;
+
+ //quadrature for
+ // -int_iElem beta * phi_iElem,0 * avg(phi_elem,j) ds
+ //for elem in {elemIn, elemOut}
+ (*Base::A)[elemInIdxSLE][iElemIdxSLE] -= couplingUpdate1;
+ (*Base::A)[iElemIdxSLE][elemInIdxSLE] -= couplingUpdate1;
+ (*Base::A)[elemOutIdxSLE][iElemIdxSLE] -= couplingUpdate1;
+ (*Base::A)[iElemIdxSLE][elemOutIdxSLE] -= couplingUpdate1;
+
+ for (int i = 0; i < dim; i++)
{
- const Scalar bulkUpdate5 = qpGlobal[i] * qpGlobal[j] * bulkUpdate1;
- const Scalar bulkUpdate6 = qpGlobal[i] * qpGlobal[j] * bulkUpdate2;
-
- (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
- bulkUpdate5;
- (*Base::A)[elemInIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
- bulkUpdate5;
- (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
- bulkUpdate5;
- (*Base::A)[elemOutIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
- bulkUpdate5;
-
- (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
- bulkUpdate6;
- (*Base::A)[elemOutIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
- bulkUpdate6;
- (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
- bulkUpdate6;
- (*Base::A)[elemInIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
- bulkUpdate6;
- }
- }
- }
-
- //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 bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
+ const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
+ const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
- const Scalar couplingUpdate1 = quadratureFactor * betaEvaluation2;
+ if (i != dim - 1)
+ {
+ const Scalar interfaceUpdate =
+ weight * betaEvaluation * (iFrame[i] * qpGlobal);
+ const Scalar couplingUpdate3 =
+ (iFrame[i] * qpGlobal) * couplingUpdate1;
+
+ //quadrature for
+ // int_iElem beta * phi_iElem,0 * phi_iElem,i ds
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] += interfaceUpdate;
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] += interfaceUpdate;
+
+ //quadrature for
+ // -int_iElem beta * phi_iElem,i * avg(phi_elem,0) ds
+ //for elem in {elemIn, elemOut} using a quadrature rule
+ (*Base::A)[elemInIdxSLE][iElemIdxSLE + i + 1] -= couplingUpdate3;
+ (*Base::A)[iElemIdxSLE + i + 1][elemInIdxSLE] -= couplingUpdate3;
+ (*Base::A)[elemOutIdxSLE][iElemIdxSLE + i + 1] -= couplingUpdate3;
+ (*Base::A)[iElemIdxSLE + i + 1][elemOutIdxSLE] -= couplingUpdate3;
+ }
- (*Base::A)[elemInIdxSLE][iElemIdxSLE] -= couplingUpdate1;
- (*Base::A)[iElemIdxSLE][elemInIdxSLE] -= couplingUpdate1;
- (*Base::A)[elemOutIdxSLE][iElemIdxSLE] -= couplingUpdate1;
- (*Base::A)[iElemIdxSLE][elemOutIdxSLE] -= couplingUpdate1;
+ //quadrature for
+ // int_iElem K_perp * 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}
+ (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate3;
+ (*Base::A)[elemInIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate3;
+ (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate3;
+ (*Base::A)[elemOutIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate3;
+ (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate4;
+ (*Base::A)[elemOutIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate4;
+ (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate4;
+ (*Base::A)[elemInIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate4;
+
+ //quadrature for
+ // -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
+ //for elem in {elemIn, elemOut}
+ (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE] -= couplingUpdate2;
+ (*Base::A)[iElemIdxSLE][elemInIdxSLE + i + 1] -= couplingUpdate2;
+ (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE] -= couplingUpdate2;
+ (*Base::A)[iElemIdxSLE][elemOutIdxSLE + i + 1] -= couplingUpdate2;
+ }
+ };
- for (int i = 0; i < dim - 1; i++)
+ //lambda to evaluate the coupling terms
+ // int_iElem beta * phi_iElem,i * phi_iElem,j ds
+ //and
+ // int_iElem K_perp * jump(phi_elem1,k) * jump(phi_elem2,l) ds
+ //and
+ // int_iElem beta * avg(phi_elem1,k) * avg(phi_elem2,l) ds
+ //and
+ // -int_iElem beta * phi_iElem,i * avg(phi_elem1,k) ds
+ //for i,j = 1,...,dim-1 and k,l = 1,...,dim
+ //and for elem1,elem2 in {elemIn, elemOut} using a quadrature rule
+ //TODO symmetrize more efficiently
+ const auto lambda_KperpPlus1 =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
{
- const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
- const Scalar couplingUpdate3 =
- (iFrame[i] * qpGlobal) * couplingUpdate1;
-
- (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE] -= couplingUpdate2;
- (*Base::A)[iElemIdxSLE][elemInIdxSLE + i + 1] -= couplingUpdate2;
- (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE] -= couplingUpdate2;
- (*Base::A)[iElemIdxSLE][elemOutIdxSLE + i + 1] -= couplingUpdate2;
-
- (*Base::A)[elemInIdxSLE][iElemIdxSLE + i + 1] -= couplingUpdate3;
- (*Base::A)[iElemIdxSLE + i + 1][elemInIdxSLE] -= couplingUpdate3;
- (*Base::A)[elemOutIdxSLE][iElemIdxSLE + i + 1] -= couplingUpdate3;
- (*Base::A)[iElemIdxSLE + i + 1][elemOutIdxSLE] -= couplingUpdate3;
- }
-
- //i = dim - 1
- 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
- {
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = iGeo.global(qpLocal);
+ const Scalar betaEvaluation = beta(qpGlobal);
+ const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
- const Scalar quadratureFactor =
- ip.weight() * iGeo.integrationElement(qpLocal);
- const Scalar betaEvaluation2 = 0.5 * beta(qpGlobal);
+ const Scalar bulkUpdate1 =
+ weight * (KEvaluation + 0.25 * betaEvaluation);
+ const Scalar bulkUpdate2 =
+ weight * (-KEvaluation + 0.25 * betaEvaluation);
- for (int i = 0; i < dim; i++)
- {
- for (int j = 0; j < dim - 1; j++)
+ for (int i = 0; i < dim; i++)
{
- 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;
+ const Scalar qpI = iFrame[i] * qpGlobal;
+ const Scalar interfaceUpdate1 = weight * betaEvaluation * qpI;
+ const Scalar couplingUpdate1 =
+ 0.5 * weight * betaEvaluation * qpGlobal[i];
+ const Scalar couplingUpdate2 =
+ (iFrame[i] * qpGlobal) * couplingUpdate1;
+ const Scalar bulkUpdate3 = qpGlobal[i] * qpGlobal[i] * bulkUpdate1;
+ const Scalar bulkUpdate4 = qpGlobal[i] * qpGlobal[i] * bulkUpdate2;
+
+ if (i != dim - 1)
+ {
+ //quadrature for
+ // int_iElem beta * (phi_iElem,i)^2 ds
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+ interfaceUpdate1 * qpI;
+
+ //quadrature for
+ // -int_iElem beta * phi_iElem,i * avg(phi_elem,i) ds
+ //for elem in {elemIn, elemOut}
+ (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE + i + 1] -=
+ couplingUpdate2;
+ (*Base::A)[iElemIdxSLE + i + 1][elemInIdxSLE + i + 1] -=
+ couplingUpdate2;
+ (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE + i + 1] -=
+ couplingUpdate2;
+ (*Base::A)[iElemIdxSLE + i + 1][elemOutIdxSLE + i + 1] -=
+ couplingUpdate2;
+ }
+
+ //quadrature for
+ // int_iElem K_perp * jump(phi_elem1,i) * jump(phi_elem2,i) ds
+ //and
+ // int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,i) ds
+ //for elem1, elem2 in {elemIn, elemOut}
+ (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
+ bulkUpdate3;
+ (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
+ bulkUpdate3;
+ (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
+ bulkUpdate4;
+ (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
+ bulkUpdate4;
+
+ for (int j = 0; j < i; j++)
+ {
+ const Scalar couplingUpdate3 =
+ (iFrame[j] * qpGlobal) * couplingUpdate1;
+ const Scalar bulkUpdate5 =
+ qpGlobal[i] * qpGlobal[j] * bulkUpdate1;
+ const Scalar bulkUpdate6 =
+ qpGlobal[i] * qpGlobal[j] * bulkUpdate2;
+
+ if (i != dim - 1)
+ {
+ const Scalar interfaceUpdate2 =
+ interfaceUpdate1 * (iFrame[j] * qpGlobal);
+
+ //quadrature for
+ // int_iElem beta * phi_iElem,i * phi_iElem,j ds
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+ interfaceUpdate2;
+ (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ interfaceUpdate2;
+
+ //quadrature for
+ // -int_iElem beta * phi_iElem,i * avg(phi_elem,j) ds
+ //for elem in {elemIn, elemOut}
+ (*Base::A)[elemInIdxSLE + j + 1][iElemIdxSLE + i + 1] -=
+ couplingUpdate3;
+ (*Base::A)[iElemIdxSLE + i + 1][elemInIdxSLE + j + 1] -=
+ couplingUpdate3;
+ (*Base::A)[elemOutIdxSLE + j + 1][iElemIdxSLE + i + 1] -=
+ couplingUpdate3;
+ (*Base::A)[iElemIdxSLE + i + 1][elemOutIdxSLE + j + 1] -=
+ couplingUpdate3;
+ }
+
+ //quadrature for
+ // int_iElem K_perp * 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}
+ (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
+ bulkUpdate5;
+ (*Base::A)[elemInIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
+ bulkUpdate5;
+ (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
+ bulkUpdate5;
+ (*Base::A)[elemOutIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
+ bulkUpdate5;
+ (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
+ bulkUpdate6;
+ (*Base::A)[elemOutIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
+ bulkUpdate6;
+ (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
+ bulkUpdate6;
+ (*Base::A)[elemInIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
+ bulkUpdate6;
+
+ //quadrature for
+ // -int_iElem beta * phi_iElem,j * avg(phi_elem,i) ds
+ //for elem in {elemIn, elemOut}
+ (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
+ couplingUpdate3;
+ (*Base::A)[iElemIdxSLE + j + 1][elemInIdxSLE + i + 1] -=
+ couplingUpdate3;
+ (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
+ couplingUpdate3;
+ (*Base::A)[iElemIdxSLE + j + 1][elemOutIdxSLE + i + 1] -=
+ couplingUpdate3;
+ }
}
- }
- }
+ };
// === interface entries ===
- //fill load vector b with interface entries using a quadrature rule
- for (const auto& ip : rule_qInterface)
- {
- //quadrature point in local and global coordinates
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = iGeo.global(qpLocal);
-
- const Scalar quadratureFactor =
- ip.weight() * iGeo.integrationElement(qpLocal);
- const Scalar qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal));
+ //lambda to fill the load vector b with interface entries using a
+ //quadrature rule
+ const auto lambda_qInterface =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
+ {
+ const Scalar
+ qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal));
- //quadrature for int_elem q*phi_elem,0 dV
- Base::b[iElemIdxSLE] += quadratureFactor * qInterfaceEvaluation;
+ //quadrature for int_elem q*phi_elem,0 dV
+ Base::b[iElemIdxSLE] += weight * qInterfaceEvaluation;
- //quadrature for int_elem q*phi_elem,i dV
- for (int i = 0; i < dim - 1; i++)
- {
- Base::b[iElemIdxSLE + i + 1] += quadratureFactor *
- (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
- }
- }
+ //quadrature for int_elem q*phi_elem,i dV
+ for (int i = 0; i < dim - 1; i++)
+ {
+ Base::b[iElemIdxSLE + i + 1] +=
+ weight * (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
+ }
+ };
- //evaluate
+ //lambda to evaluate
// int_iElem (Kparallel * grad(d phi_iElem,i)) * grad(d phi_iElem,j) ds
// = int_iElem (Kpar * tau_i) * tau_j ds //TODO: terms with grad(d)
//using a quadrature rule for i,j = 1,...,dim-1
- for (const auto& ip : rule_Kpar)
- {
- const auto& qpLocal = ip.position();
- const auto& qpGlobal = iGeo.global(qpLocal);
-
- const Scalar quadratureFactor =
- ip.weight() * iGeo.integrationElement(qpLocal);
- Scalar dEvaluation2 = Base::problem_.aperture(qpGlobal);
- dEvaluation2 *= dEvaluation2;
-
- for (int i = 0; i < dim-1 ; i++)
+ const auto lambda_Kpar =
+ [&](const Coordinate& qpGlobal, const Scalar weight)
{
- Coordinate KparDotTau_i;
- Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], KparDotTau_i);
-
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
- quadratureFactor * dEvaluation2 * ( KparDotTau_i * iFrame[i] );
+ Scalar dEvaluation2 = Base::problem_.aperture(qpGlobal);
+ dEvaluation2 *= dEvaluation2;
- for (int j = 0; j < i; j++)
+ for (int i = 0; i < dim-1 ; i++)
{
- const Scalar interfaceUpdate =
- quadratureFactor * dEvaluation2 * ( KparDotTau_i * iFrame[j] );
+ Coordinate KparDotTau_i;
+ Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], KparDotTau_i);
+
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+ weight * dEvaluation2 * ( KparDotTau_i * iFrame[i] );
+
+ for (int j = 0; j < i; j++)
+ {
+ const Scalar interfaceUpdate =
+ weight * dEvaluation2 * ( KparDotTau_i * iFrame[j] );
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
- interfaceUpdate;
- (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
- interfaceUpdate;
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+ interfaceUpdate;
+ (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ interfaceUpdate;
+ }
}
- }
- }
+ };
+
+ //apply quadrature rules on iElem
+ Base::applyQuadratureRule(iGeo, rule_Kperp, lambda_Kperp);
+ Base::applyQuadratureRule(iGeo, rule_KperpPlus1, lambda_KperpPlus1);
+ Base::applyQuadratureRule(iGeo, rule_qInterface, lambda_qInterface);
+ Base::applyQuadratureRule(iGeo, rule_Kpar, lambda_Kpar);
//iterate over all intersection with the boundary of iElem
//NOTE: only works for 2d!, otherwise quadrature required
@@ -491,7 +481,7 @@ private:
// * jump(phi_iElem,0) dr
(*Base::A)[iElemIdxSLE][iElemIdxSLE] += mu * dEvaluation2;
- if (intersct.neighbor()) //inner interface edge
+ if (intersct.neighbor()) //interior interface edge
{
//evaluate the interface terms //TODO terms with grad(d)
// int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_iElem,j) dr