From e62bdc7488b650104bb0b8903c9ab5d6f6d0e28e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Maximilian=20H=C3=B6rl?=
 <maximilian.hoerl@mathematik.uni-stuttgart.de>
Date: Sat, 11 Apr 2020 16:14:08 +0200
Subject: [PATCH] change problem mmdg6 to allow arbitrary xi

---
 dune/mmdg/problems/mmdgproblem6.hh | 18 +++++++++++++++---
 src/mmdg.cc                        |  4 ++--
 2 files changed, 17 insertions(+), 5 deletions(-)

diff --git a/dune/mmdg/problems/mmdgproblem6.hh b/dune/mmdg/problems/mmdgproblem6.hh
index 36966a5..c191e5a 100644
--- a/dune/mmdg/problems/mmdgproblem6.hh
+++ b/dune/mmdg/problems/mmdgproblem6.hh
@@ -12,7 +12,10 @@ public:
   static constexpr int dim = Coordinate::dimension;
 
   //constructor
-  MMDGProblem6 (const Scalar d0) : d0_(d0) {}
+  MMDGProblem6 (const Scalar d0) : d0_(d0), xi_(1.0) {}
+
+  //constructor
+  MMDGProblem6 (const Scalar d0, const Scalar xi) : d0_(d0), xi_(xi) {}
 
   //the exact bulk solution at position pos
   Scalar exactSolution (const Coordinate& pos) const
@@ -26,7 +29,11 @@ public:
   //the exact solution on the interface at position pos
   Scalar exactInterfaceSolution (const Coordinate& pos) const
   {
-    return std::cos(4.0 * pos[1]);
+    constexpr Scalar coshsqrt2 = std::cosh(sqrt(2.0));
+    constexpr Scalar sinhsqrt2 = std::sinh(sqrt(2.0));
+
+    return std::cos(4.0 * pos[1]) * (coshsqrt2 * coshsqrt2
+      - (2 * xi_ - 1) * sinhsqrt2 * sinhsqrt2);
   }
 
   //indicates whether an exact solution is implemented for the problem
@@ -49,8 +56,12 @@ public:
   {
     const Scalar cos4y = std::cos(4.0 * pos[1]);
     constexpr Scalar sqrt2 = 2.0 * sqrt(2.0);
+    constexpr Scalar coshsqrt2 = std::cosh(sqrt(2.0));
+    constexpr Scalar sinhsqrt2 = std::sinh(sqrt(2.0));
+    const Scalar hyperbolicFactor =
+      (coshsqrt2 * coshsqrt2 - (2 * xi_ - 1) * sinhsqrt2 * sinhsqrt2);
 
-    return -48.0 * d0_ * d0_ * std::exp(-8.0 * pos[1]) * ( cos4y
+    return -48.0 * d0_ * d0_ * std::exp(-8.0 * pos[1]) * hyperbolicFactor * ( cos4y
       + 3.0 * std::sin(4.0 * pos[1]) ) - sqrt2 * std::sinh(sqrt2) * cos4y;
   }
 
@@ -117,6 +128,7 @@ public:
 
 private:
   Scalar d0_; //prefactor of the aperture
+  Scalar xi_; //coupling constant
 
 };
 
diff --git a/src/mmdg.cc b/src/mmdg.cc
index 706813e..11a081b 100644
--- a/src/mmdg.cc
+++ b/src/mmdg.cc
@@ -95,9 +95,9 @@ int main(int argc, char** argv)
     }
     else if (pt["problem"] == "mmdg6")
     {
-      problem = new MMDGProblem6<Coordinate>(aperture);
+      xi = (pt.hasKey("xi")) ? std::stod(pt["xi"]) : 1.0;
+      problem = new MMDGProblem6<Coordinate>(aperture, xi);
       gridType = "mmdg2_5e-2";
-      xi = 1.0;
     }
     else
     {
-- 
GitLab