new assemble idea

mlmc/src/assemble/DGEllipticAssemble.hpp

@@ -3,9 +3,9 @@
 
 #include "discretization/DGDiscretization.hpp"
 #include "elements/Elements.h"
-#include "IEllipticAssemble.hpp"
+#include "IStochasticEllipticAssemble.hpp"
 
-class DGEllipticAssemble : public IEllipticAssemble {
+class DGEllipticAssemble : public IStochasticEllipticAssemble {
 public:
     DGDiscretization *disc;
 
@@ -13,7 +13,7 @@ public:
     int sign = 1;
 
     DGEllipticAssemble(IDiscretization *disc, IStochasticEllipticProblem *problem)
-        : IEllipticAssemble(problem), disc(dynamic_cast<DGDiscretization *> (disc)) {
+        : IStochasticEllipticAssemble(problem), disc(dynamic_cast<DGDiscretization *> (disc)) {
         this->disc = dynamic_cast<DGDiscretizationT<> *>(disc);
         config.get("penalty", penalty);
         config.get("sign", sign);

mlmc/src/assemble/HybridEllipticAssemble.cpp

@@ -545,3 +545,4 @@ void HybridEllipticAssemble::AssembleTransfer(TransferMatrix &TM) const {
         }
     }
 }
+

mlmc/src/assemble/HybridEllipticAssemble.hpp

@@ -10,6 +10,14 @@ public:
     HybridEllipticAssemble(IDiscretization *disc, IStochasticEllipticProblem *problem)
         : MixedEllipticAssemble(disc, problem) {}
 
+    virtual void DrawSample(SampleID id) override {
+        problem->DrawSample(id);
+    };
+
+    virtual IDiscretization *GetDiscretization() override {
+        return disc;
+    }
+
     const char *Name() const override;
 
     void Initialize(Vector &u) const override;

mlmc/src/assemble/IStochasticAssemble.hpp

+#ifndef ISTOCHASTICASSEMBLE_HPP
+#define ISTOCHASTICASSEMBLE_HPP
+
+#include "Assemble.h"
+#include "problem/StochasticProblem.hpp"
+
+// Todo Refactor TAssemble and use ITimeAssemble
+class IStochasticAssemble : public IAssemble {
+protected:
+    int verbose = 1;
+
+public:
+    IStochasticAssemble() {
+        config.get("AssembleVerbose", verbose);
+    }
+
+    const char *Name() const override = 0;
+
+    void Initialize(Vector &u) const override = 0;
+
+    double Energy(const Vector &u) const override = 0;
+
+    void Residual(const cell &c, const Vector &u, Vector &r) const override = 0;
+
+    void Jacobi(const cell &c, const Vector &u, Matrix &A) const override = 0;
+
+    virtual void DrawSample(SampleID id) = 0;
+
+    virtual IDiscretization *GetDiscretization() = 0;
+
+    virtual IStochasticProblem *GetStochasticProblem() = 0;
+
+    void PrintInfo() const override = 0;
+};
+
+#endif //ISTOCHASTICASSEMBLE_HPP

mlmc/src/assemble/IEllipticAssemble.hpp → mlmc/src/assemble/IStochasticEllipticAssemble.hpp

-#ifndef IELLIPTICASSEMBLE_HPP
-#define IELLIPTICASSEMBLE_HPP
+#ifndef ISTOCHASTICELLIPTICASSEMBLE_HPP
+#define ISTOCHASTICELLIPTICASSEMBLE_HPP
 
-#include "Assemble.h"
+#include "assemble/IStochasticAssemble.hpp"
 #include "problem/StochasticEllipticProblem.hpp"
 
 
 typedef std::pair<double, double> FluxPair;
 
-class IEllipticAssemble : public IAssemble {
-protected:
-    int verbose = 1;
-
+class IStochasticEllipticAssemble : public IStochasticAssemble {
 public:
     IStochasticEllipticProblem *problem;
 
-    explicit IEllipticAssemble(IStochasticEllipticProblem *problem) : problem(problem) {
-        config.get("EllipticVerbose", verbose);
+    explicit IStochasticEllipticAssemble(IStochasticEllipticProblem *problem) :
+        IStochasticAssemble(), problem(problem) {
+    };
+
+    void DrawSample(SampleID id) override {
+        problem->DrawSample(id);
+    }
+
+    IStochasticProblem *GetStochasticProblem() override {
+        return problem;
     };
 
-    virtual ~IEllipticAssemble() = default;;
+    virtual ~IStochasticEllipticAssemble() = default;;
 
     const char *Name() const override = 0;
 
@@ -32,6 +37,7 @@ public:
 
     virtual double EnergyError(const Vector &u) const = 0;
 
+    // Todo make class out of it SpacialMeasures see EmpiricalMeasures
     virtual double L2(const Vector &u) const = 0;
 
     virtual double H1(const Vector &u) const = 0;
@@ -104,4 +110,4 @@ public:
     }
 };
 
-#endif //IELLIPTICASSEMBLE_HPP
+#endif //ISTOCHASTICELLIPTICASSEMBLE_HPP

mlmc/src/assemble/LagrangeEllipticAssemble.hpp

@@ -3,14 +3,18 @@
 
 #include "discretization/LagrangeDiscretization.hpp"
 #include "elements/Elements.h"
-#include "IEllipticAssemble.hpp"
+#include "IStochasticEllipticAssemble.hpp"
 
-class LagrangeEllipticAssemble : public IEllipticAssemble {
+class LagrangeEllipticAssemble : public IStochasticEllipticAssemble {
 public:
     IDiscretization *disc;
 
     LagrangeEllipticAssemble(IDiscretization *disc, IStochasticEllipticProblem *problem)
-        : IEllipticAssemble(problem), disc(disc) {}
+        : IStochasticEllipticAssemble(problem), disc(disc) {}
+
+    IDiscretization *GetDiscretization() override {
+        return disc;
+    };
 
     const char *Name() const override;
 

mlmc/src/assemble/MixedEllipticAssemble.hpp

@@ -4,15 +4,20 @@
 #include "discretization/RTLagrangeDiscretization.hpp"
 #include "discretization/LagrangeDiscretization.hpp"
 #include "elements/Elements.h"
-#include "IEllipticAssemble.hpp"
+#include "IStochasticEllipticAssemble.hpp"
 
 
-class MixedEllipticAssemble : public IEllipticAssemble {
-public:
+class MixedEllipticAssemble : public IStochasticEllipticAssemble {
+protected:
     IDiscretization *disc;
 
+public:
     MixedEllipticAssemble(IDiscretization *disc, IStochasticEllipticProblem *problem) :
-        IEllipticAssemble(problem), disc(disc) {};
+        IStochasticEllipticAssemble(problem), disc(disc) {};
+
+    IDiscretization *GetDiscretization() override {
+        return disc;
+    };
 
     const char *Name() const override;
 
@@ -53,44 +58,6 @@ public:
     virtual void SetFlux(const Vector &u, Vector &flux) override;
 
     void AssembleTransfer(TransferMatrix &TM) const override;
-
-//protected:
-//    Element *createElement(const cell &c, const Vector &u) const override {
-//        return new RT0_LagrangeElementT<>(*disc, u, c);
-//    }
-//
-//    FaceElement *createFaceElement(const cell &c, const Vector &u,
-//                                                 int face) const override {
-//        return new RT0FaceElement(*disc, u, c, face);
-//    }
-//
-//    Scalar getValue(const cell &c, const Vector &u, Element *elem, int q) const override {
-//        return u(c(), 0);
-//    }
-//
-//    Scalar getFaceValue(const cell &c, const Vector &u, FaceElement *faceElem,
-//                        int q, int face) const override {
-//        return u(c(), 0);
-//    }
-//
-//    VectorField getFlux(const cell &c, const Vector &u,
-//                        Element *elem, int q) const override {
-//        auto rt0Elem = dynamic_cast<RT0_LagrangeElementT<> *>(elem);
-//        return rt0Elem->VelocityField(q, u);
-//    }
-//
-//    VectorField getFaceFlux(const cell &c, const Vector &u,
-//                            FaceElement *faceElem, int q,
-//                            int face) const override {
-//        auto rt0FaceElem = dynamic_cast<RT0FaceElement *>(faceElem);
-//        return rt0FaceElem->VectorValue(q, u, face);
-//    }
-//
-//    VectorField getGradient(const cell &c, const Vector &u,
-//                            Element *elem, int q) const override {
-//        Tensor IK = problem->InversePermeability(c.Subdomain(), elem->QPoint(q));
-//        return IK * getFlux(c, u, elem, q);
-//    }
 };
 
 #endif

mlmc/src/main/MainProgram.cpp

@@ -102,26 +102,26 @@ IStochasticProblem *MainProgram::createStochasticProblem() {
         return new StochasticLaplace1D(*meshes);
     if (problemName == "StochasticLaplace2D")
         return new StochasticLaplace2D(*meshes);
-    if (problemName == "StochasticPollution1D")
-        return new StochasticPollution1D(*meshes);
-    if (problemName == "StochasticPollutionCosHat1D")
-        return new StochasticPollutionCosHat1D(*meshes);
-    if (problemName == "StochasticPollution2D")
-        return new StochasticPollution2D(*meshes);
-    if (problemName == "StochasticPollutionCosHat2D")
-        return new StochasticPollutionCosHat2D(*meshes);
-    if (problemName == "StochasticPollutionMollifiedBar2D")
-        return new StochasticPollutionMollifiedBar2D(*meshes);
-    if (problemName == "StochasticInitialConditions2D")
-        return nullptr;
+//    if (problemName == "StochasticPollution1D")
+//        return new StochasticPollution1D(*meshes);
+//    if (problemName == "StochasticPollutionCosHat1D")
+//        return new StochasticPollutionCosHat1D(*meshes);
+//    if (problemName == "StochasticPollution2D")
+//        return new StochasticPollution2D(*meshes);
+//    if (problemName == "StochasticPollutionCosHat2D")
+//        return new StochasticPollutionCosHat2D(*meshes);
+//    if (problemName == "StochasticPollutionMollifiedBar2D")
+//        return new StochasticPollutionMollifiedBar2D(*meshes);
+//    if (problemName == "StochasticInitialConditions2D")
+//        return nullptr;
     Exit("\nProblem field not found in " + problemName + "\n")
 }
 
-IAssemble *MainProgram::createAssemble() {
+IStochasticAssemble * MainProgram::createAssemble() {
     IStochasticEllipticProblem *ellipticProblem =
         dynamic_cast<IStochasticEllipticProblem * >(problem);
-    IStochasticTransportProblem *transportProblem =
-        dynamic_cast<IStochasticTransportProblem * >(problem);
+//    IStochasticTransportProblem *transportProblem =
+//        dynamic_cast<IStochasticTransportProblem * >(problem);
 
     if (ellipticProblem != nullptr) {
     if (modelName == "LagrangeFEM")
@@ -134,12 +134,12 @@ IAssemble *MainProgram::createAssemble() {
 //            return new DGEllipticAssemble(disc, ellipticProblem);
     }
 
-    if (transportProblem != nullptr) {
-        if (modelName == "DGTransport")
-            return new DGTransportAssemble(disc, transportProblem, nullptr);
-        if (modelName == "PGTransport")
-            return nullptr; // Todo open research question
-    }
+////    if (transportProblem != nullptr) {
+////        if (modelName == "DGTransport")
+////            return new DGTransportAssemble(disc, transportProblem, nullptr);
+////        if (modelName == "PGTransport")
+////            return nullptr; // Todo open research question
+//    }
 
     Exit(("\nAssembling for Problem=" + problemName + " with Model=" +
         modelName + " not implemented\n"))

mlmc/src/main/MainProgram.hpp

@@ -36,7 +36,7 @@ public:
     Meshes *meshes;
     IDiscretization *disc;
     IStochasticProblem *problem;
-    IAssemble *assemble;
+    IStochasticAssemble *assemble;
     MultilevelMonteCarlo *mlmc;
 
     MainProgram() {
@@ -83,7 +83,7 @@ private:
 
     IStochasticProblem *createStochasticProblem();
 
-    IAssemble *createAssemble();
+    IStochasticAssemble * createAssemble();
 };
 
 #endif //MLMC_MLMCMAIN_H