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

add interface entries to mmdg

parent 420ab854
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
dune/mmdg/mmdg.hh
+ 236
− 83
View file @ 53ab1652
......@@ -11,7 +11,7 @@ public:
static constexpr int dim = GridView::dimension;
using Base = DG<Grid, GridView, Mapper, Problem>;
using Scalar = typename Base::Scalar;
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>;
......@@ -105,60 +105,6 @@ private:
//basis functions (0, phi_iElem,i)
const int iElemIdxSLE = (dim + 1)*Base::gridView_.size(0) + dim*iElemIdx;
// === 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;
}
}
//TODO: fill stiffness matrix A with interface entries
//evaluate
// int_iElem d * (K_par*grad_par(phi_iElem,i)) * grad_par(phi_iElem,j) ds
// = int_iElem d * (K_par * 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 K_parDotTau_i;
Base::problem_.Kparallel(qpGlobal).mv(iFrame[i], K_parDotTau_i);
for (int j = 0; j <= i; j++) //TODO: symmetrize
{
Base::A->entry(iElemIdxSLE + i + 1, iElemIdxSLE + j + 1) +=
quadratureFactor * Base::problem_.aperture(qpGlobal)
* ( K_parDotTau_i * iFrame[j] );
}
}
}
//TODO: check interface bilinear form
// === coupling entries ===
//evaluate the coupling term
......@@ -373,9 +319,239 @@ private:
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] ]
......@@ -399,39 +575,16 @@ private:
tau[i] += normal[j];
tau[j] = -normal[i];
}
if constexpr (dim != 2)
{
tau /= tau.two_norm();
}
iBasis[i] = tau;
}
return iBasis;
/* //NOTE: which version is better?
if constexpr (dim == 2)
{
Coordinate tau;
tau[0] = normal[1];
tau[1] = -normal[0];
interfaceFrame = {tau};
}
else if constexpr (dim == 3)
{
Coordinate tau1;
tau1[0] = normal[1] + normal[2];
tau1[1] = -normal[0];
tau1[2] = -normal[0];
Coordinate tau2;
tau2[0] = -normal[1];
tau2[1] = normal[0] + normal[2];
tau2[2] = -normal[1];
interfaceFrame = {tau1, tau2};
}
else
{
DUNE_THROW(Dune::Exception,
"Invalid dimension: dim = " + std::to_string(dim));
}*/
}
};
......
dune/mmdg/problems/coupleddgproblem.hh
+ 1
− 1
View file @ 53ab1652
......@@ -61,7 +61,7 @@ class CoupledDGProblem : public DGProblem<Vector, Scalar>
}
//returns the recommended quadrature order to compute an integral
//over d(x) * Kparallel(x)
//over Kparallel(x)
virtual int quadratureOrder_Kparallel () const
{
return 1;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment