diff --git a/dune/mmdg/mmdg.hh b/dune/mmdg/mmdg.hh
index c2040e3d811ba8a29363b9fcb3628431e6aa6754..32c59532a04664f0a9b057bca265792270b46730 100644
--- a/dune/mmdg/mmdg.hh
+++ b/dune/mmdg/mmdg.hh
@@ -186,7 +186,6 @@ private:
// -int_iElem beta * phi_iElem,j * avg(phi_elem,0) ds
//for i = 0,...,dim and j = 0,...,dim-1
//and for elem in {elemIn, elemOut} using a quadrature rule
- //TODO symmetrize more efficiently
const auto lambda_Kperp =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
@@ -292,7 +291,6 @@ private:
// -int_iElem beta * phi_iElem,i * avg(phi_elem1,k) ds
//for i,j = 1,...,dim-1 and k,l = 1,...,dim
//and for elem1,elem2 in {elemIn, elemOut} using a quadrature rule
- //TODO symmetrize more efficiently
const auto lambda_KperpPlus1 =
[&](const Coordinate& qpGlobal, const Scalar weight)
{
diff --git a/dune/mmdg/problems/coupleddgproblem.hh b/dune/mmdg/problems/coupleddgproblem.hh
index ed9ad18be7f05b3c29ccf7ee1486869bc1ce7098..9f7774ff6e3611f92dfcf8b16d21727815fda430 100644
--- a/dune/mmdg/problems/coupleddgproblem.hh
+++ b/dune/mmdg/problems/coupleddgproblem.hh
@@ -57,7 +57,8 @@ public:
return permeability;
}
- //permeability per aperture of the fracture in normal direction at position pos
+ //permeability per aperture of the fracture in normal direction
+ //at position pos
virtual Scalar Kperp (const Coordinate& pos) const
{
return Scalar(1.0);
@@ -71,21 +72,26 @@ public:
}
//returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
+ //over x * d^2 * interfaceBoundary(x)
+ //and x * d * interfaceBoundary(x) * (Kpar * grad(d))
+ //and d^2 * interfaceBoundary(x) * Kpar
virtual int quadratureOrderInterfaceBoundary () const
{
return 1;
}
- //returns the recommended quadrature order to compute an integral
- //over Kparallel(x)
+ //returns the recommended quadrature order to compute integrals
+ //over x^2 * (Kpar * grad(d) * grad(d))
+ //and x^2 * d^2
+ //and x^2 * d * (Kpar * grad(d))
+ //and x * d^2 * Kpar
virtual int quadratureOrder_Kparallel () const
{
- return 1;
+ return 2;
}
//returns the recommended quadrature order to compute an integral
- //over x * Kperp(x) / d(x)
+ //over x * Kperp(x)
virtual int quadratureOrder_Kperp () const
{
return 1;
diff --git a/dune/mmdg/problems/poisson.hh b/dune/mmdg/problems/dgproblem1.hh
similarity index 92%
rename from dune/mmdg/problems/poisson.hh
rename to dune/mmdg/problems/dgproblem1.hh
index 135dc58591255f517e196005698246866a69b32a..39f67481129a76e9217813072d0bf0cf96c54d10 100644
--- a/dune/mmdg/problems/poisson.hh
+++ b/dune/mmdg/problems/dgproblem1.hh
@@ -1,11 +1,11 @@
-#ifndef __DUNE_MMDG_POISSON_HH__
-#define __DUNE_MMDG_POISSON_HH__
+#ifndef __DUNE_MMDG_DGPROBLEM1_HH__
+#define __DUNE_MMDG_DGPROBLEM1_HH__
#include <dune/mmdg/problems/dgproblem.hh>
//a Poisson problem
template<class Coordinate, class Scalar = double>
-class Poisson : public DGProblem<Coordinate, Scalar>
+class DGProblem1 : public DGProblem<Coordinate, Scalar>
{
public:
static constexpr int dim = Coordinate::dimension;
diff --git a/dune/mmdg/problems/inhomogeneousbulkproblem.hh b/dune/mmdg/problems/dgproblem2.hh
similarity index 90%
rename from dune/mmdg/problems/inhomogeneousbulkproblem.hh
rename to dune/mmdg/problems/dgproblem2.hh
index 601d0bb799d03688d75797b12c9fbc29c87402b9..5d8240452bc414b30a6f886b73bdc64d0068709e 100644
--- a/dune/mmdg/problems/inhomogeneousbulkproblem.hh
+++ b/dune/mmdg/problems/dgproblem2.hh
@@ -1,10 +1,10 @@
-#ifndef __DUNE_MMDG_INHOMOGENEOUSBULKPROBLEM_HH__
-#define __DUNE_MMDG_INHOMOGENEOUSBULKPROBLEM_HH__
+#ifndef __DUNE_MMDG_DGPROBLEM2_HH__
+#define __DUNE_MMDG_DGPROBLEM2_HH__
#include <dune/mmdg/problems/dgproblem.hh>
template<class Coordinate, class Scalar = double>
-class InhomogeneousBulkProblem : public DGProblem<Coordinate, Scalar>
+class DGProblem2 : public DGProblem<Coordinate, Scalar>
{
public:
static constexpr int dim = Coordinate::dimension;
diff --git a/dune/mmdg/problems/zeropoisson.hh b/dune/mmdg/problems/dgproblem3.hh
similarity index 88%
rename from dune/mmdg/problems/zeropoisson.hh
rename to dune/mmdg/problems/dgproblem3.hh
index 56525ab2777cd925573f79cd30c426a80cb3018b..29421380f614ec5941cd579769ca6e409179f5c8 100644
--- a/dune/mmdg/problems/zeropoisson.hh
+++ b/dune/mmdg/problems/dgproblem3.hh
@@ -1,15 +1,14 @@
-#ifndef __DUNE_MMDG_ZEROPOISSON_HH__
-#define __DUNE_MMDG_ZEROPOISSON_HH__
+#ifndef __DUNE_MMDG_DGPROBLEM3_HH__
+#define __DUNE_MMDG_DGPROBLEM3_HH__
#include <cmath>
#include <string>
#include <dune/mmdg/problems/dgproblem.hh>
-
//Poisson problem with right side q = -1
template<class Coordinate, class Scalar = double>
-class ZeroPoisson : public DGProblem<Coordinate, Scalar>
+class DGProblem3 : public DGProblem<Coordinate, Scalar>
{
public:
static constexpr int dim = Coordinate::dimension;
@@ -19,7 +18,7 @@ class ZeroPoisson : public DGProblem<Coordinate, Scalar>
{
if constexpr (dim == 1)
{
- return 0.5*pos[0]*(pos[0] - 1);
+ return 0.5 * pos[0] * (pos[0] - 1);
}
if constexpr (dim == 2) //approximate series
diff --git a/dune/mmdg/problems/mmdgproblem1.hh b/dune/mmdg/problems/mmdgproblem1.hh
index 6ea0cd89881315afa5d63d895939cc4ff5b7ba17..f35a61ff999ac3eb073e87408dc67cf92c07a1b1 100644
--- a/dune/mmdg/problems/mmdgproblem1.hh
+++ b/dune/mmdg/problems/mmdgproblem1.hh
@@ -47,13 +47,6 @@ public:
return 2;
}
- //returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
- int quadratureOrderInterfaceBoundary () const
- {
- return 2;
- }
-
private:
Scalar d_; //aperture
};
diff --git a/dune/mmdg/problems/mmdgproblem3.hh b/dune/mmdg/problems/mmdgproblem3.hh
index a6914da6ffd9d6038bfbc3751818f82bae36e14f..d3f284c1c2a69d257bdf3acd9a11583477eb876d 100644
--- a/dune/mmdg/problems/mmdgproblem3.hh
+++ b/dune/mmdg/problems/mmdgproblem3.hh
@@ -79,20 +79,6 @@ public:
return 10;
}
- //returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
- int quadratureOrderInterfaceBoundary () const
- {
- return 10;
- }
-
- //returns the recommended quadrature order to compute an integral
- //over x * qInterface(x)
- int quadratureOrder_qInterface () const
- {
- return 10;
- }
-
private:
Scalar d_; //aperture
};
diff --git a/dune/mmdg/problems/mmdgproblem4.hh b/dune/mmdg/problems/mmdgproblem4.hh
index 1cae12485e35479a48b48575b82bb13a786abf9e..cf09ef6b727e7b946854d364c29605e582e9b6ba 100644
--- a/dune/mmdg/problems/mmdgproblem4.hh
+++ b/dune/mmdg/problems/mmdgproblem4.hh
@@ -76,28 +76,16 @@ public:
//over x * boundary(x)
int quadratureOrderBoundary () const
{
- return 10;
+ return 3;
}
//returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
+ //over x * d^2 * interfaceBoundary(x)
+ //and x * d * interfaceBoundary(x) * (Kpar * grad(d))
+ //and d^2 * interfaceBoundary(x) * Kpar
int quadratureOrderInterfaceBoundary () const
{
- return 10;
- }
-
- //returns the recommended quadrature order to compute an integral
- //over x * qInterface(x)
- int quadratureOrder_qInterface () const
- {
- return 10;
- }
-
- //returns the recommended quadrature order to compute an integral
- //over x * q(x)
- int quadratureOrder_q () const
- {
- return 10;
+ return 3;
}
private:
diff --git a/dune/mmdg/problems/mmdgproblem5.hh b/dune/mmdg/problems/mmdgproblem5.hh
index 34ade2f4875a6579c6ac17878c542c8c87384712..ee614898bdfdaadf48f5e3b1e97aa5986101799d 100644
--- a/dune/mmdg/problems/mmdgproblem5.hh
+++ b/dune/mmdg/problems/mmdgproblem5.hh
@@ -100,7 +100,9 @@ public:
}
//returns the recommended quadrature order to compute an integral
- //over x * interfaceBoundary(x)
+ //over x * d^2 * interfaceBoundary(x)
+ //and x * d * interfaceBoundary(x) * (Kpar * grad(d))
+ //and d^2 * interfaceBoundary(x) * Kpar
int quadratureOrderInterfaceBoundary () const
{
return 10;
@@ -120,6 +122,23 @@ public:
return 10;
}
+ //returns the recommended quadrature order to compute integrals
+ //over x^2 * (Kpar * grad(d) * grad(d))
+ //and x^2 * d^2
+ //and x^2 * d * (Kpar * grad(d))
+ //and x * d^2 * Kpar
+ int quadratureOrder_Kparallel () const
+ {
+ return 5;
+ }
+
+ //returns the recommended quadrature order to compute an integral
+ //over an entry K[i][j] of the permeability tensor
+ int quadratureOrder_K () const
+ {
+ return 10;
+ }
+
private:
Scalar d_; //aperture
diff --git a/src/dg.cc b/src/dg.cc
index dcc46594c5300a78bcf52a631e2577a08d81a328..1b2149213438ff800941d6d8554d135a1ad61070 100644
--- a/src/dg.cc
+++ b/src/dg.cc
@@ -21,9 +21,9 @@
#include <dune/mmdg/dg.hh>
#include <dune/mmdg/clockhelper.hh>
#include <dune/mmdg/problems/dgproblem.hh>
-#include <dune/mmdg/problems/zeropoisson.hh>
-#include <dune/mmdg/problems/poisson.hh>
-#include <dune/mmdg/problems/inhomogeneousbulkproblem.hh>
+#include <dune/mmdg/problems/dgproblem1.hh>
+#include <dune/mmdg/problems/dgproblem2.hh>
+#include <dune/mmdg/problems/dgproblem3.hh>
int main(int argc, char** argv)
{
@@ -65,23 +65,29 @@ int main(int argc, char** argv)
Problem* problem; //problem to be solved
//determine problem type
- if (pt["problem"] == "zeropoisson")
+ if (pt["problem"] == "dg1")
{
- problem = new ZeroPoisson<Coordinate>();
+ problem = new DGProblem1<Coordinate>();
}
- else if (pt["problem"] == "poisson")
+ else if (pt["problem"] == "dg2")
{
- problem = new Poisson<Coordinate>();
+ problem = new DGProblem2<Coordinate>();
}
- else if (pt["problem"] == "inhomogeneousbulk")
+ else if (pt["problem"] == "dg3")
{
- problem = new InhomogeneousBulkProblem<Coordinate>();
+ problem = new DGProblem3<Coordinate>();
+
+ if constexpr (dim == 2)
+ {
+ std::cout << "NOTE: The exact solution of your problem is given by "
+ "an infinite series; increase the cut-off if necessary!\n\n";
+ }
}
else
{
DUNE_THROW(Dune::Exception,
"Invalid problem type or parameter 'problem' not specified in file "
- "parameterDG.ini (use 'poisson' or TODO)");
+ "parameterDG.ini (use 'dg1' to 'dg3')");
}
DG dg(gridView, mapper, *problem);
diff --git a/src/parameterDG.ini b/src/parameterDG.ini
index 10896d307fc1d477c906a7d92644cee815001212..e3e1b98f436304cfd37913f8d0d2d0a8e7de9c36 100644
--- a/src/parameterDG.ini
+++ b/src/parameterDG.ini
@@ -1,2 +1,2 @@
mu0 = 1000
-problem = poisson
+problem = dg1