refactored HybridFluxGenerator

mlmc/src/stochastics/HybridFluxGenerator.C

@@ -2,38 +2,45 @@
 
 
 void HybridFluxGenerator::generateFineSample(Vector &fineField) {
-    Vector &normalFlux = fineField;
-    Vector flux(cellMg3[l - meshes.pLevel()]);
-    Vector perm(cellMg1[l - meshes.pLevel()]);
-    Vector u(faceMg[l - meshes.pLevel()]);
+    initialize(true);
+    permeabilityGenerator->GetFineSample(*finePerm);
+    assemble->problem->LoadNewSample(finePerm);
 
-    permeabilityGenerator->GetFineSample(perm);
-    assemble->problem->LoadNewSample(&perm);
-
-    Preconditioner *pc = GetPC(faceMg, *assemble, "SuperLU");
+    Preconditioner *pc = GetPC("SuperLU");
     Solver solver(pc, "GMRES");
     Newton newton(solver);
 
-    newton(*assemble, u);
-    assemble->setFlux(u, flux);
-    assemble->SetNormalFlux(u, normalFlux);
+    (newton)(*assemble, *u);
+    assemble->setFlux(*u, *flux);
+    assemble->SetNormalFlux(*u, fineField);
 
-    if (plotting) generatePlot(flux);
+    if (plotting) generatePlot(*flux);
 }
 
 void HybridFluxGenerator::generateCoarseSample(const Vector &fineField,
                                                Vector &coarseField) {
-    coarseField = 0;
-    for (cell c = coarseField.cells(); c != coarseField.cells_end(); c++) {
-        row coarseRow = coarseField.find_row(c());
-        for (int k = 0; k < c.Children(); ++k) {
-            row fineRow = fineField.find_row(c.Child(k));
-            for (int i = 0; i < coarseRow.n(); ++i)
-                coarseField(coarseRow, i) += fineField(fineRow, i);
-        }
-        for (int i = 0; i < coarseRow.n(); ++i)
-            coarseField(coarseRow, i) /= c.Children();
+    initialize(false);
+    permeabilityGenerator->GetCoarseSample(*finePerm, *coarsePerm);
+    assemble->problem->LoadNewSample(coarsePerm);
+
+    Preconditioner *pc = GetPC("SuperLU");
+    Solver solver(pc, "GMRES");
+    Newton newton(solver);
+
+    (newton)(*assemble, *u);
+    assemble->setFlux(*u, *flux);
+    assemble->SetNormalFlux(*u, coarseField);
+
+    if (plotting) generatePlot(*flux);
+}
+
+void HybridFluxGenerator::initialize(bool fineField) {
+    if (fineField) {
+        flux = new Vector(cellMatrixGraphs[l - meshes.pLevel()]);
+        u = new Vector(faceMatrixGraphs[l - meshes.pLevel()]);
+    } else {
+        flux = new Vector(cellMatrixGraphs[l - meshes.pLevel() - 1]);
+        u = new Vector(faceMatrixGraphs[l - meshes.pLevel() - 1]);
     }
-    // TOD0
-//    if (plotting) generatePlot(flux);
-}
+    *flux = 0.0, *u = 0.0;
+}
\ No newline at end of file

mlmc/src/stochastics/HybridFluxGenerator.h

@@ -8,9 +8,19 @@
 
 class HybridFluxGenerator : public SampleGenerator {
 private:
-    MatrixGraphs cellMg1;
-    MatrixGraphs cellMg3;
-    MatrixGraphs faceMg;
+    MatrixGraphs cellMatrixGraphs;
+    MatrixGraphs faceMatrixGraphs;
+
+//    Preconditioner *pc;
+//    Solver *solver;
+//    Newton *newton;
+
+    Vector *finePerm = nullptr;
+    Vector *coarsePerm = nullptr;
+    Vector *flux = nullptr;
+    Vector *u = nullptr;
+
+    void initialize(bool fineField);
 
 protected:
     Discretization *disc = nullptr;
@@ -31,9 +41,13 @@ protected:
 public:
     explicit HybridFluxGenerator(int l, Meshes &meshes)
         : SampleGenerator(l, meshes),
-          cellMg1(MatrixGraphs(meshes, dof("cell", 1))),
-          cellMg3(MatrixGraphs(meshes, dof("cell", 3))),
-          faceMg(MatrixGraphs(meshes, dof("face", 1))) {
+          cellMatrixGraphs(MatrixGraphs(meshes, dof("cell", 3))),
+          faceMatrixGraphs(MatrixGraphs(meshes, dof("face", 1))) {
+        finePerm = new Vector(cellMatrixGraphs[l - meshes.pLevel()]);
+        if (l != meshes.pLevel())
+            coarsePerm = new Vector(cellMatrixGraphs[l - meshes.pLevel() - 1]);
+        else
+            coarsePerm = nullptr;
     }
 
     ~HybridFluxGenerator() {
@@ -41,7 +55,12 @@ public:
         delete problem;
         delete assemble;
         delete permeabilityGenerator;
-    };
+
+        delete finePerm;
+        delete coarsePerm;
+        delete flux;
+        delete u;
+    }
 
     void SetPlot(Plot *plotFromMC) override {
         permeabilityGenerator->SetPlot(plotFromMC);