diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 044da9845ae966c3801c3f9f6d71545ccd1ae126..5892fa85f95ce721a010292df2ff7a1f2f01be1f 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -11,7 +11,7 @@ public:
static constexpr int dim = GridView::dimension;
using Base = DG<Grid, GridView, Mapper, Problem>;
- using Scalar = typename Base::Scalar;
+ using Scalar = typename Base::Scalar; //NOTE using Base::Scalar
using Matrix = typename Base::Matrix;
using Vector = typename Base::Vector;
using Coordinate = Dune::FieldVector<Scalar, dim>;
@@ -105,60 +105,6 @@ 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)
- {
- //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));
-
- //quadrature for int_elem q*phi_elem,0 dV
- 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] += quadratureFactor *
- (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
- }
- }
-
- //TODO: fill stiffness matrix A with interface entries
-
- //evaluate
- // int_iElem d * (K_par*grad_par(phi_iElem,i)) * grad_par(phi_iElem,j) ds
- // = int_iElem d * (K_par * tau_i) * tau_j ds
- //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);
-
- for (int i = 0; i < dim-1 ; i++)
- {
- Coordinate K_parDotTau_i;
- Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], K_parDotTau_i);
-
- for (int j = 0; j <= i; j++) //TODO: symmetrize
- {
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
- quadratureFactor * Base::problem_.aperture(qpGlobal)
- * ( K_parDotTau_i * iFrame[j] );
- }
- }
- }
-
- //TODO: check interface bilinear form
-
// === coupling entries ===
//evaluate the coupling term
@@ -373,9 +319,239 @@ private:
Base::A->entry(elemOutIdxSLE + dim + 1, iElemIdx) -= couplingUpdate4;
Base::A->entry(iElemIdx, elemOutIdxSLE + dim + 1) -= couplingUpdate4;
}
+
+ // === 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));
+
+ //quadrature for int_elem q*phi_elem,0 dV
+ 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] += quadratureFactor *
+ (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
+ }
+ }
+
+ //evaluate
+ // int_iElem (Kparallel * grad(phi_iElem,i)) * grad(phi_iElem,j) ds
+ // = int_iElem (Kpar * tau_i) * tau_j ds
+ //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);
+
+ 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->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ quadratureFactor * ( KparDotTau_i * iFrame[j] );
+
+ if (i != j)
+ {
+ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ quadratureFactor * ( KparDotTau_i * iFrame[j] );
+ }
+ }
+ }
+ }
+
+ //iterate over all intersection with the boundary of iElem
+ //NOTE: only works for 2d!, otherwise quadrature required
+ for (const auto& intersct : intersections(iGridView_, iElem))
+ {
+ const auto& center = intersct.geometry().center();
+ const auto& normal = intersct.centerUnitOuterNormal();
+ Coordinate KparDotTau_i;
+
+ //evaluate the interface term
+ // int_intersct mu^Gamma * jump(phi_iElem,0) * jump(phi_iElem,0) dr
+ Base::A->entry(iElemIdxSLE, iElemIdxSLE) += mu;
+
+ if (intersct.neighbor()) //inner interface edge
+ {
+ //evaluate the interface terms
+ // int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_iElem,j) dr
+ //and
+ // -int_intersct avg(Kparallel * grad(phi_iElem,i))
+ // * jump(phi_iElem,j) dr
+ //for i,j = 0,...,dim-1 and (i != 0 or j != 0)
+ for (int i = 0; i < dim - 1; i++)
+ {
+ Base::problem_.Kparallel(center).mv(iFrame[i], KparDotTau_i);
+ const Scalar centerI = center * iFrame[i];
+ const Scalar interfaceUpdate1 = mu * centerI;
+ const Scalar interfaceUpdate2 = KparDotTau_i * normal;
+ const Scalar interfaceUpdate3 =
+ interfaceUpdate1 - 0.5 * 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 - 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
+ - 0.5 * (interfaceUpdate2 * (center * iFrame[j])
+ + (KparDotTau_j * normal) * centerI);
+
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ interfaceUpdate4;
+ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ interfaceUpdate4;
+ }
+ }
+
+ //index of the neighboring element
+ const int neighborIdx = iMapper_.index(intersct.outside());
+
+ if (neighborIdx > iElemIdx)
+ { //make sure that each interface edge is considered only once
+ continue;
+ }
+
+ const int neighborIdxSLE =
+ (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
+ //and
+ // -int_intersct avg(Kparallel * grad(phi_iElem,i))
+ // * jump(phi_neighbor,j) dr
+ //resp.
+ // -int_intersct avg(Kparallel * grad(phi_neighbor,i))
+ // * jump(phi_iElem,j) dr
+ //for i,j = 0,...,dim-1
+ Base::A->entry(iElemIdxSLE, neighborIdxSLE) += -mu;
+ Base::A->entry(neighborIdxSLE, iElemIdxSLE) += -mu;
+
+ for (int i = 0; i < dim - 1; i++)
+ {
+ Base::problem_.Kparallel(center).mv(iFrame[i], KparDotTau_i);
+ const Scalar centerI = center * iFrame[i];
+ const Scalar interfaceUpdate1 = mu * centerI;
+ const Scalar interfaceUpdate2 = KparDotTau_i * normal;
+ const Scalar interfaceUpdate3 =
+ -interfaceUpdate1 + 0.5 * interfaceUpdate2;
+ const Scalar interfaceUpdate4 =
+ -interfaceUpdate1 - 0.5 * interfaceUpdate2;
+
+ Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE) +=
+ interfaceUpdate3;
+ Base::A->entry(neighborIdxSLE, iElemIdxSLE + i + 1) +=
+ interfaceUpdate3;
+ Base::A->entry(neighborIdxSLE + i + 1, iElemIdxSLE) +=
+ -interfaceUpdate4;
+ Base::A->entry(iElemIdxSLE, neighborIdxSLE + i + 1) +=
+ -interfaceUpdate4;
+
+ Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + i + 1) +=
+ -centerI * interfaceUpdate1;
+
+ for (int j = 0; j < i; j++)
+ { //NOTE: only relevant for n=3, requires quadrature rule
+ Coordinate KparDotTau_j;
+ Base::problem_.Kparallel(center).mv(iFrame[j], KparDotTau_j);
+ const Scalar interfaceUpdate5 =
+ (center * iFrame[j]) * interfaceUpdate1;
+ const Scalar interfaceUpdate6 = interfaceUpdate2 *
+ (center * iFrame[j]) - (KparDotTau_j * normal) * centerI;
+ const Scalar interfaceUpdate7 =
+ -interfaceUpdate5 - 0.5 * interfaceUpdate6;
+ const Scalar interfaceUpdate8 =
+ -interfaceUpdate5 + 0.5 * interfaceUpdate6;
+
+ Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
+ interfaceUpdate7;
+ Base::A->entry(neighborIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ interfaceUpdate7;
+ Base::A->entry(iElemIdxSLE + j + 1, neighborIdxSLE + i + 1) +=
+ interfaceUpdate8;
+ Base::A->entry(neighborIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ interfaceUpdate8;
+ }
+ }
+ }
+ else //boundary interface edge
+ {
+ //evaluate the interface terms
+ // int_intersct mu^Gamma * phi_iElem,i * phi_iElem,j dr
+ //and
+ // -int_intersct phi_iElem,j * (Kparallel * grad(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
+ const Scalar dgGamma = Base::problem_.aperture(center)
+ * Base::problem_.boundary(center);
+ Base::b[iElemIdxSLE] += mu * dgGamma;
+
+ for (int i = 0; i < dim - 1; i++)
+ {
+ Base::problem_.Kparallel(center).mv(iFrame[i], KparDotTau_i);
+ const Scalar centerI = center * iFrame[i];
+ const Scalar interfaceUpdate1 = mu * centerI;
+ const Scalar interfaceUpdate2 = KparDotTau_i * normal;
+ const Scalar interfaceUpdate3 =
+ interfaceUpdate1 - interfaceUpdate2;
+
+ Base::b[iElemIdxSLE] += dgGamma * (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);
+
+ 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
+ - (interfaceUpdate2 * (center * iFrame[j])
+ + (KparDotTau_j * normal) * centerI);
+
+ Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
+ interfaceUpdate4;
+ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
+ interfaceUpdate4;
+ }
+ }
+ }
+ }
}
}
+//NOTE: MMesh ensures unique intersection, i.e. unique normal, add assert index(inside) < index(outside)
//returns a basis of the interface coordinate system
//2d: tau = [ normal[1],
// -normal[0] ]
@@ -399,39 +575,16 @@ private:
tau[i] += normal[j];
tau[j] = -normal[i];
}
- tau /= tau.two_norm();
+
+ if constexpr (dim != 2)
+ {
+ tau /= tau.two_norm();
+ }
iBasis[i] = tau;
}
return iBasis;
-
-/* //NOTE: which version is better?
- if constexpr (dim == 2)
- {
- Coordinate tau;
- tau[0] = normal[1];
- tau[1] = -normal[0];
-
- interfaceFrame = {tau};
- }
- else if constexpr (dim == 3)
- {
- Coordinate tau1;
- tau1[0] = normal[1] + normal[2];
- tau1[1] = -normal[0];
- tau1[2] = -normal[0];
- Coordinate tau2;
- tau2[0] = -normal[1];
- tau2[1] = normal[0] + normal[2];
- tau2[2] = -normal[1];
- interfaceFrame = {tau1, tau2};
- }
- else
- {
- DUNE_THROW(Dune::Exception,
- "Invalid dimension: dim = " + std::to_string(dim));
- }*/
}
};
diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh
index ad26400cc62d663d9948d01156fa6f963c736c30..6e23c898c3e04ae11b07ccb1e64f5ef8ce770a3a 100644
--- a/dune/mmdg/problems/coupleddgproblem.hh
+++ b/dune/mmdg/problems/coupleddgproblem.hh
@@ -61,7 +61,7 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
}
//returns the recommended quadrature order to compute an integral
- //over d(x) * Kparallel(x)
+ //over Kparallel(x)
virtual int quadratureOrder_Kparallel () const
{
return 1;