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

merge with back2sparse

parent 28f4d849
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 280 deletions
dune/mmdg/mmdg.hh
+ 0
− 280
View file @ a1fa6547
...@@ -172,54 +172,6 @@ private: ...@@ -172,54 +172,6 @@ private:
const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn); const int elemInIdxSLE = (dim + 1) * Base::mapper_.index(elemIn);
const int elemOutIdxSLE = (dim + 1) * Base::mapper_.index(elemOut); 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 === // === coupling entries ===
//lambda to evaluate the coupling terms //lambda to evaluate the coupling terms
...@@ -228,7 +180,6 @@ private: ...@@ -228,7 +180,6 @@ private:
// int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds // int_iElem K_perp * jump(phi_elem1,0) * jump(phi_elem2,i) ds
//and //and
// int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds // int_iElem beta * avg(phi_elem1,0) * avg(phi_elem2,i) ds
>>>>>>> back2NonSparse
//and //and
// -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds // -int_iElem beta * phi_iElem,0 * avg(phi_elem,i) ds
//and //and
...@@ -272,49 +223,9 @@ private: ...@@ -272,49 +223,9 @@ private:
for (int i = 0; i < dim; i++) 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 bulkUpdate3 = qpGlobal[i] * bulkUpdate1;
const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2; const Scalar bulkUpdate4 = qpGlobal[i] * bulkUpdate2;
const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1; const Scalar couplingUpdate2 = qpGlobal[i] * couplingUpdate1;
>>>>>>> back2NonSparse
if (i != dim - 1) if (i != dim - 1)
{ {
...@@ -511,22 +422,6 @@ private: ...@@ -511,22 +422,6 @@ private:
const Scalar const Scalar
qInterfaceEvaluation(Base::problem_.qInterface(qpGlobal)); 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 //quadrature for int_elem q*phi_elem,0 dV
Base::b[iElemIdxSLE] += weight * qInterfaceEvaluation; Base::b[iElemIdxSLE] += weight * qInterfaceEvaluation;
...@@ -537,7 +432,6 @@ private: ...@@ -537,7 +432,6 @@ private:
weight * (qpGlobal * iFrame[i]) * qInterfaceEvaluation; weight * (qpGlobal * iFrame[i]) * qInterfaceEvaluation;
} }
}; };
>>>>>>> back2NonSparse
//lambda to evaluate //lambda to evaluate
// int_iElem (Kparallel * grad(d phi_iElem,i)) * grad(d phi_iElem,j) ds // int_iElem (Kparallel * grad(d phi_iElem,i)) * grad(d phi_iElem,j) ds
...@@ -545,49 +439,6 @@ private: ...@@ -545,49 +439,6 @@ private:
const auto lambda_Kpar = const auto lambda_Kpar =
[&](const Coordinate& qpGlobal, const Scalar weight) [&](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 = const CoordinateMatrix KEvaluation =
Base::problem_.Kparallel(qpGlobal); Base::problem_.Kparallel(qpGlobal);
const Scalar dEvaluation = Base::problem_.aperture(qpGlobal); const Scalar dEvaluation = Base::problem_.aperture(qpGlobal);
...@@ -607,7 +458,6 @@ private: ...@@ -607,7 +458,6 @@ private:
Base::A->entry(iElemIdxSLE, iElemIdxSLE) += weight * interfaceUpdate1; Base::A->entry(iElemIdxSLE, iElemIdxSLE) += weight * interfaceUpdate1;
for (int i = 0; i < dim - 1 ; i++) for (int i = 0; i < dim - 1 ; i++)
>>>>>>> back2NonSparse
{ {
KEvaluation.mv(iFrame[i], KparDotTau_i); KEvaluation.mv(iFrame[i], KparDotTau_i);
...@@ -615,11 +465,7 @@ private: ...@@ -615,11 +465,7 @@ private:
const Scalar interfaceUpdate2 = const Scalar interfaceUpdate2 =
dEvaluation * (KparDotGradD * iFrame[i]); dEvaluation * (KparDotGradD * iFrame[i]);
const Scalar interfaceUpdate3 = const Scalar interfaceUpdate3 =
<<<<<<< HEAD
d2Evaluation * (interfaceUpdate1 - 0.5 * interfaceUpdate2);
=======
weight * (qpI * interfaceUpdate1 + interfaceUpdate2); weight * (qpI * interfaceUpdate1 + interfaceUpdate2);
>>>>>>> back2NonSparse
//quadrature for //quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i)) // int_iElem (Kpar * grad(d * phi_iElem,i))
...@@ -634,17 +480,6 @@ private: ...@@ -634,17 +480,6 @@ private:
weight * ( dEvaluation2 * (KparDotTau_i * iFrame[i]) weight * ( dEvaluation2 * (KparDotTau_i * iFrame[i])
+ qpI * (qpI * interfaceUpdate1 + 2.0 * interfaceUpdate2) ); + 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++) for (int j = 0; j < i; j++)
{ {
const Scalar interfaceUpdate4 = const Scalar interfaceUpdate4 =
...@@ -656,7 +491,6 @@ private: ...@@ -656,7 +491,6 @@ private:
//quadrature for //quadrature for
// int_iElem (Kpar * grad(d * phi_iElem,i)) // int_iElem (Kpar * grad(d * phi_iElem,i))
// * grad(d * phi,iElem,j) ds // * grad(d * phi,iElem,j) ds
>>>>>>> back2NonSparse
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) += Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
interfaceUpdate4; interfaceUpdate4;
Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) += Base::A->entry(iElemIdxSLE + j + 1, iElemIdxSLE + i + 1) +=
...@@ -677,119 +511,6 @@ private: ...@@ -677,119 +511,6 @@ private:
const auto& normal = intersct.centerUnitOuterNormal(); const auto& normal = intersct.centerUnitOuterNormal();
const auto& intersctGeo = intersct.geometry(); 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 //determine penalty parameter
const Scalar mu = getPenaltyParameter(mu0, intersct, iGeoCenter); const Scalar mu = getPenaltyParameter(mu0, intersct, iGeoCenter);
...@@ -1096,7 +817,6 @@ private: ...@@ -1096,7 +817,6 @@ private:
lambda_Kpar_iFacet_boundary); lambda_Kpar_iFacet_boundary);
Base::applyQuadratureRule(intersctGeo, rule_interfaceBoundary, Base::applyQuadratureRule(intersctGeo, rule_interfaceBoundary,
lambda_interfaceBoundary); lambda_interfaceBoundary);
>>>>>>> back2NonSparse
} }
} }
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment