Select Git revision
mmdg.hh 27.17 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);
using Base = DG<GridView, Mapper, Problem>;
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>;
using InterfaceBasis = std::array<Coordinate, dim-1>;
using QRuleInterface = typename Base::QRuleEdge;
using QRulesInterface = typename Base::QRulesEdge;
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 mu, const Scalar xi)
{
assembleSLE(mu, 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
{
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
for (const auto& ip : qrule)
{
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
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 += quadratureFactor * errorEvaluation * errorEvaluation;
}
}
return sqrt(bulkError * bulkError + interfaceError);
}
const InterfaceGridView& iGridView_;
const InterfaceMapper& iMapper_;
private:
//assemble stiffness matrix A and load vector b
void assembleSLE (const Scalar mu, const Scalar xi)
{
//modell parameters //NOTE: what should be the capture of alpha?
const auto alpha = [this](const Coordinate& pos) -> Scalar
{
return Base::problem_.Kperp(pos) / Base::problem_.aperture(pos);
};
const auto beta = [&alpha, &xi](const Coordinate& pos) -> Scalar
{
return 4 * alpha(pos) / (2 * xi - 1);
};
//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);
//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(mu);
for (const auto& iElem : elements(iGridView_))
{
const int iElemIdx = iMapper_.index(iElem);
const auto& iGeo = iElem.geometry();
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;
// === coupling entries ===
//evaluate the coupling term
// int_iElem beta * phi_iElem,i * phi_iElem,j ds
//using a quadrature rule for i = 0 and j = 0,...,dim-1 or vice versa //TODO symmetrize more efficiently
for (const auto& ip : rule_Kperp)
{
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar betaEvaluation = beta(qpGlobal);
Base::A->entry(iElemIdxSLE, iElemIdxSLE) +=
quadratureFactor * betaEvaluation;
for (int i = 0; i < dim - 1; i++)
{
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE) +=
quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
Base::A->entry(iElemIdxSLE, iElemIdxSLE + i + 1) +=
quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal);
}
}
//evaluate the coupling term
// int_iElem beta * phi_iElem,i * phi_iElem,j ds
//using a quadrature rule for i,j = 1,...,dim-1 //TODO symmetrize more efficiently
for (const auto& ip : rule_KperpPlus1)
{
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar betaEvaluation = beta(qpGlobal);
for (int i = 0; i < dim - 1; i++)
{
for (int j = 0; j <= i; j++)
{
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
* (iFrame[j] * qpGlobal);
if (i != j)
{ Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
quadratureFactor * betaEvaluation * (iFrame[i] * qpGlobal)
* (iFrame[j] * qpGlobal);
}
}
}
}
//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);
//evalute 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
//for elem1, elem2 in {elemIn, elemOut} and i = 0,...,dim
//using a quadrature rule
for (const auto& ip : rule_Kperp)
{
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar alphaEvaluation = alpha(qpGlobal);
const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
const Scalar bulkUpdate1 =
quadratureFactor * (alphaEvaluation + betaEvaluation4);
const Scalar bulkUpdate2 =
quadratureFactor * (-alphaEvaluation + betaEvaluation4);
Base::A->entry(elemInIdxSLE, elemInIdxSLE) += bulkUpdate1;
Base::A->entry(elemOutIdxSLE, elemOutIdxSLE) += bulkUpdate1;
Base::A->entry(elemInIdxSLE, elemOutIdxSLE) += bulkUpdate2;
Base::A->entry(elemOutIdxSLE, elemInIdxSLE) += bulkUpdate2;
for (int i = 0; i < dim; i++)
{
const Scalar bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
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;
}
}
//evalute 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
//for elem1, elem2 in {elemIn, elemOut} and i,j = 1,...,dim
//using a quadrature rule
for (const auto& ip : rule_KperpPlus1)
{
const auto& qpLocal = ip.position();
const auto& qpGlobal = iGeo.global(qpLocal);
const Scalar quadratureFactor =
ip.weight() * iGeo.integrationElement(qpLocal);
const Scalar alphaEvaluation = alpha(qpGlobal);
const Scalar betaEvaluation4 = 0.25 * beta(qpGlobal);
const Scalar bulkUpdate1 =
quadratureFactor * (alphaEvaluation + betaEvaluation4);
const Scalar bulkUpdate2 =
quadratureFactor * (-alphaEvaluation + betaEvaluation4);
for (int i = 0; i < dim; i++)
{
const Scalar bulkUpdate3 = qpGlobal[i] * qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * qpGlobal[i] * bulkUpdate2;
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 bulkUpdate5 = qpGlobal[i] * qpGlobal[j] * bulkUpdate1;
const Scalar bulkUpdate6 = qpGlobal[i] * qpGlobal[j] * bulkUpdate2;
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;
}
}
}
//evalute 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
for (const auto& ip : rule_Kperp) //TODO: symmetrize more efficiently
//TODO: apply quadrature rule only once
{
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);
const Scalar couplingUpdate1 = quadratureFactor * betaEvaluation2;
Base::A->entry(elemInIdxSLE, iElemIdx) -= couplingUpdate1;
Base::A->entry(iElemIdx, elemInIdxSLE) -= couplingUpdate1;
Base::A->entry(elemOutIdxSLE, iElemIdx) -= couplingUpdate1; Base::A->entry(iElemIdx, elemOutIdxSLE) -= couplingUpdate1;
for (int i = 0; i < dim - 1; i++)
{
const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
const Scalar couplingUpdate3 =
(iFrame[i] * qpGlobal) * couplingUpdate1;
Base::A->entry(elemInIdxSLE + i + 1, iElemIdx) -= couplingUpdate2;
Base::A->entry(iElemIdx, elemInIdxSLE + i + 1) -= couplingUpdate2;
Base::A->entry(elemOutIdxSLE + i + 1, iElemIdx) -= couplingUpdate2;
Base::A->entry(iElemIdx, elemOutIdxSLE + i + 1) -= couplingUpdate2;
Base::A->entry(elemInIdxSLE, iElemIdx + i + 1) -= couplingUpdate3;
Base::A->entry(iElemIdx + i + 1, elemInIdxSLE) -= couplingUpdate3;
Base::A->entry(elemOutIdxSLE, iElemIdx + i + 1) -= couplingUpdate3;
Base::A->entry(iElemIdx + i + 1, elemOutIdxSLE) -= couplingUpdate3;
}
//i = dim - 1
const Scalar couplingUpdate4 = qpGlobal[dim] * couplingUpdate1;
Base::A->entry(elemInIdxSLE + dim + 1, iElemIdx) -= couplingUpdate4;
Base::A->entry(iElemIdx, elemInIdxSLE + dim + 1) -= couplingUpdate4;
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] ]
//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 - 1; 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 iElemIdxSLE = bulkDOF + dim * iMapper_.index(iElem);
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[iElemIdxSLE + 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[iElemIdxSLE + 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[iElemIdxSLE + 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");
}
};
#endif