diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh index 8b05f05bd75ef8a3c12a45eef422be15fc3dada9..b6d5c222ba74ef10883cadd072bf576aa6ff228e 100644 --- a/dune/mmdg/mmdg.hh +++ b/dune/mmdg/mmdg.hh @@ -10,15 +10,19 @@ class MMDG : public DG<GridView, Mapper, Problem> public: static constexpr int dim = GridView::dimension; static constexpr auto iSimplex = Dune::GeometryTypes::simplex(dim-1); + static constexpr auto iFacet = Dune::GeometryTypes::simplex(dim-2); using Base = DG<GridView, Mapper, Problem>; using Scalar = typename Base::Scalar; using Matrix = typename Base::Matrix; using Vector = typename Base::Vector; - using Coordinate = Dune::FieldVector<Scalar, dim>; + using Coordinate = typename Base::Coordinate; + using CoordinateMatrix = typename Base::CoordinateMatrix; using InterfaceBasis = std::array<Coordinate, dim-1>; using QRuleInterface = typename Base::QRuleEdge; using QRulesInterface = typename Base::QRulesEdge; + using QRuleIFacet = Dune::QuadratureRule<Scalar, dim-2>; + using QRulesIFacet = Dune::QuadratureRules<Scalar, dim-2>; using VTKFunction = Dune::VTKFunction<InterfaceGridView>; using P1Function = Dune::NonconformingP1VTKFunction<InterfaceGridView, Dune::DynamicVector<Scalar>>; @@ -114,6 +118,10 @@ private: Base::problem_.quadratureOrder_Kperp()); const QRuleInterface& rule_KperpPlus1 = QRulesInterface::rule(iSimplex, Base::problem_.quadratureOrder_Kperp() + 1); + const QRuleIFacet& rule_Kpar_iFacet = + QRulesIFacet::rule(iFacet, Base::problem_.quadratureOrder_Kparallel()); + const QRuleIFacet& rule_interfaceBoundary = QRulesIFacet::rule(iFacet, + Base::problem_.quadratureOrderInterfaceBoundary()); //we use the bulk basis // phi_elem,0 (x) = indicator(elem); @@ -463,193 +471,297 @@ private: Base::applyQuadratureRule(iGeo, rule_Kpar, lambda_Kpar); //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; + const auto& intersctGeo = intersct.geometry(); //determine penalty parameter const Scalar mu = getPenaltyParameter(mu0, intersct, iGeoCenter); - const Scalar dEvaluation = Base::problem_.aperture(center); - const Scalar dEvaluation2 = dEvaluation * dEvaluation; - - //evaluate the interface term - // int_intersct mu^Gamma * d^2 * jump(phi_iElem,0) - // * jump(phi_iElem,0) dr - (*Base::A)[iElemIdxSLE][iElemIdxSLE] += mu * dEvaluation2; - if (intersct.neighbor()) //interior interface edge { - //evaluate the interface terms //TODO terms with grad(d) - // int_intersct mu^Gamma * jump(phi_iElem,i) * jump(phi_iElem,j) dr + //index of the neighboring element + const int neighborIdx = iMapper_.index(intersct.outside()); + + //lambda to evaluate the interface terms + // int_intersct mu * d^2 * jump(phi_iElem1,i) * jump(phi_iElem2,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 = - 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] += - dEvaluation2 * centerI * (interfaceUpdate1 - interfaceUpdate2); + // -int_intersct d * avg(Kparallel * grad(d * phi_iElem1,i)) + // * jump(phi_iElem2,j) dr + //for iElem1, iElem2 in {iElem, neighbor} and for i,j = 0,...,dim-1 + //using a quadrature rule //TODO: terms with grad(d) + const auto lambda_Kpar_iFacet = + [&](const Coordinate& qpGlobal, const Scalar weight) + { + const Scalar dEvaluation = Base::problem_.aperture(qpGlobal); + const Scalar dEvaluation2 = dEvaluation * dEvaluation; + const CoordinateMatrix KEvaluation = + Base::problem_.Kparallel(qpGlobal); - 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 = dEvaluation2 * - ((center * iFrame[j]) * interfaceUpdate1 - - 0.5 * (interfaceUpdate2 * (center * iFrame[j]) - + (KparDotTau_j * normal) * centerI)); + Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i] + Coordinate KparDotTau_j; //storage for K_parallel * iFrame[j] - (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] += - interfaceUpdate4; - (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] += - interfaceUpdate4; - } - } + //quadrature for + // int_intersct mu * d^2 * jump(phi_iElem,0) + // * jump(phi_iElem,0) dr + (*Base::A)[iElemIdxSLE][iElemIdxSLE] += + weight * mu * dEvaluation2; - //index of the neighboring element - const int neighborIdx = iMapper_.index(intersct.outside()); + for (int i = 0; i < dim - 1; i++) + { + KEvaluation.mv(iFrame[i], KparDotTau_i); - if (neighborIdx > iElemIdx) - { //make sure that each interface edge is considered only once - continue; - } + const Scalar qpI = qpGlobal * iFrame[i]; + const Scalar interfaceUpdate1 = mu * qpI; + const Scalar interfaceUpdate2 = KparDotTau_i * normal; + const Scalar interfaceUpdate3 = weight * dEvaluation2 + * (interfaceUpdate1 - 0.5 * interfaceUpdate2); - const int neighborIdxSLE = - (dim + 1) * Base::gridView_.size(0) + dim * neighborIdx; + //quadrature for + // int_intersct mu * jump(phi_iElem,0) * jump(phi_iElem,i) dr + //and + // -int_intersct avg(Kparallel * grad(d * phi_iElem,i)) + // * jump(phi_iElem,0) dr //TODO: terms with grad(d) + (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] += + interfaceUpdate3; + (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] += + interfaceUpdate3; - //get basis of the interface coordinate system of the neighboring - //interface element, for uncurved interfaces the sign of normals - //and hence the frames of iElem and its neighbor can differ - const InterfaceBasis neighborFrame = - interfaceFrame(Base::gridView_.grid().asIntersection( - intersct.outside()).centerUnitOuterNormal()); + //quadrature for + // int_intersct mu * jump(phi_iElem,i)^2 dr + //and + // -int_intersct avg(Kparallel * grad(d * phi_iElem,i)) + // * jump(phi_iElem,i) dr //TODO: terms with grad(d) + (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] += + weight * dEvaluation2 * qpI * + (interfaceUpdate1 - interfaceUpdate2); + + for (int j = 0; j < i; j++) + { + KEvaluation.mv(iFrame[j], KparDotTau_j); + + const Scalar interfaceUpdate4 = weight * dEvaluation2 * + ((qpGlobal * iFrame[j]) * interfaceUpdate1 + - 0.5 * (interfaceUpdate2 * (qpGlobal * iFrame[j]) + + (KparDotTau_j * normal) * qpI)); + + //quadrature for + // int_intersct mu * jump(phi_iElem,i) * jump(phi_iElem,j) dr + //and + // -int_intersct avg(Kparallel * grad(d * phi_iElem,i)) + // * jump(phi_iElem,j) dr //TODO: terms with grad(d) + (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] += + interfaceUpdate4; + (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] += + interfaceUpdate4; + } + } - //storage for K_par * neighborFrame[i] - Coordinate KparDotNeighborTau_i; + if (neighborIdx > iElemIdx) + { //make sure that each interface edge is considered only once + return; + } - //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)) - // * 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)[iElemIdxSLE][neighborIdxSLE] += -mu * dEvaluation2; - (*Base::A)[neighborIdxSLE][iElemIdxSLE] += -mu * dEvaluation2; + const int neighborIdxSLE = + (dim + 1) * Base::gridView_.size(0) + dim * neighborIdx; + const InterfaceBasis neighborFrame = + interfaceFrame(Base::gridView_.grid().asIntersection( + intersct.outside()).centerUnitOuterNormal()); - for (int i = 0; i < dim - 1; i++) - { - Base::problem_.Kparallel(center).mv(iFrame[i], KparDotTau_i); - Base::problem_.Kparallel(center).mv( - neighborFrame[i], KparDotNeighborTau_i); - const Scalar centerI = center * iFrame[i]; - const Scalar neighborCenterI = center * neighborFrame[i]; - const Scalar interfaceUpdate1a = mu * centerI; - const Scalar interfaceUpdate1b = mu * neighborCenterI; - const Scalar interfaceUpdate2a = KparDotTau_i * normal; - const Scalar interfaceUpdate2b = KparDotNeighborTau_i * normal; - const Scalar interfaceUpdate3a = - dEvaluation2 * (-interfaceUpdate1a + 0.5 * interfaceUpdate2a); - const Scalar interfaceUpdate3b = - dEvaluation2 * (-interfaceUpdate1b - 0.5 * interfaceUpdate2b); - - (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] += - interfaceUpdate3a; - (*Base::A)[neighborIdxSLE][iElemIdxSLE + i + 1] += - interfaceUpdate3a; - (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] += - interfaceUpdate3b; - (*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] += - interfaceUpdate3b; - - for (int j = 0; j < dim - 1; j++) - { - Coordinate KparDotNeighborTau_j; + //storage for K_par * neighborFrame[i] + Coordinate KparDotNeighborTau_i; - Base::problem_.Kparallel(center).mv( - neighborFrame[j], KparDotNeighborTau_j); + //quadrature for + // int_intersct mu * d^2 * jump(phi_iElem,0) + // * jump(phi_neighbor,0) dr + (*Base::A)[iElemIdxSLE][neighborIdxSLE] -= + weight * mu * dEvaluation2; + (*Base::A)[neighborIdxSLE][iElemIdxSLE] -= + weight * mu * dEvaluation2; + + for (int i = 0; i < dim - 1; i++) + { + KEvaluation.mv(iFrame[i], KparDotTau_i); + KEvaluation.mv(neighborFrame[i], KparDotNeighborTau_i); + + const Scalar qpI = qpGlobal * iFrame[i]; + const Scalar neighborQpI = qpGlobal * neighborFrame[i]; + const Scalar interfaceUpdate1a = mu * qpI; + const Scalar interfaceUpdate1b = mu * neighborQpI; + const Scalar interfaceUpdate2a = KparDotTau_i * normal; + const Scalar interfaceUpdate2b = KparDotNeighborTau_i * normal; + const Scalar interfaceUpdate3a = weight * dEvaluation2 * + (-interfaceUpdate1a + 0.5 * interfaceUpdate2a); + const Scalar interfaceUpdate3b = weight * dEvaluation2 * + (-interfaceUpdate1b - 0.5 * interfaceUpdate2b); - const Scalar interfaceUpdate5 = - -dEvaluation2 * ((center * neighborFrame[j]) * interfaceUpdate1a - + 0.5 * (-interfaceUpdate2a * (center * neighborFrame[j]) - + (KparDotNeighborTau_j * normal) * centerI)); + //quadrature for + // int_intersct mu * d^2 * jump(phi_iElem,i) + // * jump(phi_neighbor,0) dr + //and + // -int_intersct d * avg(Kparallel * grad(d * phi_iElem,i)) + // * jump(phi_neighbor,0) dr //TODO: terms with grad(d) + (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] += + interfaceUpdate3a; + (*Base::A)[neighborIdxSLE][iElemIdxSLE + i + 1] += + interfaceUpdate3a; - (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + j + 1] += - interfaceUpdate5; - (*Base::A)[neighborIdxSLE + j + 1][iElemIdxSLE + i + 1] += - interfaceUpdate5; - } - } + //quadrature for + // int_intersct mu * d^2 * jump(phi_iElem,0) + // * jump(phi_neighbor,i) dr + //and. + // -int_intersct d * avg(Kparallel * grad(d * phi_neighbor,i)) + // * jump(phi_iElem,0) dr //TODO: terms with grad(d) + (*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] += + interfaceUpdate3b; + (*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] += + interfaceUpdate3b; + + for (int j = 0; j < dim - 1; j++) + { + Coordinate KparDotNeighborTau_j; + + KEvaluation.mv(neighborFrame[j], KparDotNeighborTau_j); + + const Scalar interfaceUpdate5 = -weight * dEvaluation2 + * ((qpGlobal * neighborFrame[j]) * interfaceUpdate1a + + 0.5 * (-interfaceUpdate2a * (qpGlobal * neighborFrame[j]) + + (KparDotNeighborTau_j * normal) * qpI)); + + //quadrature for + // int_intersct mu * d^2 * jump(phi_iElem,i) + // * jump(phi_neighbor,j) dr + //and + // -int_intersct d * avg(Kparallel * grad(d * phi_iElem,i)) + // * jump(phi_neighbor,j) dr + //resp. + // -int_intersct d * avg(Kparallel * grad(d * phi_neighbor,i)) + // * jump(phi_iElem,j) dr + //TODO: terms with grad(d) + (*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE + j + 1] += + interfaceUpdate5; + (*Base::A)[neighborIdxSLE + j + 1][iElemIdxSLE + i + 1] += + interfaceUpdate5; + } + } + }; + + Base::applyQuadratureRule(intersctGeo, rule_Kpar_iFacet, + lambda_Kpar_iFacet); } else //boundary interface edge { - //evaluate the interface terms - // int_intersct mu^Gamma * phi_iElem,i * phi_iElem,j dr + //lambda to evaluate the interface terms on the boundary + // int_intersct mu * d^2 * 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 d2gGamma = - dEvaluation2 * Base::problem_.interfaceBoundary(center); - - Base::b[iElemIdxSLE] += mu * d2gGamma; + // -int_intersct d * (Kparallel * grad(d * phi_iElem,i)) + // * phi_iElem,j * normal dr + //for i,j = 0,...,dim-1 using a quadrature rule + //TODO: terms with grad(d) + const auto lambda_Kpar_iFacet_boundary = + [&](const Coordinate& qpGlobal, const Scalar weight) + { + const Scalar dEvaluation = Base::problem_.aperture(qpGlobal); + const Scalar dEvaluation2 = dEvaluation * dEvaluation; + const CoordinateMatrix KEvaluation = + Base::problem_.Kparallel(qpGlobal); - 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 = - dEvaluation2 * (interfaceUpdate1 - interfaceUpdate2); + Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i] + Coordinate KparDotTau_j; //storage for K_parallel * iFrame[j] - Base::b[iElemIdxSLE + i + 1] += - d2gGamma * (interfaceUpdate1 - interfaceUpdate2); + //quadrature for + // int_intersct mu * d^2 * jump(phi_iElem,0) + // * jump(phi_iElem,0) dr + (*Base::A)[iElemIdxSLE][iElemIdxSLE] += + weight * mu * dEvaluation2; - (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] += - interfaceUpdate3; - (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] += - interfaceUpdate3; - (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] += - 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 = dEvaluation2 * - ((center * iFrame[j]) * interfaceUpdate1 - - (interfaceUpdate2 * (center * iFrame[j]) - + (KparDotTau_j * normal) * centerI)); + for (int i = 0; i < dim - 1; i++) + { + KEvaluation.mv(iFrame[i], KparDotTau_i); - (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] += - interfaceUpdate4; - (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] += - interfaceUpdate4; - } - } + const Scalar qpI = qpGlobal * iFrame[i]; + const Scalar interfaceUpdate1 = mu * qpI; + const Scalar interfaceUpdate2 = KparDotTau_i * normal; + const Scalar interfaceUpdate3 = + weight * dEvaluation2 * (interfaceUpdate1 - interfaceUpdate2); + + //quadrature for + // int_intersct mu * d^2 * phi_iElem,0 * phi_iElem,i dr + //and + // -int_intersct d * (Kparallel * grad(d * phi_iElem,i)) + // * phi_iElem,0 * normal dr + (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] += + interfaceUpdate3; + (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] += + interfaceUpdate3; + + //quadrature for + // int_intersct mu * d^2 * (phi_iElem,i)^2 dr + //and + // -int_intersct d * (Kparallel * grad(d * phi_iElem,i)) + // * phi_iElem,i * normal dr + (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] += + weight * qpI * dEvaluation2 * + (interfaceUpdate1 - 2.0 * interfaceUpdate2); + + for (int j = 0; j < i; i++) + { + KEvaluation.mv(iFrame[j], KparDotTau_j); + + const Scalar interfaceUpdate4 = weight * dEvaluation2 * + ((qpGlobal * iFrame[j]) * interfaceUpdate1 + - (interfaceUpdate2 * (qpGlobal * iFrame[j]) + + (KparDotTau_j * normal) * qpI)); + + //quadrature for + // int_intersct mu * d^2 * phi_iElem,i * phi_iElem,j dr + //and + // -int_intersct d * (Kparallel * grad(d * phi_iElem,i)) + // * phi_iElem,j * normal dr + (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] += + interfaceUpdate4; + (*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] += + interfaceUpdate4; + } + } + }; + + //lambda to evaluate the interface boundary terms of the load vector b + // int_intersct mu * d^2 * g^Gamma * phi_iElem,i dr + //and + // int_intersct d * g^Gamma * K_parallel * grad(d * phi_iElem,i) dr + //for i = 0,...,dim-1 using a quadrature rule + //TODO: terms with grad(d) + const auto lambda_interfaceBoundary = + [&](const Coordinate& qpGlobal, const Scalar weight) + { + const Scalar dEvaluation = Base::problem_.aperture(qpGlobal); + const Scalar dEvaluation2 = dEvaluation * dEvaluation; + const CoordinateMatrix KEvaluation = + Base::problem_.Kparallel(qpGlobal); + const Scalar d2gGamma = weight * dEvaluation2 + * Base::problem_.interfaceBoundary(qpGlobal); + + //storage for K_parallel * iFrame[i] + Coordinate KparDotTau_i; + + Base::b[iElemIdxSLE] += mu * d2gGamma; + + for (int i = 0; i < dim - 1; i++) + { + KEvaluation.mv(iFrame[i], KparDotTau_i); + + Base::b[iElemIdxSLE + i + 1] += d2gGamma * + (mu * (qpGlobal * iFrame[i]) - KparDotTau_i * normal); + } + }; + + //apply quadrature rules on the boundary interface facet intersct + Base::applyQuadratureRule(intersctGeo, rule_Kpar_iFacet, + lambda_Kpar_iFacet_boundary); + Base::applyQuadratureRule(intersctGeo, rule_interfaceBoundary, + lambda_interfaceBoundary); } } }