Commit 86c7aec9 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactoring HybridFluxGenerator

parent a5fd0bb3
#ifndef HYBRIDFLUXGENERATOR_HPP
#define HYBRIDFLUXGENERATOR_HPP
#include "SampleGenerator.hpp"
#include "CirculantEmbedding.hpp"
#include <memory>
#include "pdesolver/PDESolver.hpp"
//class HybridFluxGenerator : public SampleGenerator<Scalar> {
//private:
// PDESolver *pdeSolver;
//
//public:
// HybridFluxGenerator(Meshes &meshes) : SampleGenerator(meshes) {
// std::string problemName = "StochasticLaplace2D";
// pdeSolver = PDESolverCreator(meshes).
// WithModel("HybridElliptic").
// WithProblem(problemName).Create();
// }
//};
//class MultilevelHybridFluxGenerator : public SampleGenerator {
//private:
// MatrixGraphs cellMGraphs;
// MatrixGraphs faceMGraphs;
//
// Vector *faceFlux = nullptr;
// Vector *faceValues = nullptr;
// Vector *cellFlux = nullptr;
//
//protected:
// IStochasticEllipticProblem *problem;
// HybridEllipticAssemble *assemble;
// RTLagrangeDiscretization *disc;
// Preconditioner *pc;
// Solver *solver;
// Newton *newton;
// Meshes &meshes;
//
// void drawSample(const SampleID &id) override {
// if (!id.coarse) {
// generateFineSample(id);
// if (plotting)
// plotter->PlotVector("flux", *cellFlux, 3,
// id.level.fine, "CellData");
// } else {
// generateCoarseSample(id);
// if (plotting)
// plotter->PlotVector("flux", *cellFlux, 3,
// id.level.coarse, "CellData");
// }
// }
//
//public:
// 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);
// }
// }
//
// virtual ~MultilevelHybridFluxGenerator() {
// delete problem;
// delete assemble;
// delete disc;
// delete solver;
// delete newton;
// }
//
// string Name() const override { return "Hybrid Flux Generator"; }
//
// void generateFineSample(SampleID id) {
// cellFlux = new Vector(cellMGraphs[id.level.mGfine]);
// faceValues = new Vector(faceMGraphs[id.level.mGfine]);
// faceFlux = new Vector(faceMGraphs[id.level.mGfine]);
//
// assemble->DrawSample(id);
// newton->operator()(*assemble, *faceValues);
// assemble->SetNormalFlux(*faceValues, *faceFlux);
// assemble->SetFlux(*faceValues, *cellFlux);
// }
//
// void generateCoarseSample(SampleID id) {
// cellFlux = new Vector(cellMGraphs[id.level.mGcoarse]);
// faceValues = new Vector(faceMGraphs[id.level.mGcoarse]);
// faceFlux = new Vector(faceMGraphs[id.level.mGcoarse]);
//
// assemble->DrawSample(id);
// newton->operator()(*assemble, *faceValues);
// assemble->SetNormalFlux(*faceValues, *faceFlux);
// assemble->SetFlux(*faceValues, *cellFlux);
// };
//
// VectorField EvalVectorFieldSample(const cell &c) override {
// RTLagrangeElement 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;
// }
//
// Scalar EvalSample(int face, const cell &c) override {
// RTLagrangeElement elem(*disc, *faceFlux, c);
// RTFaceElementT<> faceElem(*disc, *faceFlux, c, face);
// return (*faceFlux)(elem[face], 0) * elem.Sign(face) / faceElem.Area();
// }
//};
#endif //HYBRIDFLUXGENERATOR_HPP
#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;
}
#ifndef HYBRIDFLUXGENERATOR_HPP
#define HYBRIDFLUXGENERATOR_HPP
#include "generators/SampleGenerator.hpp"
#include "CirculantEmbedding.hpp"
class PDESolver;
class HybridFaceNormalFluxGenerator : public SampleGenerator<Scalar> {
private:
// Created with object
PlotMap plotMap;
// Created in constructor
MatrixGraphs *faceMGraphs;
MatrixGraphs *cellMGraphs;
void createPDESolver();
void createMGraphs();
void drawSample(const SampleID &id) override;
public:
// Created in constructor
PDESolver *pdeSolver;
// Created in drawSample
SampleSolution *solutionFaceValues;
SampleSolution *solutionFaceFlux;
SampleSolution *solutionCellFlux;
HybridFaceNormalFluxGenerator(Meshes &meshes) :
SampleGenerator(meshes) {
createPDESolver();
createMGraphs();
}
~HybridFaceNormalFluxGenerator();
string Name() const override {
return "HybridFaceNormalFluxGenerator";
}
Scalar EvalSample(int face, const cell &c) override;
};
class HybridCellFluxGenerator : public SampleGenerator<VectorField> {
const HybridFaceNormalFluxGenerator &generator;
public:
HybridCellFluxGenerator(Meshes &meshes,
const HybridFaceNormalFluxGenerator &generator) :
SampleGenerator(meshes),
generator(generator) {
}
VectorField EvalSample(const cell &c) override;
};
// MatrixGraphs cellMGraphs; cellMGraphs(MatrixGraphs(meshes, dof(new CellDoF(3)))),
// MatrixGraphs faceMGraphs; faceMGraphs(MatrixGraphs(meshes, dof(new FaceDoF(1)))),
//
// Vector *faceFlux = nullptr;
// Vector *faceValues = nullptr;
// Vector *cellFlux = nullptr;
#endif //HYBRIDFLUXGENERATOR_HPP
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment