Skip to content
Snippets Groups Projects
Commit 24d2a897 authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

use method applyQuadratureRule in mmdg

parent 5ebc1b57
Branches
Tags
No related merge requests found
...@@ -80,7 +80,7 @@ public: ...@@ -80,7 +80,7 @@ public:
} }
//compute L2 error using a quadrature rule //compute L2 error using a quadrature rule
Scalar computeL2error(const int quadratureOrder = 5) //const Scalar computeL2error(const int quadratureOrder = 5)
{ {
if (!problem_.hasExactSolution()) if (!problem_.hasExactSolution())
{ {
...@@ -97,7 +97,7 @@ public: ...@@ -97,7 +97,7 @@ public:
const QRuleElem& qrule = QRulesElem::rule(simplex, quadratureOrder); const QRuleElem& qrule = QRulesElem::rule(simplex, quadratureOrder);
//determine L2-error on element elem using a quadrature rule //determine L2-error on element elem using a quadrature rule
const auto& qlambda = const auto qlambda =
[&](const Coordinate& qpGlobal, const Scalar weight) [&](const Coordinate& qpGlobal, const Scalar weight)
{ {
Scalar numericalSolutionAtQP = d[elemIdxSLE]; Scalar numericalSolutionAtQP = d[elemIdxSLE];
......
...@@ -48,7 +48,7 @@ public: ...@@ -48,7 +48,7 @@ public:
//compute L2 error using quadrature rules //compute L2 error using quadrature rules
// norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) ) // 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); const Scalar bulkError = Base::computeL2error(quadratureOrder);
Scalar interfaceError(0.0); Scalar interfaceError(0.0);
...@@ -65,14 +65,9 @@ public: ...@@ -65,14 +65,9 @@ public:
//determine L2-error on the interface element iElem using a //determine L2-error on the interface element iElem using a
//quadrature rule //quadrature rule
for (const auto& ip : qrule) const auto qlambda =
[&](const Coordinate& qpGlobal, const Scalar weight)
{ {
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
Scalar numericalSolutionAtQP = Base::d[iElemIdxSLE]; Scalar numericalSolutionAtQP = Base::d[iElemIdxSLE];
for (int i = 0; i < dim - 1; i++) for (int i = 0; i < dim - 1; i++)
...@@ -84,8 +79,10 @@ public: ...@@ -84,8 +79,10 @@ public:
const Scalar errorEvaluation = numericalSolutionAtQP const Scalar errorEvaluation = numericalSolutionAtQP
- Base::problem_.exactInterfaceSolution(qpGlobal); - Base::problem_.exactInterfaceSolution(qpGlobal);
interfaceError += quadratureFactor * errorEvaluation * errorEvaluation; interfaceError += weight * errorEvaluation * errorEvaluation;
} };
Base::applyQuadratureRule(iGeo, qrule, qlambda);
} }
//print bulk and interface error //print bulk and interface error
...@@ -159,68 +156,9 @@ private: ...@@ -159,68 +156,9 @@ private:
//in the total system of linear equations (SLE) Ad = b, //in the total system of linear equations (SLE) Ad = b,
//the index iElemIdxSLE refers to the basis function (0, phi_iElem,0) //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 //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; 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 //get neighboring bulk grid elements of iELem
const auto elemIn = intersct.inside(); const auto elemIn = intersct.inside();
const auto elemOut = intersct.outside(); const auto elemOut = intersct.outside();
...@@ -228,81 +166,168 @@ private: ...@@ -228,81 +166,168 @@ private:
const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn); const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn);
const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut); 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 // int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds
//and //and
// int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds // int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds
//for elem1, elem2 in {elemIn, elemOut} and i = 0,...,dim //and
//using a quadrature rule // -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
for (const auto& ip : rule_Kperp) //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 auto& qpLocal = ip.position(); const Scalar betaEvaluation = beta(qpGlobal);
const auto& qpGlobal = iGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal); const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
const Scalar bulkUpdate1 = const Scalar bulkUpdate1 =
quadratureFactor * (KEvaluation + betaEvaluation4); weight * (KEvaluation + 0.25 * betaEvaluation);
const Scalar bulkUpdate2 = const Scalar bulkUpdate2 =
quadratureFactor * (-KEvaluation + betaEvaluation4); weight * (-KEvaluation + 0.25 * betaEvaluation);
const Scalar couplingUpdate1 = 0.5 * weight * betaEvaluation;
//quadrature for
// int_iElem beta * phi_iElem,0 * phi_iElem,0 ds
(*Base::A)[iElemIdxSLE][iElemIdxSLE] += weight * betaEvaluation;
//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)[elemInIdxSLE][elemInIdxSLE] += bulkUpdate1;
(*Base::A)[elemOutIdxSLE][elemOutIdxSLE] += bulkUpdate1; (*Base::A)[elemOutIdxSLE][elemOutIdxSLE] += bulkUpdate1;
(*Base::A)[elemInIdxSLE][elemOutIdxSLE] += bulkUpdate2; (*Base::A)[elemInIdxSLE][elemOutIdxSLE] += bulkUpdate2;
(*Base::A)[elemOutIdxSLE][elemInIdxSLE] += 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++) for (int i = 0; i < dim; i++)
{ {
const Scalar bulkUpdate3 = qpGlobal[i] * bulkUpdate1; const Scalar bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2; const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
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;
}
//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 + i + 1][elemInIdxSLE] += bulkUpdate3;
(*Base::A)[elemInIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate3; (*Base::A)[elemInIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate3;
(*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate3; (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate3;
(*Base::A)[elemOutIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate3; (*Base::A)[elemOutIdxSLE][elemOutIdxSLE + i + 1] += bulkUpdate3;
(*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate4; (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE] += bulkUpdate4;
(*Base::A)[elemOutIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate4; (*Base::A)[elemOutIdxSLE][elemInIdxSLE + i + 1] += bulkUpdate4;
(*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate4; (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE] += bulkUpdate4;
(*Base::A)[elemInIdxSLE][elemOutIdxSLE + i + 1] += 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;
} }
} };
//evalute the coupling terms //lambda to evaluate the coupling terms
// int_iElem K_perp * jump(phi_elem1,i) * jump(phi_elem2,j) ds // int_iElem beta * phi_iElem,i * phi_iElem,j ds
//and //and
// int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,j) ds // int_iElem K_perp * jump(phi_elem1,k) * jump(phi_elem2,l) ds
//for elem1, elem2 in {elemIn, elemOut} and i,j = 1,...,dim //and
//using a quadrature rule // int_iElem beta * avg(phi_elem1,k) * avg(phi_elem2,l) ds
for (const auto& ip : rule_KperpPlus1) //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 auto& qpLocal = ip.position(); const Scalar betaEvaluation = beta(qpGlobal);
const auto& qpGlobal = iGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal); const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
const Scalar bulkUpdate1 = const Scalar bulkUpdate1 =
quadratureFactor * (KEvaluation + betaEvaluation4); weight * (KEvaluation + 0.25 * betaEvaluation);
const Scalar bulkUpdate2 = const Scalar bulkUpdate2 =
quadratureFactor * (-KEvaluation + betaEvaluation4); weight * (-KEvaluation + 0.25 * betaEvaluation);
for (int i = 0; i < dim; i++) for (int i = 0; i < dim; i++)
{ {
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 bulkUpdate3 = qpGlobal[i] * qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * qpGlobal[i] * bulkUpdate2; 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] += (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
bulkUpdate3; bulkUpdate3;
(*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + i + 1] += (*Base::A)[elemOutIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
bulkUpdate3; bulkUpdate3;
(*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + i + 1] += (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + i + 1] +=
bulkUpdate4; bulkUpdate4;
(*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + i + 1] += (*Base::A)[elemOutIdxSLE + i + 1][elemInIdxSLE + i + 1] +=
...@@ -310,9 +335,43 @@ private: ...@@ -310,9 +335,43 @@ private:
for (int j = 0; j < i; j++) for (int j = 0; j < i; j++)
{ {
const Scalar bulkUpdate5 = qpGlobal[i] * qpGlobal[j] * bulkUpdate1; const Scalar couplingUpdate3 =
const Scalar bulkUpdate6 = qpGlobal[i] * qpGlobal[j] * bulkUpdate2; (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] += (*Base::A)[elemInIdxSLE + i + 1][elemInIdxSLE + j + 1] +=
bulkUpdate5; bulkUpdate5;
(*Base::A)[elemInIdxSLE + j + 1][elemInIdxSLE + i + 1] += (*Base::A)[elemInIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
...@@ -321,7 +380,6 @@ private: ...@@ -321,7 +380,6 @@ private:
bulkUpdate5; bulkUpdate5;
(*Base::A)[elemOutIdxSLE + j + 1][elemOutIdxSLE + i + 1] += (*Base::A)[elemOutIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
bulkUpdate5; bulkUpdate5;
(*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + j + 1] += (*Base::A)[elemInIdxSLE + i + 1][elemOutIdxSLE + j + 1] +=
bulkUpdate6; bulkUpdate6;
(*Base::A)[elemOutIdxSLE + j + 1][elemInIdxSLE + i + 1] += (*Base::A)[elemOutIdxSLE + j + 1][elemInIdxSLE + i + 1] +=
...@@ -330,124 +388,50 @@ private: ...@@ -330,124 +388,50 @@ private:
bulkUpdate6; bulkUpdate6;
(*Base::A)[elemInIdxSLE + j + 1][elemOutIdxSLE + i + 1] += (*Base::A)[elemInIdxSLE + j + 1][elemOutIdxSLE + i + 1] +=
bulkUpdate6; 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 couplingUpdate1 = quadratureFactor * betaEvaluation2;
(*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 - 1; i++) //quadrature for
{ // -int_iElem beta * phi_iElem,j * avg(phi_elem,i) ds
const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1; //for elem in {elemIn, elemOut}
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 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] -= (*Base::A)[elemInIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
couplingUpdate; couplingUpdate3;
(*Base::A)[iElemIdxSLE + j + 1][elemInIdxSLE + i + 1] -= (*Base::A)[iElemIdxSLE + j + 1][elemInIdxSLE + i + 1] -=
couplingUpdate; couplingUpdate3;
(*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE + j + 1] -= (*Base::A)[elemOutIdxSLE + i + 1][iElemIdxSLE + j + 1] -=
couplingUpdate; couplingUpdate3;
(*Base::A)[iElemIdxSLE + j + 1][elemOutIdxSLE + i + 1] -= (*Base::A)[iElemIdxSLE + j + 1][elemOutIdxSLE + i + 1] -=
couplingUpdate; couplingUpdate3;
}
} }
} }
};
// === interface entries === // === interface entries ===
//fill load vector b with interface entries using a quadrature rule //lambda to fill the load vector b with interface entries using a
for (const auto& ip : rule_qInterface) //quadrature rule
const auto lambda_qInterface =
[&](const Coordinate& qpGlobal, const Scalar weight)
{ {
//quadrature point in local and global coordinates const Scalar
const auto& qpLocal = ip.position(); qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal));
const auto& qpGlobal = iGeo.global(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 //quadrature for int_elem q*phi_elem,0 dV
Base::b[iElemIdxSLE] += quadratureFactor * qInterfaceEvaluation; Base::b[iElemIdxSLE] += weight * qInterfaceEvaluation;
//quadrature for int_elem q*phi_elem,i dV //quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim - 1; i++) for (int i = 0; i < dim - 1; i++)
{ {
Base::b[iElemIdxSLE + i + 1] += quadratureFactor * Base::b[iElemIdxSLE + i + 1] +=
(qpGlobal * iFrame[i]) * qInterfaceEvaluation; 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 (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) // = int_iElem (Kpar * tau_i) * tau_j ds //TODO: terms with grad(d)
//using a quadrature rule for i,j = 1,...,dim-1 //using a quadrature rule for i,j = 1,...,dim-1
for (const auto& ip : rule_Kpar) const auto lambda_Kpar =
[&](const Coordinate& qpGlobal, const Scalar weight)
{ {
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); Scalar dEvaluation2 = Base::problem_.aperture(qpGlobal);
dEvaluation2 *= dEvaluation2; dEvaluation2 *= dEvaluation2;
...@@ -457,12 +441,12 @@ private: ...@@ -457,12 +441,12 @@ private:
Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], KparDotTau_i); Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], KparDotTau_i);
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] += (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
quadratureFactor * dEvaluation2 * ( KparDotTau_i * iFrame[i] ); weight * dEvaluation2 * ( KparDotTau_i * iFrame[i] );
for (int j = 0; j < i; j++) for (int j = 0; j < i; j++)
{ {
const Scalar interfaceUpdate = const Scalar interfaceUpdate =
quadratureFactor * dEvaluation2 * ( KparDotTau_i * iFrame[j] ); weight * dEvaluation2 * ( KparDotTau_i * iFrame[j] );
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] += (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
interfaceUpdate; interfaceUpdate;
...@@ -470,7 +454,13 @@ private: ...@@ -470,7 +454,13 @@ private:
interfaceUpdate; 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 //iterate over all intersection with the boundary of iElem
//NOTE: only works for 2d!, otherwise quadrature required //NOTE: only works for 2d!, otherwise quadrature required
...@@ -491,7 +481,7 @@ private: ...@@ -491,7 +481,7 @@ private:
// * jump(phi_iElem,0) dr // * jump(phi_iElem,0) dr
(*Base::A)[iElemIdxSLE][iElemIdxSLE] += mu * dEvaluation2; (*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) //evaluate the interface terms //TODO terms with grad(d)
// int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_iElem,j) dr // int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_iElem,j) dr
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment