HybridFluxGenerator.hpp 3.78 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
#include "problem/StochasticEllipticProblem.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
7
8
#include "assemble/HybridEllipticAssemble.hpp"
#include "problem/StochasticHybridFlux.hpp"
9
#include "discretization/RTLagrangeDiscretization.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
10
#include "solver/Newton.hpp"
niklas.baumgarten's avatar
niklas.baumgarten committed
11
#include <memory>
12
13


niklas.baumgarten's avatar
niklas.baumgarten committed
14
class MultilevelHybridFluxGenerator : public SampleGenerator {
niklas.baumgarten's avatar
niklas.baumgarten committed
15
private:
niklas.baumgarten's avatar
niklas.baumgarten committed
16
17
    MatrixGraphs cellMGraphs;
    MatrixGraphs faceMGraphs;
18

niklas.baumgarten's avatar
niklas.baumgarten committed
19
20
21
    Vector *faceFlux = nullptr;
    Vector *faceValues = nullptr;
    Vector *cellFlux = nullptr;
niklas.baumgarten's avatar
niklas.baumgarten committed
22

23
protected:
niklas.baumgarten's avatar
niklas.baumgarten committed
24
    IStochasticEllipticProblem *problem;
25
26
27
28
29
    HybridEllipticAssemble *assemble;
    RTLagrangeDiscretization *disc;
    Preconditioner *pc;
    Solver *solver;
    Newton *newton;
niklas.baumgarten's avatar
niklas.baumgarten committed
30
31
    Meshes &meshes;

32
    void drawSample(const SampleID &id) override {
niklas.baumgarten's avatar
niklas.baumgarten committed
33
34
        if (!id.coarse) {
            generateFineSample(id);
niklas.baumgarten's avatar
niklas.baumgarten committed
35
36
37
//            if (plotting)
//                plotter->PlotVector("flux", *cellFlux, 3,
//                                    id.level.fine, "CellData");
niklas.baumgarten's avatar
niklas.baumgarten committed
38
39
        } else {
            generateCoarseSample(id);
niklas.baumgarten's avatar
niklas.baumgarten committed
40
41
42
//            if (plotting)
//                plotter->PlotVector("flux", *cellFlux, 3,
//                                    id.level.coarse, "CellData");
niklas.baumgarten's avatar
niklas.baumgarten committed
43
44
        }
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
45

46
public:
niklas.baumgarten's avatar
niklas.baumgarten committed
47
48
49
50
51
    explicit MultilevelHybridFluxGenerator(Meshes &meshes) :
        SampleGenerator(meshes),
        meshes(meshes),
        cellMGraphs(MatrixGraphs(meshes, dof(new CellDoF(3)))),
        faceMGraphs(MatrixGraphs(meshes, dof(new FaceDoF(1)))),
52
53
54
55
        disc(new RTLagrangeDiscretization(meshes, 0, 0)),
        pc(GetPC("SuperLU")),
        solver(new Solver(pc, "GMRES")),
        newton(new Newton(*solver)) {
56
        if (meshes.dim() == 1) {
niklas.baumgarten's avatar
niklas.baumgarten committed
57
            problem = new StochasticLaplace1D(meshes);
58
59
60
            assemble = new HybridEllipticAssemble(disc, problem);
        }
        if (meshes.dim() == 2) {
niklas.baumgarten's avatar
niklas.baumgarten committed
61
            problem = new StochasticLaplace2D(meshes);
62
63
            assemble = new HybridEllipticAssemble(disc, problem);
        }
niklas.baumgarten's avatar
niklas.baumgarten committed
64
    }
65

niklas.baumgarten's avatar
niklas.baumgarten committed
66
67
68
69
70
71
    virtual ~MultilevelHybridFluxGenerator() {
        delete problem;
        delete assemble;
        delete disc;
        delete solver;
        delete newton;
72
    }
73

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

niklas.baumgarten's avatar
niklas.baumgarten committed
76
    void generateFineSample(SampleID id) {
77
78
79
        cellFlux = new Vector(cellMGraphs[id.level.mGfine]);
        faceValues = new Vector(faceMGraphs[id.level.mGfine]);
        faceFlux = new Vector(faceMGraphs[id.level.mGfine]);
80

81
        assemble->DrawSample(id);
niklas.baumgarten's avatar
niklas.baumgarten committed
82
83
84
85
        newton->operator()(*assemble, *faceValues);
        assemble->SetNormalFlux(*faceValues, *faceFlux);
        assemble->SetFlux(*faceValues, *cellFlux);
    }
86

niklas.baumgarten's avatar
niklas.baumgarten committed
87
    void generateCoarseSample(SampleID id) {
88
89
90
        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
91

92
        assemble->DrawSample(id);
niklas.baumgarten's avatar
niklas.baumgarten committed
93
94
95
96
97
98
        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
99
        RTLagrangeElement elem(*disc, *faceFlux, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
100
101
102
103
104
105
106
107
108
        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;
109
    }
110

niklas.baumgarten's avatar
niklas.baumgarten committed
111
    Scalar EvalScalarSample(int face, const cell &c) override {
niklas.baumgarten's avatar
niklas.baumgarten committed
112
        RTLagrangeElement elem(*disc, *faceFlux, c);
niklas.baumgarten's avatar
niklas.baumgarten committed
113
        RTFaceElementT<> faceElem(*disc, *faceFlux, c, face);
niklas.baumgarten's avatar
niklas.baumgarten committed
114
        return (*faceFlux)(elem[face], 0) * elem.Sign(face) / faceElem.Area();
115
    }
niklas.baumgarten's avatar
niklas.baumgarten committed
116
117
};

118
#endif //HYBRIDFLUXGENERATOR_HPP