From 29afbe70fd8805b9d670f67d9faaf3fabe0ea101 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Mon, 30 Mar 2020 13:11:30 +0200
Subject: [PATCH] use quadrature rules for interface facets in mmdg
---
dune/mmdg/mmdg.hh | 428 +++++++++++++++++++++++++++++-----------------
1 file changed, 270 insertions(+), 158 deletions(-)
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 8b05f05..b6d5c22 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);
}
}
}
--
GitLab