diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh index ccaef4d5945090dba159bf9a51f19f741a2d4813..4221388e2146812bea7dc63d3f2f67d3dccfb942 100644 --- a/dune/mmdg/mmdg.hh +++ b/dune/mmdg/mmdg.hh @@ -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;