moved Problems to cpp

mluq/src/problems/StochasticTransportProblem.cpp

 #include "StochasticTransportProblem.hpp"
 
 
+class StochasticPollution1D : public IStochasticTransportProblem {
+  HybridFaceNormalFluxGenerator scalarGenerator;
+  HybridCellFluxGenerator vectorFieldGenerator;
+public:
+  explicit StochasticPollution1D(const Meshes &meshes) :
+    IStochasticTransportProblem(meshes),
+    scalarGenerator(HybridFaceNormalFluxGenerator(meshes)),
+    vectorFieldGenerator(HybridCellFluxGenerator(meshes, scalarGenerator)) {}
+
+  void DrawSample(const SampleID &id) override {
+    scalarGenerator.DrawSample(id);
+    vectorFieldGenerator.DrawSample(id);
+  }
+
+  Scalar Solution(double t, const Point &x) const override {
+    if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
+    return 1.0;
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return this->vectorFieldGenerator.EvalSample(c);
+  }
+
+  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
+                        const Point &x) const override {
+    return this->scalarGenerator.EvalSample(face, c);
+  }
+
+  string Name() const override { return "StochasticPollution1D"; }
+};
+
+class Pollution1D : public IStochasticTransportProblem {
+public:
+  explicit Pollution1D(const Meshes &meshes) : IStochasticTransportProblem(meshes) {}
+
+  void DrawSample(const SampleID &id) override {}
+
+  Scalar Solution(double t, const Point &x) const override {
+    if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
+    return 1.0;
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return VectorField(-1.0, 0.0);
+  }
+
+  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
+                        const Point &x) const override {
+    return 1.0;
+  }
+
+  string Name() const override { return "Pollution1D"; }
+};
+
+class StochasticPollutionCosHat1D : public IStochasticTransportProblem {
+  HybridFaceNormalFluxGenerator scalarGenerator;
+
+  HybridCellFluxGenerator vectorFieldGenerator;
+
+  double amplitude = 1.00;
+
+  double cc = 6.0;
+
+public:
+  explicit StochasticPollutionCosHat1D(const Meshes &meshes) :
+    IStochasticTransportProblem(meshes),
+    scalarGenerator(HybridFaceNormalFluxGenerator(meshes)),
+    vectorFieldGenerator(HybridCellFluxGenerator(meshes, scalarGenerator)) {}
+
+  void DrawSample(const SampleID &id) override {
+    scalarGenerator.DrawSample(id);
+    vectorFieldGenerator.DrawSample(id);
+  }
+
+  Scalar Solution(double t, const Point &x) const override {
+    Point midPoint = Point(0.2, 0.0);
+    double r = dist(midPoint, x);
+    if (r < 1 / cc)
+      return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
+    return 0.0;
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return this->vectorFieldGenerator.EvalSample(c);
+  }
+
+  Scalar FaceNormalFlux(const cell &c, int face,
+                        const VectorField &N, const Point &x) const override {
+    return this->scalarGenerator.EvalSample(face, c);
+  }
+
+  string Name() const override { return "StochasticPollutionCosHat1D"; }
+};
+
+class StochasticPollution2D : public IStochasticTransportProblem {
+  HybridFaceNormalFluxGenerator scalarGenerator;
+
+  HybridCellFluxGenerator vectorFieldGenerator;
+
+  double d = 1.0 / 16.0;
+
+public:
+  explicit StochasticPollution2D(const Meshes &meshes) :
+    IStochasticTransportProblem(meshes),
+    scalarGenerator(HybridFaceNormalFluxGenerator(meshes)),
+    vectorFieldGenerator(HybridCellFluxGenerator(meshes, scalarGenerator)) {}
+
+  void DrawSample(const SampleID &id) override {
+    scalarGenerator.DrawSample(id);
+    vectorFieldGenerator.DrawSample(id);
+  }
+
+  Scalar Solution(double t, const Point &x) const override {
+    if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 3 * d))
+      return 1.0;
+    return 0.0;
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return this->vectorFieldGenerator.EvalSample(c);
+  }
+
+  Scalar FaceNormalFlux(const cell &c, int face,
+                        const VectorField &N, const Point &x) const override {
+    return this->scalarGenerator.EvalSample(face, c);
+  }
+
+  string Name() const override { return "StochasticPollution2D"; }
+};
+
+class Pollution2D : public IStochasticTransportProblem {
+private:
+  double d = 1.0 / 16.0;
+
+public:
+  explicit Pollution2D(const Meshes &meshes) : IStochasticTransportProblem(meshes) {}
+
+  void DrawSample(const SampleID &id) override {}
+
+  Scalar Solution(double t, const Point &x) const override {
+    if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 4 * d))
+      return 32.0;
+    return 0.0;
+  }
+
+  static VectorField TransportFlux(const Point &x) {
+    return VectorField(0.0, -1.0);
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return TransportFlux(x);
+  }
+
+  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
+                        const Point &x) const override {
+    return TransportFlux(x) * N;
+  }
+
+  string Name() const override { return "Pollution2D"; }
+};
+
+class StochasticPollutionCosHat2D : public IStochasticTransportProblem {
+  HybridFaceNormalFluxGenerator scalarGenerator;
+
+  HybridCellFluxGenerator vectorFieldGenerator;
+
+  double amplitude = 1.00;
+
+  double cc = 6.0;
+
+public:
+  explicit StochasticPollutionCosHat2D(const Meshes &meshes) :
+    IStochasticTransportProblem(meshes),
+    scalarGenerator(HybridFaceNormalFluxGenerator(meshes)),
+    vectorFieldGenerator(HybridCellFluxGenerator(meshes, scalarGenerator)) {}
+
+  void DrawSample(const SampleID &id) override {
+    scalarGenerator.DrawSample(id);
+    vectorFieldGenerator.DrawSample(id);
+  }
+
+  double Solution(double t, const Point &x) const override {
+    Point midPoint = Point(0.5, 0.8);
+    double r = dist(midPoint, x);
+    if (r < 1 / cc)
+      return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
+    return 0.0;
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return this->vectorFieldGenerator.EvalSample(c);
+  }
+
+  Scalar FaceNormalFlux(const cell &c, int face,
+                        const VectorField &N, const Point &x) const override {
+    return this->scalarGenerator.EvalSample(face, c);
+  }
+
+  string Name() const override { return "StochasticPollutionCosHat2D"; }
+};
+
+class PollutionCosHat2D : public IStochasticTransportProblem {
+  double amplitude = 2.0 / Pi;
+public:
+  explicit PollutionCosHat2D(const Meshes &meshes)
+    : IStochasticTransportProblem(meshes) {}
+
+  void DrawSample(const SampleID &id) override {}
+
+  double Solution(double t, const Point &x) const override {
+    Point midPoint = 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t))
+      + Point(1.0, 1.0));
+    double r = dist(midPoint, x);
+    if (r < 1.0 / (2.0 * Pi)) return amplitude * pow(cos(r / 2.0), 2.0);
+    return 0.0;
+  }
+
+  static VectorField CircleVectorField(const Point &x) {
+    return 2 * Pi * VectorField(-x[1], x[0], 0.0);
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return CircleVectorField(x);
+  };
+
+  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
+                        const Point &x) const override {
+    return CircleVectorField(x) * N;
+  };
+
+  string Name() const override { return "PollutionCosHat2D"; }
+};
+
+class CosHat2D : public IStochasticTransportProblem {
+  double amplitude = 1.0;
+  double frequency = 8.0;
+public:
+  explicit CosHat2D(const Meshes &meshes) : IStochasticTransportProblem(meshes) {}
+
+  void DrawSample(const SampleID &id) override {}
+
+  double Solution(double t, const Point &x) const override {
+    Point midPoint = 0.25 * (Point(cos(2.0 * Pi * t), sin(2.0 * Pi * t)))
+      + Point(0.5, 0.5);
+    double r = dist(midPoint, x);
+    if (r < (Pi / 2.0) / frequency)
+      return amplitude * pow(cos(frequency * r), 2.0);
+    return 0.0;
+  }
+
+  VectorField CircleVectorField(const Point &x) const {
+    return 2 * Pi * VectorField(0.5 - x[1], x[0] - 0.5, 0.0);
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return CircleVectorField(x);
+  };
+
+  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
+                        const Point &x) const override {
+    return CircleVectorField(x) * N;
+  };
+
+  string Name() const override { return "CosHat2D"; }
+};
+
+class StochasticGaussHat2D : public IStochasticTransportProblem {
+private:
+  double mu = 0.0;
+  double sigma = 1.0; // (sigma squared)
+
+  VectorField vectorField;
+
+public:
+  explicit StochasticGaussHat2D(const Meshes &meshes) :
+    vectorField(VectorField(1 / sqrt(2), 1 / sqrt(2), 0)),
+    IStochasticTransportProblem(meshes) {}
+
+  void DrawSample(const SampleID &id) override {}
+
+  Scalar Solution(double t, const Point &x) const override {
+    // Todo 2D Gauss function for t = 0
+    // 1 / (sigma * sqrt(2 * PI)) exp ()
+    if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
+    return 1.0;
+  }
+
+  static VectorField CircleVectorField(const Point &x) {
+    return 2 * Pi * VectorField(-x[1], x[0], 0.0);
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return CircleVectorField(x);
+  }
+
+  Scalar FaceNormalFlux(const cell &c, int face,
+                        const VectorField &N, const Point &x) const override {
+    return CircleVectorField(x) * N;
+  }
+
+  string Name() const override { return "GaussHat2D"; }
+};
+
+class GaussHat2D : public IStochasticTransportProblem {
+private:
+  Tensor sigma;
+
+public:
+  explicit GaussHat2D(const Meshes &meshes) :
+    sigma(Tensor()),
+    IStochasticTransportProblem(meshes) {
+    sigma[0][0] = 0.005;
+    sigma[0][1] = 0.0;
+    sigma[1][0] = 0.0;
+    sigma[1][1] = 0.005;
+    sigma[2][2] = 1.0; // Todo this is needed even in 2D case which is ugly
+  }
+
+  void DrawSample(const SampleID &id) override {}
+
+  static Point mu(double t) {
+    return 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t)) + Point(1.0, 1.0));
+  }
+
+  Scalar Solution(double t, const Point &x) const override {
+    double factor = 1 / sqrt(pow(2 * Pi, 2) * sigma.det());
+    VectorField diff = (x - mu(t));
+    VectorField temp = Invert(sigma) * diff;
+    Scalar exponent = -0.5 * diff * temp;
+    return factor * exp(exponent);
+  }
+
+  static VectorField CircleVectorField(const Point &x) {
+    return 2 * Pi * VectorField(-x[1], x[0], 0.0);
+  }
+
+  VectorField CellFlux(const cell &c, const Point &x) const override {
+    return CircleVectorField(x);
+  }
+
+  Scalar FaceNormalFlux(const cell &c, int face,
+                        const VectorField &N, const Point &x) const override {
+    return CircleVectorField(x) * N;
+  }
+
+  string Name() const override { return "GaussHat2D"; }
+};
+
 IStochasticTransportProblem *
 CreateStochasticTransportProblem(const string &problemName, const Meshes &meshes) {
   if (problemName == "StochasticPollution1D")

mluq/src/problems/StochasticTransportProblem.hpp

@@ -43,354 +43,6 @@ public:
                                 const Point &x) const = 0;
 };
 
-class StochasticPollution1D : public IStochasticTransportProblem {
-  HybridFaceNormalFluxGenerator scalarGenerator;
-  HybridCellFluxGenerator vectorFieldGenerator;
-public:
-  explicit StochasticPollution1D(const Meshes &meshes) :
-    IStochasticTransportProblem(meshes),
-    scalarGenerator(HybridFaceNormalFluxGenerator(meshes)),
-    vectorFieldGenerator(HybridCellFluxGenerator(meshes, scalarGenerator)) {}
-
-  void DrawSample(const SampleID &id) override {
-    scalarGenerator.DrawSample(id);
-    vectorFieldGenerator.DrawSample(id);
-  }
-
-  Scalar Solution(double t, const Point &x) const override {
-    if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
-    return 1.0;
-  }
-
-  VectorField CellFlux(const cell &c, const Point &x) const override {
-    return this->vectorFieldGenerator.EvalSample(c);
-  }
-
-  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
-                        const Point &x) const override {
-    return this->scalarGenerator.EvalSample(face, c);
-  }
-
-  string Name() const override { return "StochasticPollution1D"; }
-};
-
-class Pollution1D : public IStochasticTransportProblem {
-public:
-  explicit Pollution1D(const Meshes &meshes) : IStochasticTransportProblem(meshes) {}
-
-  void DrawSample(const SampleID &id) override {}
-
-  Scalar Solution(double t, const Point &x) const override {
-    if (abs(1.0 * t - x[0] + 0.5) > 0.06251) return 0.0;
-    return 1.0;
-  }
-
-  VectorField CellFlux(const cell &c, const Point &x) const override {
-    return VectorField(-1.0, 0.0);
-  }
-
-  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
-                        const Point &x) const override {
-    return 1.0;
-  }
-
-  string Name() const override { return "Pollution1D"; }
-};
-
-class StochasticPollutionCosHat1D : public IStochasticTransportProblem {
-  HybridFaceNormalFluxGenerator scalarGenerator;
-
-  HybridCellFluxGenerator vectorFieldGenerator;
-
-  double amplitude = 1.00;
-
-  double cc = 6.0;
-
-public:
-  explicit StochasticPollutionCosHat1D(const Meshes &meshes) :
-    IStochasticTransportProblem(meshes),
-    scalarGenerator(HybridFaceNormalFluxGenerator(meshes)),
-    vectorFieldGenerator(HybridCellFluxGenerator(meshes, scalarGenerator)) {}
-
-  void DrawSample(const SampleID &id) override {
-    scalarGenerator.DrawSample(id);
-    vectorFieldGenerator.DrawSample(id);
-  }
-
-  Scalar Solution(double t, const Point &x) const override {
-    Point midPoint = Point(0.2, 0.0);
-    double r = dist(midPoint, x);
-    if (r < 1 / cc)
-      return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
-    return 0.0;
-  }
-
-  VectorField CellFlux(const cell &c, const Point &x) const override {
-    return this->vectorFieldGenerator.EvalSample(c);
-  }
-
-  Scalar FaceNormalFlux(const cell &c, int face,
-                        const VectorField &N, const Point &x) const override {
-    return this->scalarGenerator.EvalSample(face, c);
-  }
-
-  string Name() const override { return "StochasticPollutionCosHat1D"; }
-};
-
-class StochasticPollution2D : public IStochasticTransportProblem {
-  HybridFaceNormalFluxGenerator scalarGenerator;
-
-  HybridCellFluxGenerator vectorFieldGenerator;
-
-  double d = 1.0 / 16.0;
-
-public:
-  explicit StochasticPollution2D(const Meshes &meshes) :
-    IStochasticTransportProblem(meshes),
-    scalarGenerator(HybridFaceNormalFluxGenerator(meshes)),
-    vectorFieldGenerator(HybridCellFluxGenerator(meshes, scalarGenerator)) {}
-
-  void DrawSample(const SampleID &id) override {
-    scalarGenerator.DrawSample(id);
-    vectorFieldGenerator.DrawSample(id);
-  }
-
-  Scalar Solution(double t, const Point &x) const override {
-    if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 3 * d))
-      return 1.0;
-    return 0.0;
-  }
-
-  VectorField CellFlux(const cell &c, const Point &x) const override {
-    return this->vectorFieldGenerator.EvalSample(c);
-  }
-
-  Scalar FaceNormalFlux(const cell &c, int face,
-                        const VectorField &N, const Point &x) const override {
-    return this->scalarGenerator.EvalSample(face, c);
-  }
-
-  string Name() const override { return "StochasticPollution2D"; }
-};
-
-class Pollution2D : public IStochasticTransportProblem {
-private:
-  double d = 1.0 / 16.0;
-
-public:
-  explicit Pollution2D(const Meshes &meshes) : IStochasticTransportProblem(meshes) {}
-
-  void DrawSample(const SampleID &id) override {}
-
-  Scalar Solution(double t, const Point &x) const override {
-    if ((abs(1 - x[1] - 1.5 * d) < d / 2) && (abs(0.5 - x[0]) < 4 * d))
-      return 32.0;
-    return 0.0;
-  }
-
-  static VectorField TransportFlux(const Point &x) {
-    return VectorField(0.0, -1.0);
-  }
-
-  VectorField CellFlux(const cell &c, const Point &x) const override {
-    return TransportFlux(x);
-  }
-
-  Scalar FaceNormalFlux(const cell &c, int face, const VectorField &N,
-                        const Point &x) const override {
-    return TransportFlux(x) * N;
-  }
-
-  string Name() const override { return "Pollution2D"; }
-};
-
-class StochasticPollutionCosHat2D : public IStochasticTransportProblem {
-  HybridFaceNormalFluxGenerator scalarGenerator;
-
-  HybridCellFluxGenerator vectorFieldGenerator;
-
-  double amplitude = 1.00;
-
-  double cc = 6.0;
-
-public:
-  explicit StochasticPollutionCosHat2D(const Meshes &meshes) :
-    IStochasticTransportProblem(meshes),
-    scalarGenerator(HybridFaceNormalFluxGenerator(meshes)),
-    vectorFieldGenerator(HybridCellFluxGenerator(meshes, scalarGenerator)) {}
-
-  void DrawSample(const SampleID &id) override {
-    scalarGenerator.DrawSample(id);
-    vectorFieldGenerator.DrawSample(id);
-  }
-
-  double Solution(double t, const Point &x) const override {
-    Point midPoint = Point(0.5, 0.8);
-    double r = dist(midPoint, x);
-    if (r < 1 / cc)
-      return amplitude * pow(cos(cc * Pi * r) + 1.0, 2.0);
-    return 0.0;
-  }
-
-  VectorField CellFlux(const cell &c, const Point &x) const override {
-    return this->vectorFieldGenerator.EvalSample(c);
-  }
-
-  Scalar FaceNormalFlux(const cell &c, int face,
-                        const VectorField &N, const Point &x) const override {
-    return this->scalarGenerator.EvalSample(face, c);
-  }
-
-  string Name() const override { return "StochasticPollutionCosHat2D"; }
-};
-
-class PollutionCosHat2D : public IStochasticTransportProblem {
-  double amplitude = 2.0 / Pi;
-public:
-  explicit PollutionCosHat2D(const Meshes &meshes)
-    : IStochasticTransportProblem(meshes) {}
-
-  void DrawSample(const SampleID &id) override {}
-
-  double Solution(double t, const Point &x) const override {
-    Point midPoint = 0.25 * (Point(cos(2. * Pi * t), sin(2. * Pi * t))
-      + Point(1.0, 1.0));
-    double r = dist(midPoint, x);
-    if (r < 1.0 / (2.0 * Pi)) return amplitude * pow(cos(r / 2.0), 2.0);
-    return 0.0;
-  }
-
-  static VectorField CircleVectorField(const Point &x) {
-    return 2 * Pi * VectorField(-x[1], x[0], 0.0);
-  }
-
-  VectorField CellFlux(const cell &c, const Point &x) const override {
-    return CircleVectorField(x);
-  };