HybridFluxGenerator.hpp 4 KB
Newer Older
1
2
#ifndef HYBRIDFLUXGENERATOR_HPP
#define HYBRIDFLUXGENERATOR_HPP
3

4
5
#include "SampleGenerator.hpp"
#include "CirculantEmbedding.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
6
7
#include "problem/StochasticHybridFlux.hpp"
#include <memory>
8

9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
#include "montecarlo/PDESolver.hpp"


class HybridFluxGenerator : public SampleGenerator {
private:
    PDESolver *pdeSolver;

public:
    HybridFluxGenerator(Meshes &meshes) : SampleGenerator(meshes) {
        std::string problemName = "StochasticLaplace2D";
        pdeSolver = PDESolverCreator(meshes).
            WithModel("HybridElliptic").
            WithProblem(problemName).Create();
    }
};
24

niklas.baumgarten's avatar
niklas.baumgarten committed
25
class MultilevelHybridFluxGenerator : public SampleGenerator {
niklas.baumgarten's avatar
niklas.baumgarten committed
26
private:
niklas.baumgarten's avatar
niklas.baumgarten committed
27
28
    MatrixGraphs cellMGraphs;
    MatrixGraphs faceMGraphs;
29

niklas.baumgarten's avatar
niklas.baumgarten committed
30
31
32
    Vector *faceFlux = nullptr;
    Vector *faceValues = nullptr;
    Vector *cellFlux = nullptr;
niklas.baumgarten's avatar
niklas.baumgarten committed
33

34
protected:
niklas.baumgarten's avatar
niklas.baumgarten committed
35
    IStochasticEllipticProblem *problem;
36
37
38
39
40
    HybridEllipticAssemble *assemble;
    RTLagrangeDiscretization *disc;
    Preconditioner *pc;
    Solver *solver;
    Newton *newton;
niklas.baumgarten's avatar
niklas.baumgarten committed
41
42
    Meshes &meshes;

43
    void drawSample(const SampleID &id) override {
niklas.baumgarten's avatar
niklas.baumgarten committed
44
45
        if (!id.coarse) {
            generateFineSample(id);
niklas.baumgarten's avatar
niklas.baumgarten committed
46
47
48
//            if (plotting)
//                plotter->PlotVector("flux", *cellFlux, 3,
//                                    id.level.fine, "CellData");
niklas.baumgarten's avatar
niklas.baumgarten committed
49
50
        } else {
            generateCoarseSample(id);
niklas.baumgarten's avatar
niklas.baumgarten committed
51
52
53
//            if (plotting)
//                plotter->PlotVector("flux", *cellFlux, 3,
//                                    id.level.coarse, "CellData");
niklas.baumgarten's avatar
niklas.baumgarten committed
54
55
        }
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
56

57
public:
niklas.baumgarten's avatar
niklas.baumgarten committed
58
59
60
61
62
    explicit MultilevelHybridFluxGenerator(Meshes &meshes) :
        SampleGenerator(meshes),
        meshes(meshes),
        cellMGraphs(MatrixGraphs(meshes, dof(new CellDoF(3)))),
        faceMGraphs(MatrixGraphs(meshes, dof(new FaceDoF(1)))),
63
64
65
66
        disc(new RTLagrangeDiscretization(meshes, 0, 0)),
        pc(GetPC("SuperLU")),
        solver(new Solver(pc, "GMRES")),
        newton(new Newton(*solver)) {
67
        if (meshes.dim() == 1) {
niklas.baumgarten's avatar
niklas.baumgarten committed
68
            problem = new StochasticLaplace1D(meshes);
69
70
71
            assemble = new HybridEllipticAssemble(disc, problem);
        }
        if (meshes.dim() == 2) {
niklas.baumgarten's avatar
niklas.baumgarten committed
72
            problem = new StochasticLaplace2D(meshes);
73
74
            assemble = new HybridEllipticAssemble(disc, problem);
        }
niklas.baumgarten's avatar
niklas.baumgarten committed
75
    }
76

niklas.baumgarten's avatar
niklas.baumgarten committed
77
78
79
80
81
82
    virtual ~MultilevelHybridFluxGenerator() {
        delete problem;
        delete assemble;
        delete disc;
        delete solver;
        delete newton;
83
    }
84

niklas.baumgarten's avatar
niklas.baumgarten committed
85
    string Name() const override { return "Hybrid Flux Generator"; }
86

niklas.baumgarten's avatar
niklas.baumgarten committed
87
    void generateFineSample(SampleID id) {
88
89
90
        cellFlux = new Vector(cellMGraphs[id.level.mGfine]);
        faceValues = new Vector(faceMGraphs[id.level.mGfine]);
        faceFlux = new Vector(faceMGraphs[id.level.mGfine]);
91

92
        assemble->DrawSample(id);
niklas.baumgarten's avatar
niklas.baumgarten committed
93
94
95
96
        newton->operator()(*assemble, *faceValues);
        assemble->SetNormalFlux(*faceValues, *faceFlux);
        assemble->SetFlux(*faceValues, *cellFlux);
    }
97

niklas.baumgarten's avatar
niklas.baumgarten committed
98
    void generateCoarseSample(SampleID id) {
99
100
101
        cellFlux = new Vector(cellMGraphs[id.level.mGcoarse]);
        faceValues = new Vector(faceMGraphs[id.level.mGcoarse]);
        faceFlux = new Vector(faceMGraphs[id.level.mGcoarse]);
niklas.baumgarten's avatar
niklas.baumgarten committed
102

103
        assemble->DrawSample(id);
niklas.baumgarten's avatar
niklas.baumgarten committed
104
105
106
107
108
109
        newton->operator()(*assemble, *faceValues);
        assemble->SetNormalFlux(*faceValues, *faceFlux);
        assemble->SetFlux(*faceValues, *cellFlux);
    };

    VectorField EvalVectorFieldSample(const cell &c) override {
niklas.baumgarten's avatar
niklas.baumgarten committed
110
        RTLagrangeElement elem(*disc, *faceFlux, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
111
112
113
114
115
116
117
118
119
        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;
120
    }
121

niklas.baumgarten's avatar
niklas.baumgarten committed
122
    Scalar EvalScalarSample(int face, const cell &c) override {
niklas.baumgarten's avatar
niklas.baumgarten committed
123
        RTLagrangeElement elem(*disc, *faceFlux, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
124
        RTFaceElementT<> faceElem(*disc, *faceFlux, c, face);
niklas.baumgarten's avatar
niklas.baumgarten committed
125
        return (*faceFlux)(elem[face], 0) * elem.Sign(face) / faceElem.Area();
126
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
127
128
};

129
#endif //HYBRIDFLUXGENERATOR_HPP