From a1fa654790587f2c6d1cfb5e8d9c98b03a3c2554 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Wed, 1 Apr 2020 20:49:02 +0200
Subject: [PATCH] merge with back2sparse
---
dune/mmdg/mmdg.hh | 280 ----------------------------------------------
1 file changed, 280 deletions(-)
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 6b7c131..73890ca 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -172,54 +172,6 @@ private:
const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn);
const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut);
-<<<<<<< HEAD
- //evaluate 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;
- }
- }
-
- //evaluate the coupling terms
- // int_iElem alpha * jump(phi_elem1,i) * jump(phi_elem2,j) ds
-=======
// === coupling entries ===
//lambda to evaluate the coupling terms
@@ -228,7 +180,6 @@ private:
// 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
->>>>>>> back2NonSparse
//and
// -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
//and
@@ -272,49 +223,9 @@ private:
for (int i = 0; i < dim; i++)
{
-<<<<<<< HEAD
- 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;
- }
- }
- }
-
- //evaluate 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 bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
->>>>>>> back2NonSparse
if (i != dim - 1)
{
@@ -511,22 +422,6 @@ private:
const Scalar
qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal));
-<<<<<<< HEAD
- //evaluate
- // int_iElem (Kparallel * grad(d phi_iElem,i)) * grad(d phi_iElem,j) ds
- // = int_iElem d^2 * (Kpar * tau_i) * tau_j ds +//TODO terms with grad(d)
- //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);
- const auto KEvaluation = Base::problem_.Kparallel(qpGlobal);
- Scalar d2Evaluation = Base::problem_.aperture(qpGlobal);
- d2Evaluation *= d2Evaluation;
-=======
//quadrature for int_elem q*phi_elem,0 dV
Base::b[iElemIdxSLE] += weight * qInterfaceEvaluation;
@@ -537,7 +432,6 @@ private:
weight * (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
}
};
->>>>>>> back2NonSparse
//lambda to evaluate
// int_iElem (Kparallel * grad(d phi_iElem,i)) * grad(d phi_iElem,j) ds
@@ -545,49 +439,6 @@ private:
const auto lambda_Kpar =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
-<<<<<<< HEAD
- Coordinate KparDotTau_i;
- KEvaluation.mv(iFrame[i], KparDotTau_i);
-
- for (int j = 0; j <= i; j++)
- {
- Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
- quadratureFactor * d2Evaluation * ( KparDotTau_i * iFrame[j] );
-
- if (i != j)
- {
- Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
- quadratureFactor * d2Evaluation * ( 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;
- const Scalar dEvaluation = Base::problem_.aperture(center);
- const Scalar d2Evaluation = dEvaluation * dEvaluation;
-
- //evaluate the interface term
- // int_intersct mu^Gamma * d^2 * jump(phi_iElem,0)
- // * jump(phi_iElem,0) dr
- Base::A->entry(iElemIdxSLE, iElemIdxSLE) += mu * d2Evaluation;
-
- if (intersct.neighbor()) //inner interface edge
- {
- //evaluate the interface terms
- // int_intersct mu^Gamma * d^2 * jump(phi_iElem,i) * jump(phi_iElem,j) dr
- //and
- // -int_intersct d * avg(Kparallel * grad(d * phi_iElem,i))
- // * jump(phi_iElem,j) dr //TODO term with grad(d)
- //for i,j = 0,...,dim-1 and (i != 0 or j != 0)
- for (int i = 0; i < dim - 1; i++)
-=======
const CoordinateMatrix KEvaluation =
Base::problem_.Kparallel(qpGlobal);
const Scalar dEvaluation = Base::problem_.aperture(qpGlobal);
@@ -607,7 +458,6 @@ private:
Base::A->entry(iElemIdxSLE, iElemIdxSLE) += weight * interfaceUpdate1;
for (int i = 0; i < dim - 1 ; i++)
->>>>>>> back2NonSparse
{
KEvaluation.mv(iFrame[i], KparDotTau_i);
@@ -615,11 +465,7 @@ private:
const Scalar interfaceUpdate2 =
dEvaluation * (KparDotGradD * iFrame[i]);
const Scalar interfaceUpdate3 =
-<<<<<<< HEAD
- d2Evaluation * (interfaceUpdate1 - 0.5 * interfaceUpdate2);
-=======
weight * (qpI * interfaceUpdate1 + interfaceUpdate2);
->>>>>>> back2NonSparse
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i))
@@ -634,17 +480,6 @@ private:
weight * ( dEvaluation2 * (KparDotTau_i * iFrame[i])
+ qpI * (qpI * interfaceUpdate1 + 2.0 * interfaceUpdate2) );
-<<<<<<< HEAD
- 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 = d2Evaluation *
- ((center * iFrame[j]) * interfaceUpdate1
- - 0.5 * (interfaceUpdate2 * (center * iFrame[j])
- + (KparDotTau_j * normal) * centerI));
-
-=======
for (int j = 0; j < i; j++)
{
const Scalar interfaceUpdate4 =
@@ -656,7 +491,6 @@ private:
//quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i))
// * grad(d * phi,iElem,j) ds
->>>>>>> back2NonSparse
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate4;
Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
@@ -677,119 +511,6 @@ private:
const auto& normal = intersct.centerUnitOuterNormal();
const auto& intersctGeo = intersct.geometry();
-<<<<<<< HEAD
- //evaluate the interface terms
- // int_intersct mu^Gamma * 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
- //for i,j = 0,...,dim-1 //TODO terms with grad(d)
- 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 =
- d2Evaluation * (-interfaceUpdate1 + 0.5 * interfaceUpdate2);
- const Scalar interfaceUpdate4 =
- d2Evaluation * (-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) +=
- -d2Evaluation * 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 =
- d2Evaluation * (-interfaceUpdate5 - 0.5 * interfaceUpdate6);
- const Scalar interfaceUpdate8 =
- d2Evaluation * (-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 * d^2 phi_iElem,i * phi_iElem,j dr
- //and
- // -int_intersct d * phi_iElem,j * (Kparallel * grad(d * 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 * d * phi_iElem,i +
- // (Kparallel * grad(d * phi_i,iElem)) * normal) dr
- //for i = 0,...,dim-1 //TODO: terms with grad(d)
- const Scalar dgGamma = Base::problem_.aperture(center)
- * Base::problem_.boundary(center);
- Base::b[iElemIdxSLE] += mu * d2Evaluation * 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 =
- d2Evaluation * (interfaceUpdate1 - interfaceUpdate2);
-
- Base::b[iElemIdxSLE] +=
- dgGamma * d2Evaluation * (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) +=
- d2Evaluation * 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 = d2Evaluation *
- ((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;
- }
- }
-=======
//determine penalty parameter
const Scalar mu = getPenaltyParameter(mu0, intersct, iGeoCenter);
@@ -1096,7 +817,6 @@ private:
lambda_Kpar_iFacet_boundary);
Base::applyQuadratureRule(intersctGeo, rule_interfaceBoundary,
lambda_interfaceBoundary);
->>>>>>> back2NonSparse
}
}
}
--
GitLab