diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 824846fbba385fede657cd677757a457bf5b05ef..03077fd37829e3aa99894dafcd820aad78787c9a 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -156,6 +156,12 @@ private:
const auto& iGeo = iElem.geometry();
const auto& intersct = Base::gridView_.grid().asIntersection(iElem);
+ if ( Base::mapper_.index(intersct.inside())
+ > Base::mapper_.index(intersct.outside()) )
+ {
+ std::cout << "INDEX ERROR!\n\n";
+ }
+
//get basis of the interface coordinate system for evaluation of
//basis function phi_iElem,i with i=1,...,dim-1
const InterfaceBasis iFrame =
@@ -410,8 +416,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 (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)
{
@@ -420,22 +426,26 @@ private:
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
+ Scalar dEvaluation2 = Base::problem_.aperture(qpGlobal);
+ dEvaluation2 *= dEvaluation2;
for (int i = 0; i < dim-1 ; i++)
{
Coordinate KparDotTau_i;
Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], KparDotTau_i);
- for (int j = 0; j <= i; j++)
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
+ quadratureFactor * dEvaluation2 * ( KparDotTau_i * iFrame[i] );
+
+ for (int j = 0; j < i; j++)
{
- (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
- quadratureFactor * ( KparDotTau_i * iFrame[j] );
+ const Scalar interfaceUpdate =
+ quadratureFactor * dEvaluation2 * ( KparDotTau_i * iFrame[j] );
- if (i != j)
- {
- (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
- quadratureFactor * ( KparDotTau_i * iFrame[j] );
- }
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
+ interfaceUpdate;
+ (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
+ interfaceUpdate;
}
}
}
@@ -448,13 +458,17 @@ private:
const auto& normal = intersct.centerUnitOuterNormal();
Coordinate KparDotTau_i;
+ const Scalar dEvaluation = Base::problem_.aperture(center);
+ const Scalar dEvaluation2 = dEvaluation * dEvaluation;
+
//evaluate the interface term
- // int_intersct mu^Gamma * jump(phi_iElem,0) * jump(phi_iElem,0) dr
- (*Base::A)[iElemIdxSLE][iElemIdxSLE] += mu;
+ // int_intersct mu^Gamma * d^2 * jump(phi_iElem,0)
+ // * jump(phi_iElem,0) dr
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE] += mu * dEvaluation2;
if (intersct.neighbor()) //inner interface edge
{
- //evaluate the interface terms
+ //evaluate the interface terms //TODO terms with grad(d)
// int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_iElem,j) dr
//and
// -int_intersct avg(Kparallel * grad(phi_iElem,i))
@@ -467,23 +481,23 @@ private:
const Scalar interfaceUpdate1 = mu * centerI;
const Scalar interfaceUpdate2 = KparDotTau_i * normal;
const Scalar interfaceUpdate3 =
- interfaceUpdate1 - 0.5 * interfaceUpdate2;
+ dEvaluation2 * (interfaceUpdate1 - 0.5 * interfaceUpdate2);
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
interfaceUpdate3;
(*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
interfaceUpdate3;
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
- centerI * (interfaceUpdate1 - interfaceUpdate2);
+ dEvaluation2 * centerI * (interfaceUpdate1 - 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 = dEvaluation2 *
+ ((center * iFrame[j]) * interfaceUpdate1
- 0.5 * (interfaceUpdate2 * (center * iFrame[j])
- + (KparDotTau_j * normal) * centerI);
+ + (KparDotTau_j * normal) * centerI));
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
interfaceUpdate4;
@@ -503,7 +517,7 @@ private:
const int neighborIdxSLE =
(dim + 1) * Base::gridView_.size(0) + dim * neighborIdx;
- //evaluate the interface terms
+ //evaluate the interface terms //TODO: terms with grad(d)
// int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_neighbor,j) dr
//and
// -int_intersct avg(Kparallel * grad(phi_iElem,i))
@@ -512,8 +526,8 @@ private:
// -int_intersct avg(Kparallel * grad(phi_neighbor,i))
// * jump(phi_iElem,j) dr
//for i,j = 0,...,dim-1
- (*Base::A)[iElemIdxSLE][neighborIdxSLE] += -mu;
- (*Base::A)[neighborIdxSLE][iElemIdxSLE] += -mu;
+ (*Base::A)[iElemIdxSLE][neighborIdxSLE] += -mu * dEvaluation2;
+ (*Base::A)[neighborIdxSLE][iElemIdxSLE] += -mu * dEvaluation2;
for (int i = 0; i < dim - 1; i++)
{
@@ -522,21 +536,21 @@ private:
const Scalar interfaceUpdate1 = mu * centerI;
const Scalar interfaceUpdate2 = KparDotTau_i * normal;
const Scalar interfaceUpdate3 =
- -interfaceUpdate1 + 0.5 * interfaceUpdate2;
+ dEvaluation2 * (-interfaceUpdate1 + 0.5 * interfaceUpdate2);
const Scalar interfaceUpdate4 =
- -interfaceUpdate1 - 0.5 * interfaceUpdate2;
+ dEvaluation2 * (-interfaceUpdate1 - 0.5 * interfaceUpdate2);
(*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] +=
interfaceUpdate3;
(*Base::A)[neighborIdxSLE][iElemIdxSLE + i + 1] +=
interfaceUpdate3;
(*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] +=
- -interfaceUpdate4;
+ interfaceUpdate4;
(*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] +=
- -interfaceUpdate4;
+ interfaceUpdate4;
(*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + i + 1] +=
- -centerI * interfaceUpdate1;
+ -centerI * dEvaluation2 * interfaceUpdate1;
for (int j = 0; j < i; j++)
{ //NOTE: only relevant for n=3, requires quadrature rule
@@ -547,9 +561,9 @@ private:
const Scalar interfaceUpdate6 = interfaceUpdate2 *
(center * iFrame[j]) - (KparDotTau_j * normal) * centerI;
const Scalar interfaceUpdate7 =
- -interfaceUpdate5 - 0.5 * interfaceUpdate6;
+ dEvaluation2 * (-interfaceUpdate5 - 0.5 * interfaceUpdate6);
const Scalar interfaceUpdate8 =
- -interfaceUpdate5 + 0.5 * interfaceUpdate6;
+ dEvaluation2 * (-interfaceUpdate5 + 0.5 * interfaceUpdate6);
(*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + j + 1] +=
interfaceUpdate7;
@@ -574,10 +588,10 @@ private:
// int_intersct d * g * (mu * phi_iElem,i +
// (Kparallel * grad(phi_i,iElem)) * normal) dr
//for i = 0,...,dim-1
- const Scalar dEvaluation = Base::problem_.aperture(center);
- const Scalar dgGamma = dEvaluation * Base::problem_.boundary(center);
+ const Scalar d2gGamma =
+ dEvaluation2 * Base::problem_.interfaceBoundary(center);
- Base::b[iElemIdxSLE] += mu * dgGamma * dEvaluation;
+ Base::b[iElemIdxSLE] += mu * d2gGamma;
for (int i = 0; i < dim - 1; i++)
{
@@ -586,26 +600,27 @@ private:
const Scalar interfaceUpdate1 = mu * centerI;
const Scalar interfaceUpdate2 = KparDotTau_i * normal;
const Scalar interfaceUpdate3 =
- interfaceUpdate1 - interfaceUpdate2;
+ dEvaluation2 * (interfaceUpdate1 - interfaceUpdate2);
Base::b[iElemIdxSLE + i + 1] +=
- dgGamma * dEvaluation * (mu * centerI - interfaceUpdate2);
+ d2gGamma * (mu * centerI - interfaceUpdate2);
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
interfaceUpdate3;
(*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
interfaceUpdate3;
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
- centerI * (interfaceUpdate1 - 2.0 * interfaceUpdate2);
+ centerI * dEvaluation2 *
+ (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 = dEvaluation2 *
+ ((center * iFrame[j]) * interfaceUpdate1
- (interfaceUpdate2 * (center * iFrame[j])
- + (KparDotTau_j * normal) * centerI);
+ + (KparDotTau_j * normal) * centerI));
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
interfaceUpdate4;