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

adapt mmdg to changes in the model

parent 21a7a1a5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 47 additions and 38 deletions
dune/mmdg/mmdg.hh
+ 47
− 38
View file @ ac300083
......@@ -223,7 +223,7 @@ private:
const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn);
const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut);
//evalute the coupling terms
//evaluate 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
......@@ -267,7 +267,7 @@ private:
}
}
//evalute the coupling terms
//evaluate 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
......@@ -329,7 +329,7 @@ private:
}
}
//evalute the coupling terms
//evaluate 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
......@@ -400,8 +400,8 @@ private:
}
//evaluate
// int_iElem (Kparallel * grad(phi_iElem,i)) * grad(phi_iElem,j) ds
// = int_iElem (Kpar * tau_i) * tau_j ds
// int_iElem (Kparallel * grad(d phi_iElem,i)) * grad(d phi_iElem,j) ds
// = int_iElem d^2 * (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)
{
......@@ -410,21 +410,24 @@ private:
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
const auto KEvaluation = Base::problem_.Kparallel(qpGlobal);
Scalar d2Evaluation = Base::problem_.aperture(qpGlobal);
d2Evaluation *= d2Evaluation;
for (int i = 0; i < dim-1 ; i++)
{
Coordinate KparDotTau_i;
Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], KparDotTau_i);
KEvaluation.mv(iFrame[i], KparDotTau_i);
for (int j = 0; j <= i; j++)
{
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
quadratureFactor * ( KparDotTau_i * iFrame[j] );
quadratureFactor * d2Evaluation * ( KparDotTau_i * iFrame[j] );
if (i != j)
{
Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
quadratureFactor * ( KparDotTau_i * iFrame[j] );
quadratureFactor * d2Evaluation * ( KparDotTau_i * iFrame[j] );
}
}
}
......@@ -437,18 +440,21 @@ private:
const auto& center = intersct.geometry().center();
const auto& normal = intersct.centerUnitOuterNormal();
Coordinate KparDotTau_i;
const Scalar dEvaluation = Base::problem_.aperture(center);
const Scalar d2Evaluation = dEvaluation * dEvaluation;
//evaluate the interface term
// int_intersct mu^Gamma * jump(phi_iElem,0) * jump(phi_iElem,0) dr
Base::A->entry(iElemIdxSLE, iElemIdxSLE) += mu;
// int_intersct mu^Gamma * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
Base::A->entry(iElemIdxSLE, iElemIdxSLE) += mu * d2Evaluation;
if (intersct.neighbor()) //inner interface edge
{
//evaluate the interface terms
// int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_iElem,j) dr
// int_intersct mu^Gamma * d^2 * jump(phi_iElem,i) * jump(phi_iElem,j) dr
//and
// -int_intersct avg(Kparallel * grad(phi_iElem,i))
// * jump(phi_iElem,j) dr
// -int_intersct d * avg(Kparallel * grad(d * phi_iElem,i))
// * jump(phi_iElem,j) dr //TODO term with grad(d)
//for i,j = 0,...,dim-1 and (i != 0 or j != 0)
for (int i = 0; i < dim - 1; i++)
{
......@@ -457,7 +463,7 @@ private:
const Scalar interfaceUpdate1 = mu * centerI;
const Scalar interfaceUpdate2 = KparDotTau_i * normal;
const Scalar interfaceUpdate3 =
interfaceUpdate1 - 0.5 * interfaceUpdate2;
d2Evaluation * (interfaceUpdate1 - 0.5 * interfaceUpdate2);
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
interfaceUpdate3;
......@@ -470,10 +476,10 @@ private:
{ //NOTE: only relevant for n=3, requires quadrature rule for n=3
Coordinate KparDotTau_j;
Base::problem_.Kparallel(center).mv(iFrame[j], KparDotTau_j);
const Scalar interfaceUpdate4 =
(center * iFrame[j]) * interfaceUpdate1
const Scalar interfaceUpdate4 = d2Evaluation *
((center * iFrame[j]) * interfaceUpdate1
- 0.5 * (interfaceUpdate2 * (center * iFrame[j])
+ (KparDotTau_j * normal) * centerI);
+ (KparDotTau_j * normal) * centerI));
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate4;
......@@ -494,14 +500,15 @@ private:
(dim + 1) * Base::gridView_.size(0) + dim * neighborIdx;
//evaluate the interface terms
// int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_neighbor,j) dr
// int_intersct mu^Gamma * d^2 * jump(phi_iElem,i)
// * jump(phi_neighbor,j) dr
//and
// -int_intersct avg(Kparallel * grad(phi_iElem,i))
// -int_intersct d * avg(Kparallel * grad(d * phi_iElem,i))
// * jump(phi_neighbor,j) dr
//resp.
// -int_intersct avg(Kparallel * grad(phi_neighbor,i))
// -int_intersct d * avg(Kparallel * grad(d * phi_neighbor,i))
// * jump(phi_iElem,j) dr
//for i,j = 0,...,dim-1
//for i,j = 0,...,dim-1 //TODO terms with grad(d)
Base::A->entry(iElemIdxSLE, neighborIdxSLE) += -mu;
Base::A->entry(neighborIdxSLE, iElemIdxSLE) += -mu;
......@@ -512,9 +519,9 @@ private:
const Scalar interfaceUpdate1 = mu * centerI;
const Scalar interfaceUpdate2 = KparDotTau_i * normal;
const Scalar interfaceUpdate3 =
-interfaceUpdate1 + 0.5 * interfaceUpdate2;
d2Evaluation * (-interfaceUpdate1 + 0.5 * interfaceUpdate2);
const Scalar interfaceUpdate4 =
-interfaceUpdate1 - 0.5 * interfaceUpdate2;
d2Evaluation * (-interfaceUpdate1 - 0.5 * interfaceUpdate2);
Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE) +=
interfaceUpdate3;
......@@ -526,7 +533,7 @@ private:
-interfaceUpdate4;
Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + i + 1) +=
-centerI * interfaceUpdate1;
-d2Evaluation * centerI * interfaceUpdate1;
for (int j = 0; j < i; j++)
{ //NOTE: only relevant for n=3, requires quadrature rule
......@@ -537,9 +544,9 @@ private:
const Scalar interfaceUpdate6 = interfaceUpdate2 *
(center * iFrame[j]) - (KparDotTau_j * normal) * centerI;
const Scalar interfaceUpdate7 =
-interfaceUpdate5 - 0.5 * interfaceUpdate6;
d2Evaluation * (-interfaceUpdate5 - 0.5 * interfaceUpdate6);
const Scalar interfaceUpdate8 =
-interfaceUpdate5 + 0.5 * interfaceUpdate6;
d2Evaluation * (-interfaceUpdate5 + 0.5 * interfaceUpdate6);
Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
interfaceUpdate7;
......@@ -555,18 +562,18 @@ private:
else //boundary interface edge
{
//evaluate the interface terms
// int_intersct mu^Gamma * phi_iElem,i * phi_iElem,j dr
// int_intersct mu^Gamma * d^2 phi_iElem,i * phi_iElem,j dr
//and
// -int_intersct phi_iElem,j * (Kparallel * grad(phi_iElem,i))
// -int_intersct d * phi_iElem,j * (Kparallel * grad(d * phi_iElem,i))
// * normal dr
//for i,j = 0,...,dim-1 and (i != 0 or j != 0)
//and add interface boundary terms to the load vector b:
// int_intersct d * g * (mu * phi_iElem,i +
// (Kparallel * grad(phi_i,iElem)) * normal) dr
//for i = 0,...,dim-1
// int_intersct d * g * (mu * d * phi_iElem,i +
// (Kparallel * grad(d * phi_i,iElem)) * normal) dr
//for i = 0,...,dim-1 //TODO: terms with grad(d)
const Scalar dgGamma = Base::problem_.aperture(center)
* Base::problem_.boundary(center);
Base::b[iElemIdxSLE] += mu * dgGamma;
Base::b[iElemIdxSLE] += mu * d2Evaluation * dgGamma;
for (int i = 0; i < dim - 1; i++)
{
......@@ -575,25 +582,27 @@ private:
const Scalar interfaceUpdate1 = mu * centerI;
const Scalar interfaceUpdate2 = KparDotTau_i * normal;
const Scalar interfaceUpdate3 =
interfaceUpdate1 - interfaceUpdate2;
d2Evaluation * (interfaceUpdate1 - interfaceUpdate2);
Base::b[iElemIdxSLE] += dgGamma * (mu * centerI - interfaceUpdate2);
Base::b[iElemIdxSLE] +=
dgGamma * d2Evaluation * (mu * centerI - interfaceUpdate2);
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
interfaceUpdate3;
Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
interfaceUpdate3;
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
centerI * (interfaceUpdate1 - 2.0 * interfaceUpdate2);
d2Evaluation * centerI *
(interfaceUpdate1 - 2.0 * interfaceUpdate2);
for (int j = 0; j < i; i++)
{ //NOTE: only relevant for n=3, requires quadrature rule for n=3
Coordinate KparDotTau_j;
Base::problem_.Kparallel(center).mv(iFrame[j], KparDotTau_j);
const Scalar interfaceUpdate4 =
(center * iFrame[j]) * interfaceUpdate1
const Scalar interfaceUpdate4 = d2Evaluation *
((center * iFrame[j]) * interfaceUpdate1
- (interfaceUpdate2 * (center * iFrame[j])
+ (KparDotTau_j * normal) * centerI);
+ (KparDotTau_j * normal) * centerI));
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate4;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment