From ec72348baa798ae60462ac4c6b1dbde00d854d40 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 21:26:52 +0200
Subject: [PATCH] add terms with grad(d) to mmdg
---
dune/mmdg/mmdg.hh | 173 ++++++++++++++++---------
dune/mmdg/problems/coupleddgproblem.hh | 6 +
dune/mmdg/problems/mmdgproblem2.hh | 6 +-
dune/mmdg/problems/mmdgproblem4.hh | 10 +-
4 files changed, 125 insertions(+), 70 deletions(-)
diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index 6725182..cb59505 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -437,33 +437,66 @@ private:
//lambda to evaluate
// int_iElem (Kparallel * grad(d phi_iElem,i)) * grad(d phi_iElem,j) ds
- // = int_iElem (Kpar * tau_i) * tau_j ds //TODO: terms with grad(d)
- //using a quadrature rule for i,j = 1,...,dim-1
+ //using a quadrature rule for i,j = 0,...,dim-1
const auto lambda_Kpar =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
- CoordinateMatrix KEvaluation = Base::problem_.Kparallel(qpGlobal);
- Scalar dEvaluation2 = Base::problem_.aperture(qpGlobal);
- dEvaluation2 *= dEvaluation2;
+ const CoordinateMatrix KEvaluation =
+ Base::problem_.Kparallel(qpGlobal);
+ const Scalar dEvaluation = Base::problem_.aperture(qpGlobal);
+ const Scalar dEvaluation2 = dEvaluation * dEvaluation;
+ const Coordinate gradDEvaluation =
+ Base::problem_.gradAperture(qpGlobal);
Coordinate KparDotTau_i; //storage for Kparallel * iFrame[i]
+ Coordinate KparDotGradD; //storage for Kparallel * grad(d)
- for (int i = 0; i < dim-1 ; i++)
+ KEvaluation.mv(gradDEvaluation, KparDotGradD);
+
+ const Scalar interfaceUpdate1 = KparDotGradD * gradDEvaluation;
+
+ //quadrature for
+ // int_iElem (Kpar * grad(d * phi_iElem,0)) * grad(d * phi,iElem,0) ds
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE] += weight * interfaceUpdate1;
+
+ for (int i = 0; i < dim - 1 ; i++)
{
KEvaluation.mv(iFrame[i], KparDotTau_i);
+ const Scalar qpI = qpGlobal * iFrame[i];
+ const Scalar interfaceUpdate2 =
+ dEvaluation * (KparDotGradD * iFrame[i]);
+ const Scalar interfaceUpdate3 =
+ weight * (qpI * interfaceUpdate1 + interfaceUpdate2);
+
+ //quadrature for
+ // int_iElem (Kpar * grad(d * phi_iElem,i))
+ // * grad(d * phi,iElem,0) ds
+ (*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] += interfaceUpdate3;
+ (*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] += interfaceUpdate3;
+
+ //quadrature for
+ // int_iElem (Kpar * grad(d * phi_iElem,i))
+ // * grad(d * phi,iElem,i) ds
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
- weight * dEvaluation2 * ( KparDotTau_i * iFrame[i] );
+ weight * ( dEvaluation2 * (KparDotTau_i * iFrame[i])
+ + qpI * (qpI * interfaceUpdate1 + 2.0 * interfaceUpdate2) );
for (int j = 0; j < i; j++)
{
- const Scalar interfaceUpdate =
- weight * dEvaluation2 * ( KparDotTau_i * iFrame[j] );
+ const Scalar interfaceUpdate4 =
+ weight * ( dEvaluation2 * (KparDotTau_i * iFrame[j])
+ + dEvaluation * qpI * (KparDotGradD * iFrame[j])
+ + (qpGlobal * iFrame[j])
+ * (interfaceUpdate2 + qpI * interfaceUpdate1) );
+ //quadrature for
+ // int_iElem (Kpar * grad(d * phi_iElem,i))
+ // * grad(d * phi,iElem,j) ds
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
- interfaceUpdate;
+ interfaceUpdate4;
(*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
- interfaceUpdate;
+ interfaceUpdate4;
}
}
};
@@ -494,7 +527,7 @@ private:
// -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)
+ //using a quadrature rule
const auto lambda_Kpar_iFacet =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
@@ -503,62 +536,70 @@ private:
const CoordinateMatrix KEvaluation =
Base::problem_.Kparallel(qpGlobal);
- Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i]
- Coordinate KparDotTau_j; //storage for K_parallel * iFrame[j]
+ Coordinate KparDotTau_i; //storage for Kparallel * iFrame[i]
+ Coordinate KparDotTau_j; //storage for Kparallel * iFrame[j]
+ Coordinate KparDotGradD; //storage for Kparallel * grad(d)
+
+ KEvaluation.mv(Base::problem_.gradAperture(qpGlobal),
+ KparDotGradD);
+
+ const Scalar interfaceUpdate1 = KparDotGradD * normal;
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
(*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
- weight * mu * dEvaluation2;
+ weight * dEvaluation * (mu * dEvaluation - interfaceUpdate1);
for (int i = 0; i < dim - 1; i++)
{
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 = weight * dEvaluation2
- * (interfaceUpdate1 - 0.5 * interfaceUpdate2);
+ const Scalar interfaceUpdate2 = mu * dEvaluation * qpI;
+ const Scalar interfaceUpdate3 =
+ dEvaluation * (KparDotTau_i * normal);
+ const Scalar interfaceUpdate4 = weight * dEvaluation
+ * (interfaceUpdate2 - 0.5 * interfaceUpdate3
+ - qpI * interfaceUpdate1);
//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)
+ // * jump(phi_iElem,0) dr
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
- interfaceUpdate3;
+ interfaceUpdate4;
(*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
- interfaceUpdate3;
+ interfaceUpdate4;
//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)
+ // * jump(phi_iElem,i) dr
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + i + 1] +=
- weight * dEvaluation2 * qpI *
- (interfaceUpdate1 - interfaceUpdate2);
+ weight * dEvaluation * qpI * (interfaceUpdate2
+ - interfaceUpdate3 - qpI * interfaceUpdate1);
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));
+ const Scalar interfaceUpdate5 = weight * dEvaluation *
+ ( (qpGlobal * iFrame[j]) * ( interfaceUpdate2
+ - 0.5 * interfaceUpdate3 - qpI * interfaceUpdate1 )
+ - (KparDotTau_j * normal) * qpI * dEvaluation );
//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)
+ // * jump(phi_iElem,j) dr
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
- interfaceUpdate4;
+ interfaceUpdate5;
(*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
- interfaceUpdate4;
+ interfaceUpdate5;
}
}
@@ -605,7 +646,7 @@ private:
// * 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)
+ // * jump(phi_neighbor,0) dr
(*Base::A)[iElemIdxSLE + i + 1][neighborIdxSLE] +=
interfaceUpdate3a;
(*Base::A)[neighborIdxSLE][iElemIdxSLE + i + 1] +=
@@ -616,7 +657,7 @@ private:
// * 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)
+ // * jump(phi_iElem,0) dr
(*Base::A)[neighborIdxSLE + i + 1][iElemIdxSLE] +=
interfaceUpdate3b;
(*Base::A)[iElemIdxSLE][neighborIdxSLE + i + 1] +=
@@ -642,7 +683,6 @@ private:
//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] +=
@@ -662,33 +702,39 @@ private:
// -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);
Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i]
Coordinate KparDotTau_j; //storage for K_parallel * iFrame[j]
+ Coordinate KparDotGradD; //storage for K_parallel * grad(d)
+
+ KEvaluation.mv(Base::problem_.gradAperture(qpGlobal),
+ KparDotGradD);
+
+ const Scalar interfaceUpdate1 = 2.0 * (KparDotGradD * normal);
//quadrature for
// int_intersct mu * d^2 * jump(phi_iElem,0)
// * jump(phi_iElem,0) dr
(*Base::A)[iElemIdxSLE][iElemIdxSLE] +=
- weight * mu * dEvaluation2;
+ weight * dEvaluation * (mu * dEvaluation - interfaceUpdate1);
for (int i = 0; i < dim - 1; i++)
{
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 interfaceUpdate2 = mu * dEvaluation * qpI;
const Scalar interfaceUpdate3 =
- weight * dEvaluation2 * (interfaceUpdate1 - interfaceUpdate2);
+ dEvaluation * (KparDotTau_i * normal);
+ const Scalar interfaceUpdate4 = weight * dEvaluation *
+ (interfaceUpdate2 - interfaceUpdate3
+ - qpI * interfaceUpdate1);
//quadrature for
// int_intersct mu * d^2 * phi_iElem,0 * phi_iElem,i dr
@@ -696,9 +742,9 @@ private:
// -int_intersct d * (Kparallel * grad(d * phi_iElem,i))
// * phi_iElem,0 * normal dr
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE] +=
- interfaceUpdate3;
+ interfaceUpdate4;
(*Base::A)[iElemIdxSLE][iElemIdxSLE + i + 1] +=
- interfaceUpdate3;
+ interfaceUpdate4;
//quadrature for
// int_intersct mu * d^2 * (phi_iElem,i)^2 dr
@@ -706,17 +752,19 @@ private:
// -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);
+ weight * qpI * dEvaluation *
+ (interfaceUpdate2 - 2.0 * interfaceUpdate3
+ - qpI * qpI * interfaceUpdate1);
for (int j = 0; j < i; j++)
{
KEvaluation.mv(iFrame[j], KparDotTau_j);
- const Scalar interfaceUpdate4 = weight * dEvaluation2 *
- ((qpGlobal * iFrame[j]) * interfaceUpdate1
- - (interfaceUpdate2 * (qpGlobal * iFrame[j])
- + (KparDotTau_j * normal) * qpI));
+ const Scalar qpJ = qpGlobal * iFrame[j];
+ const Scalar interfaceUpdate5 =
+ weight * dEvaluation * (qpJ * (interfaceUpdate2
+ - interfaceUpdate3 - qpI * interfaceUpdate1)
+ - (KparDotTau_j * normal) * qpI * dEvaluation);
//quadrature for
// int_intersct mu * d^2 * phi_iElem,i * phi_iElem,j dr
@@ -724,9 +772,9 @@ private:
// -int_intersct d * (Kparallel * grad(d * phi_iElem,i))
// * phi_iElem,j * normal dr
(*Base::A)[iElemIdxSLE + i + 1][iElemIdxSLE + j + 1] +=
- interfaceUpdate4;
+ interfaceUpdate5;
(*Base::A)[iElemIdxSLE + j + 1][iElemIdxSLE + i + 1] +=
- interfaceUpdate4;
+ interfaceUpdate5;
}
}
};
@@ -736,28 +784,33 @@ private:
//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;
+ Coordinate KparDotTau_i; //storage for K_parallel * iFrame[i]
+ Coordinate KparDotGradD; //storage for K_parallel * grad(d)
+
+ KEvaluation.mv(Base::problem_.gradAperture(qpGlobal),
+ KparDotGradD);
+
+ const Scalar interfaceUpdate1 = weight * dEvaluation *
+ Base::problem_.interfaceBoundary(qpGlobal);
+ const Scalar interfaceUpdate2 =
+ mu * dEvaluation - KparDotGradD * normal;
- Base::b[iElemIdxSLE] += mu * d2gGamma;
+ Base::b[iElemIdxSLE] += interfaceUpdate1 * interfaceUpdate2;
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);
+ Base::b[iElemIdxSLE + i + 1] +=
+ interfaceUpdate1 * ((qpGlobal * iFrame[i]) * interfaceUpdate2
+ - dEvaluation * (KparDotTau_i * normal));
}
};
diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh
index 10488bf..4cc93b5 100644
--- a/dune/mmdg/problems/coupleddgproblem.hh
+++ b/dune/mmdg/problems/coupleddgproblem.hh
@@ -32,6 +32,12 @@ public:
return Scalar(1.0);
}
+ //tangential gradient of the aperture d of the fracture at position pos
+ virtual Coordinate gradAperture (const Coordinate& pos) const
+ {
+ return Coordinate(0.0);
+ }
+
virtual Scalar interfaceBoundary(const Coordinate& pos) const
{
return exactInterfaceSolution(pos);
diff --git a/dune/mmdg/problems/mmdgproblem2.hh b/dune/mmdg/problems/mmdgproblem2.hh
index 624cefc..7b074e7 100644
--- a/dune/mmdg/problems/mmdgproblem2.hh
+++ b/dune/mmdg/problems/mmdgproblem2.hh
@@ -17,7 +17,8 @@ public:
//the exact bulk solution at position pos
Scalar exactSolution (const Coordinate& pos) const
{
- Scalar solution = (pos[0] < 0.5) ? std::sin(4 * pos[0]) : std::cos(4 * pos[0]);
+ Scalar solution =
+ (pos[0] < 0.5) ? std::sin(4 * pos[0]) : std::cos(4 * pos[0]);
return solution * std::cos(M_PI * pos[1]);
}
@@ -38,7 +39,8 @@ public:
Scalar q (const Coordinate& pos) const
{
constexpr Scalar sourcefactor = 16.0 + M_PI * M_PI;
- Scalar source = (pos[0] < 0.5) ? std::sin(4 * pos[0]) : std::cos(4 * pos[0]);
+ Scalar source =
+ (pos[0] < 0.5) ? std::sin(4 * pos[0]) : std::cos(4 * pos[0]);
return sourcefactor * source * std::cos(M_PI * pos[1]);
}
diff --git a/dune/mmdg/problems/mmdgproblem4.hh b/dune/mmdg/problems/mmdgproblem4.hh
index 61d19dc..1cae124 100644
--- a/dune/mmdg/problems/mmdgproblem4.hh
+++ b/dune/mmdg/problems/mmdgproblem4.hh
@@ -32,7 +32,7 @@ public:
//the exact solution on the interface at position pos
Scalar exactInterfaceSolution (const Coordinate& pos) const
{
- return pos[0] * pos[0] - pos[1] * pos[1] + 0.5 * d_;
+ return 0.25 - pos[1] * pos[1] + 0.5 * d_;
}
//indicates whether an exact solution is implemented for the problem
@@ -41,12 +41,6 @@ public:
return true;
};
- //source term at position pos
- Scalar q (const Coordinate& pos) const
- {
- return 0.0;
- }
-
//interface source term at position pos
Scalar qInterface (const Coordinate& pos) const
{
@@ -88,7 +82,7 @@ public:
//returns the recommended quadrature order to compute an integral
//over x * interfaceBoundary(x)
int quadratureOrderInterfaceBoundary () const
- { //NOTE: only relevant for 3d, not implemented!
+ {
return 10;
}
--
GitLab