From 45af3c7bb87bbe751de27c1b0d060b79989bccd0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
<maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Tue, 9 Jun 2020 23:50:40 +0200
Subject: [PATCH] [WIP] add problem dg5
---
dune/mmdg/problems/dgproblem5.hh | 139 +++++++++++++++++++++++++++++++
1 file changed, 139 insertions(+)
create mode 100644 dune/mmdg/problems/dgproblem5.hh
diff --git a/dune/mmdg/problems/dgproblem5.hh b/dune/mmdg/problems/dgproblem5.hh
new file mode 100644
index 0000000..1a6e4d5
--- /dev/null
+++ b/dune/mmdg/problems/dgproblem5.hh
@@ -0,0 +1,139 @@
+#ifndef __DUNE_MMDG_DGPROBLEM5_HH__
+#define __DUNE_MMDG_DGPROBLEM5_HH__
+
+#include <dune/mmdg/problems/dgproblem.hh>
+#include <cmath>
+
+template<class Coordinate, class Scalar = double>
+class DGProblem5 : public DGProblem<Coordinate, Scalar>
+{
+ public:
+ static constexpr int dim = Coordinate::dimension;
+ using Base = DGProblem<Coordinate, Scalar>;
+ using Matrix = typename Base::Matrix;
+
+ static constexpr Scalar dmin = 0.05; //minimum aperture of the fracture
+ static constexpr Scalar dmax = 0.1; //maximum aperture of the fracture
+ static constexpr Scalar Kbulk = 1.0; //bulk permeability
+ static constexpr Scalar K_par = 3.0; //tangential fracture permeability
+ static constexpr Scalar K_perp = 0.5; //normal fracture permeability
+
+ //require dim == 2
+ DGProblem5 ()
+ {
+ if (dim != 2)
+ {
+ DUNE_THROW(Dune::Exception,
+ "Problem dg4 is not implemented for dimension = " +
+ std::to_string(dim) + "!");
+ }
+ }
+
+ //the exact solution at position pos
+ Scalar exactSolution (const Coordinate& pos) const
+ {
+ if ( pos[0] < 0.5 * (1.0 - aperture(pos)) )
+ { //pos in Omega_1 (left of the fracture)
+ return -0.5 * ( pos[0] - 0.5 ) * pos[1] * aperture(pos);
+ }
+
+ if ( pos[0] > 0.5 * (1.0 + aperture(pos)) )
+ { //pos in Omega_2 (right of the fracture)
+ return 0.5 * ( pos[0] - 0.5 ) * pos[1] * aperture(pos);
+ }
+
+ //pos in Omega_R (inside the fracture)
+ return ( pos[0] - 0.5 ) * ( pos[0] - 0.5 ) * pos[1];
+ }
+
+ //exact solution is implemented
+ bool hasExactSolution() const
+ {
+ return true;
+ };
+
+ //source term at position pos
+ Scalar q (const Coordinate& pos) const
+ {
+ if ( pos[0] < 0.5 * (1.0 - aperture(pos)) )
+ { //pos in Omega_1 (left of the fracture)
+ return -4.0 * M_PI * Kbulk * ( dmax - dmin ) *
+ ( pos[0] - 0.5 ) * ( std::sin(8.0 * M_PI * pos[1])
+ + 4.0 * M_PI * pos[1] * std::cos(8.0 * M_PI * pos[1]) );
+ }
+
+ if ( pos[0] > 0.5 * (1.0 + aperture(pos)) )
+ { //pos in Omega_2 (right of the fracture)
+ return 4.0 * M_PI * Kbulk * ( dmax - dmin ) *
+ ( pos[0] - 0.5 ) * ( std::sin(8.0 * M_PI * pos[1])
+ + 4.0 * M_PI * pos[1] * std::cos(8.0 * M_PI * pos[1]) );
+ }
+
+ //pos in Omega_R (inside the fracture)
+ return -2.0 * K_perp * pos[1];
+ // return -2.0 * K_perp * pos[1] - ( pos[0] - 0.5 ) * ( pos[0] - 0.5 ) *
+ // ( -1.0 + 4.0 * M_PI * ( dmax - dmin ) * ( 8.0 * M_PI *
+ // std::cos(8.0 * M_PI * pos[1]) - std::sin(8.0 * M_PI * pos[1]) )
+ // / aperture(pos) );
+ }
+
+ //permeability tensor at position pos
+ Matrix K (const Coordinate& pos) const
+ {
+ Matrix permeability(0.0);
+
+ if ( pos[0] < 0.5 * (1.0 - aperture(pos)) ||
+ pos[0] > 0.5 * (1.0 + aperture(pos)))
+ { //pos in Omega_1 or Omega_2 (outside the fracture)
+ for (int i = 0; i < dim; i++)
+ {
+ permeability[i][i] = Kbulk;
+ // aperture(pos) / ( aperture(pos) - 4.0 * pos[1] * M_PI * (dmax - dmin) * std::sin(8.0 * M_PI * pos[1]) );
+
+ if (permeability[i][i] <= 0.0)
+ {
+ std::cout << "WARNING!" << std::endl;
+ }
+ }
+ }
+ else
+ { //pos in Omega_R (inside the fracture)
+ permeability[0][0] = K_perp;
+
+ for (int i = 1; i < dim; i++)
+ {
+ permeability[i][i] = K_par;
+ // 1.0 - pos[1] * 4.0 * M_PI * ( dmax - dmin ) *
+ // std::sin(8.0 * M_PI * pos[1]) / aperture(pos);
+
+ if (permeability[i][i] <= 0.0)
+ {
+ std::cout << "WARNING" << std::endl;
+ }
+ }
+ }
+
+ return permeability;
+ }
+
+ private:
+ //aperture of the fracture
+ const Scalar aperture (const Coordinate& pos) const
+ {
+ return dmin + 0.5 * ( dmax - dmin ) *
+ ( 1.0 + std::cos(8.0 * M_PI * pos[1]) );
+ //return dmin + ( dmax - dmin ) * std::cos(8.0 * M_PI * pos[1]);
+
+ // const Scalar y = pos[1] - 0.25 * floor(4.0 * pos[1]);
+ //
+ // if ( y <= 0.125 )
+ // { //0 <= y <= 0.125
+ // return dmin + ( 1.0 - 8.0 * y ) * ( dmax - dmin );
+ // }
+ //
+ // //0.125 <= y <= 0.25
+ // return dmin + ( 8.0 * y - 1.0 ) * ( dmax - dmin );
+ }
+};
+
+#endif
--
GitLab