try to feed PDESolver with Assemble

mlmc/src/montecarlo/MonteCarlo.cpp

@@ -27,5 +27,5 @@ void MonteCarlo::findSampleSolution(int m, SampleID &id, SampleSolution &solutio
     solution.id = id;
     plotter->SetDirectory(id);
     assemble->DrawSample(id);
-    pdeSolver->Run(assemble, solution);
+    pdeSolver->Run(solution);
 }

mlmc/src/montecarlo/PDESolver.cpp

@@ -2,18 +2,20 @@
 #include "assemble/LagrangeEllipticAssemble.hpp"
 
 
-void EllipticPDESolver::run(IAssemble *assemble, SampleSolution &solution) {
+void EllipticPDESolver::run(SampleSolution &solution) {
     newton(*assemble, solution.U);
 }
 
-void EllipticPDESolver::computeQ(IAssemble *assemble, SampleSolution &solution) {
+void EllipticPDESolver::computeQ(SampleSolution &solution) {
     IStochasticEllipticAssemble *ellAssemble =
         dynamic_cast<IStochasticEllipticAssemble *>(assemble);
     if (quantity == "L2") solution.Q = ellAssemble->L2(solution.U);
     else if (quantity == "H1") solution.Q = ellAssemble->H1(solution.U);
     else if (quantity == "Energy") solution.Q = ellAssemble->Energy(solution.U);
-    else if (quantity == "Inflow") solution.Q = ellAssemble->InflowOutflow(solution.U).first;
-    else if (quantity == "Outflow") solution.Q = ellAssemble->InflowOutflow(solution.U).second;
+    else if (quantity == "Inflow")
+        solution.Q = ellAssemble->InflowOutflow(solution.U).first;
+    else if (quantity == "Outflow")
+        solution.Q = ellAssemble->InflowOutflow(solution.U).second;
     else Exit("Quantity of interest not implemented")
 }
 
@@ -26,12 +28,12 @@ void EllipticPDESolver::computeCost(SampleSolution &solution) {
 void EllipticPDESolver::plotSolution(SampleSolution &solution) {
 }
 
-void TransportPDESolver::run(IAssemble *assemble, SampleSolution &solution) {
+void TransportPDESolver::run(SampleSolution &solution) {
     TimeSeries timeSeries;
-    timeInt(dynamic_cast<TAssemble *>(assemble), timeSeries, solution.U);
+//    timeInt(assemble, timeSeries, solution.U);
 }
 
-void TransportPDESolver::computeQ(IAssemble *assemble, SampleSolution &solution) {
+void TransportPDESolver::computeQ(SampleSolution &solution) {
 }
 
 void TransportPDESolver::computeCost(SampleSolution &solution) {

mlmc/src/montecarlo/PDESolver.hpp

@@ -5,6 +5,8 @@
 #include "solver/Solver.h"
 #include "solver/Newton.h"
 #include "timestepping/TimeIntegrator.hpp"
+#include "assemble/IStochasticAssemble.hpp"
+#include "assemble/IStochasticEllipticAssemble.hpp"
 
 
 class PDESolver {
@@ -15,9 +17,9 @@ protected:
 
     string costMeasure = "size";
 
-    virtual void run(IAssemble *assemble, SampleSolution &solution) = 0;
+    virtual void run(SampleSolution &solution) = 0;
 
-    virtual void computeQ(IAssemble *assemble, SampleSolution &solution) = 0;
+    virtual void computeQ(SampleSolution &solution) = 0;
 
     virtual void computeCost(SampleSolution &solution) = 0;
 
@@ -30,16 +32,22 @@ public:
         config.get("CostMeasure", costMeasure);
     }
 
-    void Run(IAssemble *assemble, SampleSolution &solution) {
+    void Run(SampleSolution &solution) {
         mout.StartBlock("PDESolver");
         vout(1) << solution.id.Str() << endl;
-        run(assemble, solution);
-        computeQ(assemble, solution);
+        run(solution);
+        computeQ(solution);
         computeCost(solution);
         plotSolution(solution);
         vout(1) << "Q=" << solution.Q << " cost=" << solution.Cost << endl;
         mout.EndBlock(verbose == 0);
     }
+
+    virtual IDiscretization *GetDiscretization() = 0;
+
+    virtual void DrawSample(SampleID id) = 0;
+
+    virtual IStochasticAssemble *GetAssemble() = 0;
 };
 
 class EllipticPDESolver : public PDESolver {
@@ -48,20 +56,35 @@ private:
 
     Newton newton;
 
+    IStochasticEllipticAssemble *assemble;
+
 protected:
-    void run(IAssemble *assemble, SampleSolution &solution) override;
+    void run(SampleSolution &solution) override;
 
-    void computeQ(IAssemble *assemble, SampleSolution &solution) override;
+    void computeQ(SampleSolution &solution) override;
 
     void computeCost(SampleSolution &solution) override;
 
     void plotSolution(SampleSolution &solution) override;
 
 public:
-    EllipticPDESolver() :
+    EllipticPDESolver(IStochasticEllipticAssemble *stochasticAssemble) :
+        assemble(stochasticAssemble),
         PDESolver(),
         solver(Solver(GetPC("SuperLU"), "GMRES")),
         newton(Newton(solver)) {}
+
+    IDiscretization *GetDiscretization() override {
+        return assemble->GetDiscretization();
+    }
+
+    void DrawSample(SampleID id) override {
+        assemble->DrawSample(id);
+    }
+
+    IStochasticAssemble *GetAssemble() override {
+        return assemble;
+    }
 };
 
 class TransportPDESolver : public PDESolver {
@@ -70,17 +93,20 @@ private:
 
     TimeIntegrator timeInt;
 
+    IStochasticAssemble *assemble;
+
 protected:
-    void run(IAssemble *assemble, SampleSolution &solution) override;
+    void run(SampleSolution &solution) override;
 
-    void computeQ(IAssemble *assemble, SampleSolution &solution) override;
+    void computeQ(SampleSolution &solution) override;
 
     void computeCost(SampleSolution &solution) override;
 
     void plotSolution(SampleSolution &solution) override;
 
 public:
-    TransportPDESolver() :
+    TransportPDESolver(IStochasticAssemble *stochasticAssemble) :
+        assemble(stochasticAssemble),
         PDESolver(),
         solver(Solver(GetPC("PointBlockGaussSeidel"), "GMRES")),
         timeInt(TimeIntegrator(solver)) {}