diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh index 03077fd37829e3aa99894dafcd820aad78787c9a..102d0e254d6286ee7eaa6ed4ee3328a3a9d38823 100644 --- a/dune/mmdg/mmdg.hh +++ b/dune/mmdg/mmdg.hh @@ -384,11 +384,44 @@ private: } //i = dim - 1 - const Scalar couplingUpdate4 = qpGlobal[dim] * couplingUpdate1; - (*Base::A)[elemInIdxSLE + dim + 1][iElemIdxSLE] -= couplingUpdate4; - (*Base::A)[iElemIdxSLE][elemInIdxSLE + dim + 1] -= couplingUpdate4; - (*Base::A)[elemOutIdxSLE + dim + 1][iElemIdxSLE] -= couplingUpdate4; - (*Base::A)[iElemIdxSLE][elemOutIdxSLE + dim + 1] -= couplingUpdate4; + 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 + { //TODO TODO TODO + 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] -= + 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; + } + } } // === interface entries ===