HybridFluxGenerator.cpp 2.16 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
#include "HybridFluxGenerator.hpp"
#include "pdesolver/PDESolver.hpp"


void HybridFaceNormalFluxGenerator::createPDESolver() {
    string problemName;
    if (meshes.dim() == 1)
        problemName = "StochasticLaplace1D";
    if (meshes.dim() == 2)
        problemName = "StochasticLaplace2D";

    pdeSolver = PDESolverCreator(meshes).
        WithModel("HybridElliptic").
        WithProblem(problemName).Create();
}

void HybridFaceNormalFluxGenerator::createMGraphs() {
    faceMGraphs = pdeSolver->CreateSolutionMatrixGraphs(meshes);
    cellMGraphs = new MatrixGraphs(meshes, dof(new CellDoF(3)));
}

void HybridFaceNormalFluxGenerator::drawSample(const SampleID &id) {
    pdeSolver->DrawSample(id);
    solutionFaceValues = new SampleSolution(faceMGraphs, id);
    solutionFaceFlux = new SampleSolution(faceMGraphs, id);
    solutionCellFlux = new SampleSolution(cellMGraphs, id, "Flux");
    pdeSolver->Run(*solutionFaceValues);
    auto assemble = dynamic_cast<HybridEllipticAssemble *>(pdeSolver->GetAssemble());
    assemble->SetNormalFlux(solutionFaceValues->U, solutionFaceFlux->U);
    assemble->SetFlux(solutionFaceValues->U, solutionCellFlux->U);
    if (plotting) plotMap.VtkPlot(*solutionCellFlux);
}

Scalar HybridFaceNormalFluxGenerator::EvalSample(int face, const cell &c) {
    RTLagrangeElement elem(*pdeSolver->GetDisc(), solutionFaceFlux->U, c);
    RTFaceElement faceElem(*pdeSolver->GetDisc(), solutionFaceFlux->U, c, face);
    return (solutionFaceFlux->U)(elem[face], 0) * elem.Sign(face) / faceElem.Area();
}

HybridFaceNormalFluxGenerator::~HybridFaceNormalFluxGenerator() {
    delete faceMGraphs;
    delete cellMGraphs;
    delete pdeSolver;
    delete solutionFaceValues;
    delete solutionFaceFlux;
    delete solutionCellFlux;
}

VectorField HybridCellFluxGenerator::EvalSample(const cell &c) {
    RTLagrangeElement elem(*generator.pdeSolver->GetDisc(),
                           generator.solutionFaceFlux->U, 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, generator.solutionFaceFlux->U);
    }
    F *= (1 / area);
    return F;
}