Select Git revision
mmdg.hh 38.86 KiB
#ifndef DUNE_MMDG_HH
#define DUNE_MMDG_HH
#include <dune/mmdg/dg.hh>
template<class GridView, class Mapper, class InterfaceGridView,
class InterfaceMapper, class Problem>
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 = 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>>;
//constructor
MMDG (const GridView& gridView, const Mapper& mapper,
const InterfaceGridView& iGridView, const InterfaceMapper& iMapper,
const Problem& problem) :
Base(gridView, mapper, problem,
(1 + dim) * gridView.size(0) + dim * iGridView.size(0)),
iGridView_(iGridView), iMapper_(iMapper) {}
const void operator() (const Scalar mu0, const Scalar xi)
{
assembleSLE(mu0, xi);
Base::A->compress();
Dune::InverseOperatorResult result;
Dune::UMFPack<Matrix> solver(*Base::A);
solver.apply(Base::d, Base::b, result);
//write solution to a vtk file
writeVTKoutput();
}
//compute L2 error using quadrature rules
// norm = sqrt( ||.||^2_L2(bulk) + ||.||^2_L2(interface) )
Scalar computeL2error(const int quadratureOrder = 5)
{
const Scalar bulkError = Base::computeL2error(quadratureOrder);
Scalar interfaceError(0.0);
for (const auto& iElem : elements(iGridView_))
{
const auto& iGeo = iElem.geometry();
const int iElemIdxSLE =
(dim + 1) * Base::gridView_.size(0) + dim * iMapper_.index(iElem);
const QRuleInterface& qrule =
QRulesInterface::rule(iSimplex, quadratureOrder);
const InterfaceBasis iFrame = interfaceFrame(
Base::gridView_.grid().asIntersection(iElem).centerUnitOuterNormal());
//determine L2-error on the interface element iElem using a
//quadrature rule
const auto qlambda = [&](const Coordinate& qpGlobal, const Scalar weight)
{
Scalar numericalSolutionAtQP = Base::d[iElemIdxSLE];
for (int i = 0; i < dim - 1; i++)
{
numericalSolutionAtQP +=
Base::d[iElemIdxSLE + i + 1] * (qpGlobal * iFrame[i]);
}
const Scalar errorEvaluation = numericalSolutionAtQP
- Base::problem_.exactInterfaceSolution(qpGlobal);
interfaceError += weight * errorEvaluation * errorEvaluation;
};
Base::applyQuadratureRule(iGeo, qrule, qlambda);
}
//print bulk and interface error
std::cout << "ERROR bulk: " << bulkError << "\t ERROR interface: " <<
sqrt(interfaceError) << "\n";
return sqrt(bulkError * bulkError + interfaceError);
}
const InterfaceGridView& iGridView_;
const InterfaceMapper& iMapper_;
private:
//assemble stiffness matrix A and load vector b
void assembleSLE (const Scalar mu0, const Scalar xi)
{
//modell parameter
const auto beta = [this, &xi](const Coordinate& pos) -> Scalar
{
return 4.0 * Base::problem_.Kperp(pos) / (2.0 * xi - 1.0);
};
//quadrature rules
const QRuleInterface& rule_qInterface = QRulesInterface::rule(iSimplex,
Base::problem_.quadratureOrder_qInterface());
const QRuleInterface& rule_Kpar = QRulesInterface::rule(iSimplex,
Base::problem_.quadratureOrder_Kparallel());
const QRuleInterface& rule_Kperp = QRulesInterface::rule(iSimplex,
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);
// phi_elem,i (x) = x[i]*indicator(elem);
//for elem in elements(gridView) and i = 1,...,dim
//we use the interface basis
// phi_iElem,0 (x) = indicator(elem);
// phi_iElem,i (x) = (x * tau_i)*indicator(elem);
//for iElem in elements(iGridView) and i = 1,...,dim-1
//where {tau_i | i=1,...,dim-1} is the interface coordinate system given
//by the method interfaceFrame(normal)
//in total we use the basis given by the union of
// { (phi_elem,i , 0) | elem in elements(gridView), i = 0,...,dim }
//and
// { (0 , phi_iElem,i) | iElem in elements(iGridView), i = 0,...,dim-1 }
//thus the stiffness matrix A and the load vector b exhibit the following
//block structure
// A = [ bulk, coupling b = [ bulk,
// coupling, interface], interface ]
//call base class method for the assignment to fill stiffness matrix A
//and load vector b with bulk entries
Base::assembleSLE(mu0);
for (const auto& iElem : elements(iGridView_))
{
const int iElemIdx = iMapper_.index(iElem);
const auto& iGeo = iElem.geometry();
const auto& iGeoCenter = iGeo.center();
const auto& intersct = Base::gridView_.grid().asIntersection(iElem);
//get basis of the interface coordinate system for evaluation of
//basis function phi_iElem,i with i=1,...,dim-1
const InterfaceBasis iFrame =
interfaceFrame(intersct.centerUnitOuterNormal());
//in the total system of linear equations (SLE) Ad = b,
//the index iElemIdxSLE refers to the basis function (0, phi_iElem,0)
//and the indices iElemIdxSLE + i + 1, i = 1,...,dim-1, refer to the
//basis functions (0, phi_iElem,i)
const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;
//get neighboring bulk grid elements of iELem
const auto elemIn = intersct.inside();
const auto elemOut = intersct.outside();
const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn);
const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut);
// === coupling entries ===
//lambda to evaluate the coupling terms
// int_iElem beta * phi_iElem,0 * phi_iElem,j ds
//and
// int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds
//and
// int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds
//and
// -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
//and
// -int_iElem beta * phi_iElem,j * avg(phi_elem,0) ds
//for i = 0,...,dim and j = 0,...,dim-1
//and for elem in {elemIn, elemOut} using a quadrature rule
//TODO symmetrize more efficiently
const auto lambda_Kperp =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
const Scalar betaEvaluation = beta(qpGlobal);
const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
const Scalar bulkUpdate1 =
weight * (KEvaluation + 0.25 * betaEvaluation);
const Scalar bulkUpdate2 =
weight * (-KEvaluation + 0.25 * betaEvaluation);
const Scalar couplingUpdate1 = 0.5 * weight * betaEvaluation;
//quadrature for
// int_iElem beta * phi_iElem,0 * phi_iElem,0 ds
Base::A->entry(iElemIdxSLE, iElemIdxSLE) += weight * betaEvaluation;
//quadrature for
// int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds
//and
// int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,0) ds
//for elem1, elem2 in {elemIn, elemOut} Base::A->entry(elemInIdxSLE, elemInIdxSLE) += bulkUpdate1;
Base::A->entry(elemOutIdxSLE, elemOutIdxSLE) += bulkUpdate1;
Base::A->entry(elemInIdxSLE, elemOutIdxSLE) += bulkUpdate2;
Base::A->entry(elemOutIdxSLE, elemInIdxSLE) += bulkUpdate2;
//quadrature for
// -int_iElem beta * phi_iElem,0 * avg(phi_elem,j) ds
//for elem in {elemIn, elemOut}
Base::A->entry(elemInIdxSLE, iElemIdxSLE) -= couplingUpdate1;
Base::A->entry(iElemIdxSLE, elemInIdxSLE) -= couplingUpdate1;
Base::A->entry(elemOutIdxSLE, iElemIdxSLE) -= couplingUpdate1;
Base::A->entry(iElemIdxSLE, elemOutIdxSLE) -= couplingUpdate1;
for (int i = 0; i < dim; i++)
{
const Scalar bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
if (i != dim - 1)
{
const Scalar interfaceUpdate =
weight * betaEvaluation * (iFrame[i] * qpGlobal);
const Scalar couplingUpdate3 =
(iFrame[i] * qpGlobal) * couplingUpdate1;
//quadrature for
// int_iElem beta * phi_iElem,0 * phi_iElem,i ds
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) += interfaceUpdate;
Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) += interfaceUpdate;
//quadrature for
// -int_iElem beta * phi_iElem,i * avg(phi_elem,0) ds
//for elem in {elemIn, elemOut} using a quadrature rule
Base::A->entry(elemInIdxSLE, iElemIdxSLE + i + 1) -= couplingUpdate3;
Base::A->entry(iElemIdxSLE + i + 1, elemInIdxSLE) -= couplingUpdate3;
Base::A->entry(elemOutIdxSLE, iElemIdxSLE + i + 1) -= couplingUpdate3;
Base::A->entry(iElemIdxSLE + i + 1, elemOutIdxSLE) -= couplingUpdate3;
}
//quadrature for
// int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds
//and
// int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds
//for elem1, elem2 in {elemIn, elemOut}
Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE) += bulkUpdate3;
Base::A->entry(elemInIdxSLE, elemInIdxSLE + i + 1) += bulkUpdate3;
Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE) += bulkUpdate3;
Base::A->entry(elemOutIdxSLE, elemOutIdxSLE + i + 1) += bulkUpdate3;
Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE) += bulkUpdate4;
Base::A->entry(elemOutIdxSLE, elemInIdxSLE + i + 1) += bulkUpdate4;
Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE) += bulkUpdate4;
Base::A->entry(elemInIdxSLE, elemOutIdxSLE + i + 1) += bulkUpdate4;
//quadrature for
// -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
//for elem in {elemIn, elemOut}
Base::A->entry(elemInIdxSLE + i + 1, iElemIdxSLE) -= couplingUpdate2;
Base::A->entry(iElemIdxSLE, elemInIdxSLE + i + 1) -= couplingUpdate2;
Base::A->entry(elemOutIdxSLE + i + 1, iElemIdxSLE) -= couplingUpdate2;
Base::A->entry(iElemIdxSLE, elemOutIdxSLE + i + 1) -= couplingUpdate2;
}
};
//lambda to evaluate the coupling terms
// int_iElem beta * phi_iElem,i * phi_iElem,j ds
//and
// int_iElem K_perp * jump(phi_elem1,k) * jump(phi_elem2,l) ds
//and
// int_iElem beta * avg(phi_elem1,k) * avg(phi_elem2,l) ds //and
// -int_iElem beta * phi_iElem,i * avg(phi_elem1,k) ds
//for i,j = 1,...,dim-1 and k,l = 1,...,dim
//and for elem1,elem2 in {elemIn, elemOut} using a quadrature rule
//TODO symmetrize more efficiently
const auto lambda_KperpPlus1 =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
const Scalar betaEvaluation = beta(qpGlobal);
const Scalar KEvaluation = Base::problem_.Kperp(qpGlobal);
const Scalar bulkUpdate1 =
weight * (KEvaluation + 0.25 * betaEvaluation);
const Scalar bulkUpdate2 =
weight * (-KEvaluation + 0.25 * betaEvaluation);
for (int i = 0; i < dim; i++)
{
const Scalar qpI = iFrame[i] * qpGlobal;
const Scalar interfaceUpdate1 = weight * betaEvaluation * qpI;
const Scalar couplingUpdate1 =
0.5 * weight * betaEvaluation * qpGlobal[i];
const Scalar couplingUpdate2 =
(iFrame[i] * qpGlobal) * couplingUpdate1;
const Scalar bulkUpdate3 = qpGlobal[i] * qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * qpGlobal[i] * bulkUpdate2;
if (i != dim - 1)
{
//quadrature for
// int_iElem beta * (phi_iElem,i)^2 ds
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate1 * qpI;
//quadrature for
// -int_iElem beta * phi_iElem,i * avg(phi_elem,i) ds
//for elem in {elemIn, elemOut}
Base::A->entry(elemInIdxSLE + i + 1, iElemIdxSLE + i + 1) -=
couplingUpdate2;
Base::A->entry(iElemIdxSLE + i + 1, elemInIdxSLE + i + 1) -=
couplingUpdate2;
Base::A->entry(elemOutIdxSLE + i + 1, iElemIdxSLE + i + 1) -=
couplingUpdate2;
Base::A->entry(iElemIdxSLE + i + 1, elemOutIdxSLE + i + 1) -=
couplingUpdate2;
}
//quadrature for
// int_iElem K_perp * jump(phi_elem1,i) * jump(phi_elem2,i) ds
//and
// int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,i) ds
//for elem1, elem2 in {elemIn, elemOut}
Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
bulkUpdate3;
Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
bulkUpdate3;
Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + i + 1) +=
bulkUpdate4;
Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + i + 1) +=
bulkUpdate4;
for (int j = 0; j < i; j++)
{
const Scalar couplingUpdate3 =
(iFrame[j] * qpGlobal) * couplingUpdate1;
const Scalar bulkUpdate5 =
qpGlobal[i] * qpGlobal[j] * bulkUpdate1;
const Scalar bulkUpdate6 =
qpGlobal[i] * qpGlobal[j] * bulkUpdate2;
if (i != dim - 1)
{
const Scalar interfaceUpdate2 =
interfaceUpdate1 * (iFrame[j] * qpGlobal);
const Scalar couplingUpdate4 =
0.5 * weight * betaEvaluation * qpGlobal[j] * qpI;
//quadrature for
// int_iElem beta * phi_iElem,i * phi_iElem,j ds
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate2;
Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate2;
//quadrature for
// -int_iElem beta * phi_iElem,i * avg(phi_elem,j) ds
//for elem in {elemIn, elemOut}
Base::A->entry(elemInIdxSLE + j + 1, iElemIdxSLE + i + 1) -=
couplingUpdate4;
Base::A->entry(iElemIdxSLE + i + 1, elemInIdxSLE + j + 1) -=
couplingUpdate4;
Base::A->entry(elemOutIdxSLE + j + 1, iElemIdxSLE + i + 1) -=
couplingUpdate4;
Base::A->entry(iElemIdxSLE + i + 1, elemOutIdxSLE + j + 1) -=
couplingUpdate4;
}
//quadrature for
// int_iElem K_perp * jump(phi_elem1,i) * jump(phi_elem2,j) ds
//and
// int_iElem beta * avg(phi_elem1,i) * avg(phi_elem2,j) ds
//for elem1, elem2 in {elemIn, elemOut}
Base::A->entry(elemInIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
bulkUpdate5;
Base::A->entry(elemInIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
bulkUpdate5;
Base::A->entry(elemOutIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
bulkUpdate5;
Base::A->entry(elemOutIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
bulkUpdate5;
Base::A->entry(elemInIdxSLE + i + 1, elemOutIdxSLE + j + 1) +=
bulkUpdate6;
Base::A->entry(elemOutIdxSLE + j + 1, elemInIdxSLE + i + 1) +=
bulkUpdate6;
Base::A->entry(elemOutIdxSLE + i + 1, elemInIdxSLE + j + 1) +=
bulkUpdate6;
Base::A->entry(elemInIdxSLE + j + 1, elemOutIdxSLE + i + 1) +=
bulkUpdate6;
//quadrature for
// -int_iElem beta * phi_iElem,j * avg(phi_elem,i) ds
//for elem in {elemIn, elemOut}
Base::A->entry(elemInIdxSLE + i + 1, iElemIdxSLE + j + 1) -=
couplingUpdate3;
Base::A->entry(iElemIdxSLE + j + 1, elemInIdxSLE + i + 1) -=
couplingUpdate3;
Base::A->entry(elemOutIdxSLE + i + 1, iElemIdxSLE + j + 1) -=
couplingUpdate3;
Base::A->entry(iElemIdxSLE + j + 1, elemOutIdxSLE + i + 1) -=
couplingUpdate3;
}
}
};
// === interface entries ===
//lambda to fill the load vector b with interface entries using a
//quadrature rule
const auto lambda_qInterface =
[&](const Coordinate& qpGlobal, const Scalar weight) {
const Scalar
qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal));
//quadrature for int_elem q*phi_elem,0 dV
Base::b[iElemIdxSLE] += weight * qInterfaceEvaluation;
//quadrature for int_elem q*phi_elem,i dV
for (int i = 0; i < dim - 1; i++)
{
Base::b[iElemIdxSLE + i + 1] +=
weight * (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
}
};
//lambda to evaluate
// int_iElem (Kparallel * grad(d phi_iElem,i)) * grad(d phi_iElem,j) ds
//using a quadrature rule for i,j = 0,...,dim-1
const auto lambda_Kpar =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
const CoordinateMatrix KEvaluation =
Base::problem_.Kparallel(qpGlobal);
const Scalar dEvaluation = Base::problem_.aperture(qpGlobal);
const Scalar dEvaluation2 = dEvaluation * dEvaluation;
const Coordinate gradDEvaluation =
Base::problem_.gradAperture(qpGlobal);
Coordinate KparDotTau_i; //storage for Kparallel * iFrame[i]
Coordinate KparDotGradD; //storage for Kparallel * grad(d)
KEvaluation.mv(gradDEvaluation, KparDotGradD);
const Scalar interfaceUpdate1 = KparDotGradD * gradDEvaluation;
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,0)) * grad(d * phi,iElem,0) ds
Base::A->entry(iElemIdxSLE, iElemIdxSLE) += weight * interfaceUpdate1;
for (int i = 0; i < dim - 1 ; i++)
{
KEvaluation.mv(iFrame[i], KparDotTau_i);
const Scalar qpI = qpGlobal * iFrame[i];
const Scalar interfaceUpdate2 =
dEvaluation * (KparDotGradD * iFrame[i]);
const Scalar interfaceUpdate3 =
weight * (qpI * interfaceUpdate1 + interfaceUpdate2);
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i))
// * grad(d * phi,iElem,0) ds
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) += interfaceUpdate3;
Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) += interfaceUpdate3;
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i))
// * grad(d * phi,iElem,i) ds
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
weight * ( dEvaluation2 * (KparDotTau_i * iFrame[i])
+ qpI * (qpI * interfaceUpdate1 + 2.0 * interfaceUpdate2) );
for (int j = 0; j < i; j++)
{
const Scalar interfaceUpdate4 =
weight * ( dEvaluation2 * (KparDotTau_i * iFrame[j])
+ dEvaluation * qpI * (KparDotGradD * iFrame[j])
+ (qpGlobal * iFrame[j])
* (interfaceUpdate2 + qpI * interfaceUpdate1) );
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i))
// * grad(d * phi,iElem,j) ds
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate4;
Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate4;
}
}
};
//apply quadrature rules on iElem
Base::applyQuadratureRule(iGeo, rule_Kperp, lambda_Kperp);
Base::applyQuadratureRule(iGeo, rule_KperpPlus1, lambda_KperpPlus1);
Base::applyQuadratureRule(iGeo, rule_qInterface, lambda_qInterface);
Base::applyQuadratureRule(iGeo, rule_Kpar, lambda_Kpar);
//iterate over all intersection with the boundary of iElem
for (const auto& intersct : intersections(iGridView_, iElem))
{
const auto& normal = intersct.centerUnitOuterNormal();
const auto& intersctGeo = intersct.geometry();
//determine penalty parameter
const Scalar mu = getPenaltyParameter(mu0, intersct, iGeoCenter);
if (intersct.neighbor()) //interior interface edge
{
//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 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
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);
Coordinate KparDotTau_i; //storage for Kparallel * iFrame[i]
Coordinate KparDotTau_j; //storage for Kparallel * iFrame[j]
Coordinate KparDotGradD; //storage for Kparallel * grad(d)
KEvaluation.mv(Base::problem_.gradAperture(qpGlobal),
KparDotGradD);
const Scalar interfaceUpdate1 = KparDotGradD * normal;
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
weight * dEvaluation * (mu * dEvaluation - interfaceUpdate1);
for (int i = 0; i < dim - 1; i++)
{
KEvaluation.mv(iFrame[i], KparDotTau_i);
const Scalar qpI = qpGlobal * iFrame[i];
const Scalar interfaceUpdate2 = mu * dEvaluation * qpI;
const Scalar interfaceUpdate3 =
dEvaluation * (KparDotTau_i * normal);
const Scalar interfaceUpdate4 = weight * dEvaluation * (interfaceUpdate2 - 0.5 * interfaceUpdate3
- qpI * interfaceUpdate1);
//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
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
interfaceUpdate4;
Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
interfaceUpdate4;
//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
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
weight * dEvaluation * qpI * (interfaceUpdate2
- interfaceUpdate3 - qpI * interfaceUpdate1);
for (int j = 0; j < i; j++)
{
KEvaluation.mv(iFrame[j], KparDotTau_j);
const Scalar interfaceUpdate5 = weight * dEvaluation *
( (qpGlobal * iFrame[j]) * ( interfaceUpdate2
- 0.5 * interfaceUpdate3 - qpI * interfaceUpdate1 )
- 0.5 * (KparDotTau_j * normal) * qpI * dEvaluation );
//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
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate5;
Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate5;
}
}
if (neighborIdx > iElemIdx)
{ //make sure that each interface edge is considered only once
return;
}
const int neighborIdxSLE =
(dim + 1) * Base::gridView_.size(0) + dim * neighborIdx;
const InterfaceBasis neighborFrame =
interfaceFrame(Base::gridView_.grid().asIntersection(
intersct.outside()).centerUnitOuterNormal());
//storage for K_par * neighborFrame[i]
Coordinate KparDotNeighborTau_i;
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_neighbor,0) dr
Base::A->entry(iElemIdxSLE, neighborIdxSLE) -=
weight * mu * dEvaluation2;
Base::A->entry(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);
//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
Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE) +=
interfaceUpdate3a;
Base::A->entry(neighborIdxSLE, iElemIdxSLE + i + 1) +=
interfaceUpdate3a;
//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
Base::A->entry(neighborIdxSLE + i + 1, iElemIdxSLE) +=
interfaceUpdate3b;
Base::A->entry(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
Base::A->entry(iElemIdxSLE + i + 1, neighborIdxSLE + j + 1) +=
interfaceUpdate5;
Base::A->entry(neighborIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate5;
}
}
};
Base::applyQuadratureRule(intersctGeo, rule_Kpar_iFacet,
lambda_Kpar_iFacet);
}
else //boundary interface edge
{
//lambda to evaluate the interface terms on the boundary
// 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
//for i,j = 0,...,dim-1 using a quadrature rule
const auto lambda_Kpar_iFacet_boundary =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
const Scalar dEvaluation = Base::problem_.aperture(qpGlobal);
const CoordinateMatrix KEvaluation =
Base::problem_.Kparallel(qpGlobal);
Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i]
Coordinate KparDotTau_j; //storage for K_parallel * iFrame[j]
Coordinate KparDotGradD; //storage for K_parallel * grad(d)
KEvaluation.mv(Base::problem_.gradAperture(qpGlobal),
KparDotGradD);
const Scalar interfaceUpdate1 = 2.0 * (KparDotGradD * normal);
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
weight * dEvaluation * (mu * dEvaluation - interfaceUpdate1);
for (int i = 0; i < dim - 1; i++)
{
KEvaluation.mv(iFrame[i], KparDotTau_i);
const Scalar qpI = qpGlobal * iFrame[i];
const Scalar interfaceUpdate2 = mu * dEvaluation * qpI;
const Scalar interfaceUpdate3 =
dEvaluation * (KparDotTau_i * normal);
const Scalar interfaceUpdate4 = weight * dEvaluation *
(interfaceUpdate2 - interfaceUpdate3
- qpI * interfaceUpdate1);
//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->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
interfaceUpdate4;
Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
interfaceUpdate4;
//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->entry(iElemIdxSLE + i + 1, iElemIdxSLE + i + 1) +=
weight * qpI * dEvaluation *
(interfaceUpdate2 - 2.0 * interfaceUpdate3
- qpI * interfaceUpdate1);
for (int j = 0; j < i; j++)
{
KEvaluation.mv(iFrame[j], KparDotTau_j);
const Scalar qpJ = qpGlobal * iFrame[j];
const Scalar interfaceUpdate5 =
weight * dEvaluation * (qpJ * (interfaceUpdate2
- interfaceUpdate3 - qpI * interfaceUpdate1)
- (KparDotTau_j * normal) * qpI * dEvaluation);
//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->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate5;
Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
interfaceUpdate5;
}
}
};
//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
const auto lambda_interfaceBoundary =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
const Scalar dEvaluation = Base::problem_.aperture(qpGlobal);
const CoordinateMatrix KEvaluation =
Base::problem_.Kparallel(qpGlobal);
Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i]
Coordinate KparDotGradD; //storage for K_parallel * grad(d)
KEvaluation.mv(Base::problem_.gradAperture(qpGlobal),
KparDotGradD);
const Scalar interfaceUpdate1 = weight * dEvaluation *
Base::problem_.interfaceBoundary(qpGlobal);
const Scalar interfaceUpdate2 =
mu * dEvaluation - KparDotGradD * normal;
Base::b[iElemIdxSLE] += interfaceUpdate1 * interfaceUpdate2;
for (int i = 0; i < dim - 1; i++)
{
KEvaluation.mv(iFrame[i], KparDotTau_i);
Base::b[iElemIdxSLE + i + 1] +=
interfaceUpdate1 * ((qpGlobal * iFrame[i]) * interfaceUpdate2
- dEvaluation * (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);
}
}
}
}
//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] ]
//3d: tau_1 ~ [ normal[1] + normal[2], tau_2 ~ [ -normal[1],
// -normal[0], normal[0] + normal[2]
// -normal[0] ], -normal[1] ]
InterfaceBasis interfaceFrame(const Coordinate& normal) const
{
InterfaceBasis iBasis;
for (int i = 0; i < dim - 1; i++)
{
Coordinate tau(0.0);
for (int j = 0; j < dim; j++)
{ if (j == i)
{
continue;
}
tau[i] += normal[j];
tau[j] = -normal[i];
}
if constexpr (dim != 2)
{
tau /= tau.two_norm();
}
iBasis[i] = tau;
}
return iBasis;
}
void writeVTKoutput () const
{
//degrees of freedom of bulk and interface grids
const int bulkDOF = (dim + 1) * Base::gridView_.size(0);
const int interfaceDOF = dim * iGridView_.size(0);
//=== write bulk solution to file ===
Base::writeVTKoutput(bulkDOF);
//=== write interface solution to file ===
//storage for pressure data, for each interface element we store the
//pressure at the corners of the interface element, the VTKFunction will
//create the output using a linear interpolation for each interface element
Dune::DynamicVector<Scalar> iPressure(interfaceDOF, 0.0);
Dune::DynamicVector<Scalar> exactIPressure(interfaceDOF, 0.0);
for (const auto& iElem : elements(iGridView_))
{
const int iElemIdxOutput = dim * iMapper_.index(iElem);
const int iElemIdxSLE = bulkDOF + iElemIdxOutput;
const auto& iGeo = iElem.geometry();
//get basis of the interface coordinate system for evaluation of
//basis function phi_iElem,i with i=1,...,dim-1
const InterfaceBasis iFrame = interfaceFrame(
Base::gridView_.grid().asIntersection(iElem).centerUnitOuterNormal());
for (int k = 0; k < iGeo.corners(); k++)
{
if (Base::problem_.hasExactSolution())
{
exactIPressure[iElemIdxOutput + k] =
Base::problem_.exactInterfaceSolution(iGeo.corner(k));
}
//contribution of the basis function
// phi_iElem,0 (x) = indicator(iElem);
//at the kth corner of iElem
iPressure[iElemIdxOutput + k] = Base::d[iElemIdxSLE];
for (int i = 0; i < dim - 1; i++)
{
//contribution of the basis function
// phi_iElem,i (x) = (x * tau[i]) * indicator(iElem);
//at the kth corner of iElem
iPressure[iElemIdxOutput + k] +=
Base::d[iElemIdxSLE + i + 1] * (iGeo.corner(k) * iFrame[i]);
}
} }
Dune::VTKWriter<InterfaceGridView>
vtkWriter(iGridView_, Dune::VTK::nonconforming);
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(iGridView_, iPressure, "interfacePressure")));
if (Base::problem_.hasExactSolution())
{
vtkWriter.addVertexData( std::shared_ptr<const VTKFunction>(
new P1Function(iGridView_, exactIPressure, "exactInterfacePressure")));
}
vtkWriter.pwrite("interfacePressureData", "", "data");
}
//returns the interface penalty parameter mu^Gamma for an interface facet
//given as intersection of the interface grid where center is the center of
//intersection.inside()
const Scalar getPenaltyParameter (const Scalar mu0, const auto& intersection,
const auto& center) const
{
Scalar mu = Base::getLargestEigenvalue(Base::problem_.Kparallel(center))
/ Base::getDiameter(intersection.inside());
if (intersection.neighbor())
{
const auto& neighbor = intersection.outside();
mu = std::max(mu, Base::getLargestEigenvalue(
Base::problem_.Kparallel(neighbor.geometry().center()))
/ Base::getDiameter(neighbor));
}
return mu * mu0 * 2.0 * (1.0 + dim);
}
};
#endif