new HybridFluxGenerator

mlmc/src/stochastics/HybridFluxGenerator.cpp

 #include "HybridFluxGenerator.hpp"
-
-
-template<typename Sample>
-void HybridFluxGenerator<Sample>::generateFineSample(SampleID id,
-                                                     Vector *fineSample,
-                                                     Vector *coarseSample) {
-    initialize(true);
-    assemble->problem->DrawSample(id);
-    newton->operator()(*assemble, *u);
-    assemble->SetFlux(*u, *flux);
-    assemble->SetNormalFlux(*u, *fineSample);
-//    if (plotting)
-//        plotter->PlotVector("flux", *flux,
-//                            3, l, "CellData");
-}
-
-template<typename Sample>
-void HybridFluxGenerator<Sample>::generateCoarseSample(SampleID id,
-                                                       Vector *fineSample,
-                                                       Vector *coarseSample) {
-    initialize(true);
-    assemble->problem->DrawSample(id);
-    newton->operator()(*assemble, *u);
-    assemble->SetFlux(*u, *flux);
-    assemble->SetNormalFlux(*u, *coarseSample);
-//    if (plotting)
-//        plotter->PlotVector("flux", *flux,
-//                            3, l - 1, "CellData");
-}
-
-template<typename Sample>
-void HybridFluxGenerator<Sample>::initialize(bool fineField) {
-    if (fineField) {
-        flux = new Vector(cellMatrixGraphs[this->l - this->meshes.pLevel()]);
-        u = new Vector(faceMatrixGraphs[this->l - this->meshes.pLevel()]);
-    } else {
-        flux = new Vector(cellMatrixGraphs[this->l - this->meshes.pLevel() - 1]);
-        u = new Vector(faceMatrixGraphs[this->l - this->meshes.pLevel() - 1]);
-    }
-    *flux = 0.0, *u = 0.0;
-}
-
-
-//HybridMainForPollutionProblem::HybridMainForPollutionProblem(Meshes &meshes) {
-//    config.get("HybridVerbose", verbose);
-//    Initialize(meshes);
-//    PrintInfo();
-//    ResetVectors();
-//    Solve();
-//}
-//
-//void HybridMainForPollutionProblem::Initialize(Meshes &meshes) {
-//    disc = std::make_shared<RTLagrangeDiscretization>(meshes, 0, 0);
-//    problem = std::make_shared<Laplace2D>();
-//    assemble = std::make_shared<HybridEllipticAssemble>(disc.get(), problem.get());
-//    mGraphs = std::make_shared<MatrixGraphs>(meshes, dof("face", 1));
-//    pc = std::shared_ptr<Preconditioner>(GetPC("SuperLU"));
-//    solver = std::make_shared<Solver>(pc.get(), "GMRES");
-//    newton = std::make_shared<Newton>(*solver);
-//    u = std::make_shared<Vector>(mGraphs->fine());
-//    normalFlux = std::make_shared<Vector>(mGraphs->fine());
-//}
-//
-//void HybridMainForPollutionProblem::PrintInfo() const {
-//    if (verbose == 0) return;
-//    assemble->PrintInfo();
-//    solver->PrintInfo();
-//}
-//
-//void HybridMainForPollutionProblem::ResetVectors() const {
-//    (*u) = 0.0;
-//    (*normalFlux) = 0.0;
-//}
-//
-//void HybridMainForPollutionProblem::Solve() const {
-//    mout.StartBlock("Hybrid Flux Problem");
-//    mout << "Solve" << endl;
-//    (*newton)(*assemble, *u);
-//    assemble->SetNormalFlux(*u, *normalFlux);
-//    mout.EndBlock();
-//    mout << endl;
-//}
-//
-//void HybridMainForPollutionProblem::EquipHybridFluxTransport(
-//    StochasticHybridFluxTransport *hybridFluxTransport) {
-//    hybridFluxTransport->SetHybridEllipticAssemble(assemble.get());
-//    hybridFluxTransport->SetFlux(normalFlux.get());
-//}

mlmc/src/stochastics/HybridFluxGenerator.hpp

@@ -6,80 +6,110 @@
 #include "problem/StochasticEllipticProblem.hpp"
 #include "assemble/HybridEllipticAssemble.hpp"
 #include "problem/StochasticHybridFlux.hpp"
-#include "problem/StochasticTransportProblem.hpp"
 #include "discretization/RTLagrangeDiscretization.hpp"
 #include "solver/Newton.h"
 #include <memory>
 
 
-template<typename Sample>
-class HybridFluxGenerator : public SampleGenerator<Sample> {
+class MultilevelHybridFluxGenerator : public SampleGenerator {
 private:
-    MatrixGraphs cellMatrixGraphs;
-    MatrixGraphs faceMatrixGraphs;
+    MatrixGraphs cellMGraphs;
+    MatrixGraphs faceMGraphs;
 
-    Vector *finePerm;
-    Vector *coarsePerm;
-    Vector *flux;
-    Vector *u;
-
-    void initialize(bool fineField);
+    Vector *faceFlux = nullptr;
+    Vector *faceValues = nullptr;
+    Vector *cellFlux = nullptr;
 
 protected:
-    IStochasticEllipticProblemT<Sample> *problem;
+    IStochasticEllipticProblem *problem;
     HybridEllipticAssemble *assemble;
     RTLagrangeDiscretization *disc;
     Preconditioner *pc;
     Solver *solver;
     Newton *newton;
-
-    void generateFineSample(SampleID id,
-                            Vector *fineSample,
-                            Vector *coarseSample) override;
-
-    void generateCoarseSample(SampleID id,
-                              Vector *fineSample,
-                              Vector *coarseSample) override;
+    Meshes &meshes;
+
+    void drawSample(SampleID id) override {
+        if (!id.coarse) {
+            generateFineSample(id);
+            if (plotting)
+                plotter->PlotVector("flux", *cellFlux, 3,
+                                    id.level, "CellData");
+        } else {
+            generateCoarseSample(id);
+            if (plotting)
+                plotter->PlotVector("flux", *cellFlux, 3,
+                                    id.level, "CellData");
+        }
+    }
 
 public:
-    explicit HybridFluxGenerator(int l, Meshes &meshes) :
-        SampleGenerator<Sample>(l, meshes),
-        cellMatrixGraphs(MatrixGraphs(meshes, dof("cell", 3))),
-        faceMatrixGraphs(MatrixGraphs(meshes, dof("face", 1))),
+    explicit MultilevelHybridFluxGenerator(Meshes &meshes) :
+        SampleGenerator(meshes),
+        meshes(meshes),
+        cellMGraphs(MatrixGraphs(meshes, dof(new CellDoF(3)))),
+        faceMGraphs(MatrixGraphs(meshes, dof(new FaceDoF(1)))),
         disc(new RTLagrangeDiscretization(meshes, 0, 0)),
         pc(GetPC("SuperLU")),
         solver(new Solver(pc, "GMRES")),
         newton(new Newton(*solver)) {
+        if (meshes.dim() == 1)
+            problem = new StochasticLaplace1D(meshes);
+        assemble = new HybridEllipticAssemble(disc, problem);
+        if (meshes.dim() == 2)
+            problem = new StochasticLaplace2D(meshes);
+        assemble = new HybridEllipticAssemble(disc, problem);
+    }
 
-        finePerm = new Vector(cellMatrixGraphs[l - meshes.pLevel()]);
-        if (l != meshes.pLevel())
-            coarsePerm = new Vector(cellMatrixGraphs[l - meshes.pLevel() - 1]);
+    virtual ~MultilevelHybridFluxGenerator() {
+        delete problem;
+        delete assemble;
+        delete disc;
+        delete solver;
+        delete newton;
     }
 
-    virtual ~HybridFluxGenerator() {
+    string Name() const override { return "Hybrid Flux Generator"; }
 
-    }
+    void generateFineSample(SampleID id) {
+        cellFlux = new Vector(cellMGraphs[id.level - this->meshes.pLevel()]);
+        faceValues = new Vector(faceMGraphs[id.level - this->meshes.pLevel()]);
+        faceFlux = new Vector(faceMGraphs[id.level - this->meshes.pLevel()]);
 
-    void PrintInfo() const;
-};
+        assemble->problem->DrawSample(id);
+        newton->operator()(*assemble, *faceValues);
+        assemble->SetNormalFlux(*faceValues, *faceFlux);
+        assemble->SetFlux(*faceValues, *cellFlux);
+    }
 
-template<typename Sample=Vector>
-class HybridFluxGenerator1D : public HybridFluxGenerator<Sample> {
-public:
-    explicit HybridFluxGenerator1D(int l, Meshes &meshes)
-        : HybridFluxGenerator<Sample>(l, meshes) {
-        this->problem = new StochasticLaplace1D(meshes);
-        this->assemble = new HybridEllipticAssemble(this->disc, this->problem);
+    void generateCoarseSample(SampleID id) {
+        cellFlux = new Vector(cellMGraphs[id.level - 1 - this->meshes.pLevel()]);
+        faceValues = new Vector(faceMGraphs[id.level - 1 - this->meshes.pLevel()]);
+        faceFlux = new Vector(faceMGraphs[id.level - 1 - this->meshes.pLevel()]);
+
+        assemble->problem->DrawSample(id);
+        newton->operator()(*assemble, *faceValues);
+        assemble->SetNormalFlux(*faceValues, *faceFlux);
+        assemble->SetFlux(*faceValues, *cellFlux);
+    };
+
+    VectorField EvalVectorFieldSample(const cell &c) override {
+        RT0_LagrangeElementT<> elem(*disc, *faceFlux, c);
+        VectorField F = zero;
+        double area = 0;
+        for (int q = 0; q < elem.nQ(); ++q) {
+            double w = elem.QWeight(q);
+            area += w;
+            F += w * elem.VelocityField(q, *faceFlux);
+        }
+        F *= (1 / area);
+        return F;
     }
-};
 
-template<typename Sample=Vector>
-class HybridFluxGenerator2D : public HybridFluxGenerator<Sample> {
-public:
-    explicit HybridFluxGenerator2D(int l, Meshes &meshes)
-        : HybridFluxGenerator<Sample>(l, meshes) {
-        this->problem = new StochasticLaplace2D(meshes);
-        this->assemble = new HybridEllipticAssemble(this->disc, this->problem);
+    Scalar EvalScalarSample(int face, const cell &c) override {
+        RT0_LagrangeElementT<> elem(*disc, *faceFlux, c);
+        RT0_LagrangeFaceElementT<> faceElem(*disc, *faceFlux, c, face);
+        return (*faceFlux)(elem[face], 0) * elem.Sign(face) / faceElem.Area();
     }
 };