Skip to content
Snippets Groups Projects
Commit 29afbe70 authored by Hörl, Maximilian's avatar Hörl, Maximilian
Browse files

use quadrature rules for interface facets in mmdg

parent 24d2a897
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 270 additions and 158 deletions
dune/mmdg/mmdg.hh
+ 270
− 158
View file @ 29afbe70
......@@ -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,57 +471,86 @@ 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);
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 //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);
//evaluate the interface term
// int_intersct mu^Gamma * d^2 * jump(phi_iElem,0)
Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i]
Coordinate KparDotTau_j; //storage for K_parallel * iFrame[j]
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
(*Base::A)[iElemIdxSLE][iElemIdxSLE] += mu * dEvaluation2;
(*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
weight * 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
//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;
KEvaluation.mv(iFrame[i], KparDotTau_i);
const Scalar qpI = qpGlobal * iFrame[i];
const Scalar interfaceUpdate1 = mu * qpI;
const Scalar interfaceUpdate2 = KparDotTau_i * normal;
const Scalar interfaceUpdate3 =
dEvaluation2 * (interfaceUpdate1 - 0.5 * interfaceUpdate2);
const Scalar interfaceUpdate3 = weight * dEvaluation2
* (interfaceUpdate1 - 0.5 * interfaceUpdate2);
//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;
//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] +=
dEvaluation2 * centerI * (interfaceUpdate1 - interfaceUpdate2);
weight * dEvaluation2 * qpI *
(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 = dEvaluation2 *
((center * iFrame[j]) * interfaceUpdate1
- 0.5 * (interfaceUpdate2 * (center * iFrame[j])
+ (KparDotTau_j * normal) * centerI));
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] +=
......@@ -521,20 +558,13 @@ private:
}
}
//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;
return;
}
const int neighborIdxSLE =
(dim + 1) * Base::gridView_.size(0) + dim * neighborIdx;
//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());
......@@ -542,38 +572,47 @@ private:
//storage for K_par * neighborFrame[i]
Coordinate KparDotNeighborTau_i;
//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;
//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++)
{
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;
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 =
dEvaluation2 * (-interfaceUpdate1a + 0.5 * interfaceUpdate2a);
const Scalar interfaceUpdate3b =
dEvaluation2 * (-interfaceUpdate1b - 0.5 * interfaceUpdate2b);
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 //TODO: terms with grad(d)
(*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] +=
interfaceUpdate3a;
(*Base::A)[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 //TODO: terms with grad(d)
(*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] +=
interfaceUpdate3b;
(*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] +=
......@@ -583,73 +622,146 @@ private:
{
Coordinate KparDotNeighborTau_j;
Base::problem_.Kparallel(center).mv(
neighborFrame[j], KparDotNeighborTau_j);
KEvaluation.mv(neighborFrame[j], KparDotNeighborTau_j);
const Scalar interfaceUpdate5 =
-dEvaluation2 * ((center * neighborFrame[j]) * interfaceUpdate1a
+ 0.5 * (-interfaceUpdate2a * (center * neighborFrame[j])
+ (KparDotNeighborTau_j * normal) * centerI));
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);
// -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);
Base::b[iElemIdxSLE] += mu * d2gGamma;
Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i]
Coordinate KparDotTau_j; //storage for K_parallel * iFrame[j]
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
(*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
weight * mu * dEvaluation2;
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;
KEvaluation.mv(iFrame[i], KparDotTau_i);
const Scalar qpI = qpGlobal * iFrame[i];
const Scalar interfaceUpdate1 = mu * qpI;
const Scalar interfaceUpdate2 = KparDotTau_i * normal;
const Scalar interfaceUpdate3 =
dEvaluation2 * (interfaceUpdate1 - interfaceUpdate2);
Base::b[iElemIdxSLE + i + 1] +=
d2gGamma * (interfaceUpdate1 - interfaceUpdate2);
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] +=
centerI * dEvaluation2 *
weight * qpI * 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));
{
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);
}
}
}
......
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please to comment