Commit 266a8c90 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

adapted CirculantEmbedding

parent 04e29167
...@@ -53,20 +53,16 @@ void CirculantEmbedding1D::generateCoarseSample(const SampleID &id, ...@@ -53,20 +53,16 @@ void CirculantEmbedding1D::generateCoarseSample(const SampleID &id,
} }
ComplexField1D CirculantEmbedding1D::generateField(const SampleID &id) { ComplexField1D CirculantEmbedding1D::generateField(const SampleID &id) {
generator.DrawSample(id);
ComplexField1D normalDist = generator.EvalComplexField1DSample();
ComplexField1D Xf; ComplexField1D Xf;
ComplexField1D normallyDistributed;
normallyDistributed.resize(numberOfCellsEmbedded);
fineComplexField.clear(); fineComplexField.clear();
Xf.resize(numberOfCellsEmbedded); Xf.resize(numberOfCellsEmbedded);
for (int i = 0; i < numberOfCellsEmbedded; i++) {
normalDist.DrawSample(id);
normallyDistributed[i] = normalDist.EvalComplexSample();
}
ComplexField1D product; ComplexField1D product;
product.resize(numberOfCellsEmbedded); product.resize(numberOfCellsEmbedded);
for (int i = 0; i < numberOfCellsEmbedded; i++) { for (int i = 0; i < numberOfCellsEmbedded; i++) {
product[i] = sqrtEigenvalues[i] * normallyDistributed[i]; product[i] = sqrtEigenvalues[i] * normalDist[i];
} }
fft(product, Xf); fft(product, Xf);
double divisor = sqrt(numberOfCellsEmbedded); double divisor = sqrt(numberOfCellsEmbedded);
......
...@@ -5,10 +5,9 @@ ...@@ -5,10 +5,9 @@
#include "basics/Utilities.hpp" #include "basics/Utilities.hpp"
#include "Algebra.hpp" #include "Algebra.hpp"
#include "SampleGenerator.hpp" #include "SampleGenerator.hpp"
#include "dof/BasicDoFs.hpp"
#include "NormalDistribution.hpp" #include "NormalDistribution.hpp"
#include "basics/PlotMap.hpp"
typedef std::vector<std::complex<double>> ComplexField1D;
typedef std::vector<double> EigenValues1D; typedef std::vector<double> EigenValues1D;
...@@ -20,8 +19,6 @@ typedef std::vector<double> ToeplitzColumn; ...@@ -20,8 +19,6 @@ typedef std::vector<double> ToeplitzColumn;
typedef std::vector<double> CirculantRow; typedef std::vector<double> CirculantRow;
typedef std::vector<std::vector<std::complex<double>>> ComplexField2D;
typedef std::vector<std::vector<double>> EigenValues2D; typedef std::vector<std::vector<double>> EigenValues2D;
typedef std::vector<std::vector<double>> SqrtEigenValues2D; typedef std::vector<std::vector<double>> SqrtEigenValues2D;
...@@ -52,7 +49,7 @@ public: ...@@ -52,7 +49,7 @@ public:
explicit CirculantEmbedding(Meshes &meshes) : explicit CirculantEmbedding(Meshes &meshes) :
meshes(meshes), meshes(meshes),
normalDist(meshes), normalDist(meshes),
cellMGraphs(meshes, dof(new CellDoF(1))) { cellMGraphs(meshes, dof(new LagrangeDoF(0))) {
config.get("evtol", evtol); config.get("evtol", evtol);
config.get("StochasticField", fieldType); config.get("StochasticField", fieldType);
...@@ -72,10 +69,13 @@ public: ...@@ -72,10 +69,13 @@ public:
class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D { class CirculantEmbedding1D : public CirculantEmbedding, public CovarianceFunction1D {
public: public:
NormalDistributionComplex1DField generator;
ComplexField1D fineComplexField; ComplexField1D fineComplexField;
SqrtEigenValues1D sqrtEigenvalues; SqrtEigenValues1D sqrtEigenvalues;
int numberOfCellsEmbedded = 0; int numberOfCellsEmbedded = 0;
int numberOfCells = 0; int numberOfCells = 0;
...@@ -83,10 +83,12 @@ public: ...@@ -83,10 +83,12 @@ public:
double domain_end = 1; double domain_end = 1;
explicit CirculantEmbedding1D(int l, Meshes &meshes) : explicit CirculantEmbedding1D(int l, Meshes &meshes) :
CirculantEmbedding(meshes) { CirculantEmbedding(meshes),
generator(meshes) {
numberOfCells = meshes[l].CellCount(); numberOfCells = meshes[l].CellCount();
sqrtEigenvalues = computeSqrtEV(); sqrtEigenvalues = computeSqrtEV();
numberOfCellsEmbedded = sqrtEigenvalues.size(); numberOfCellsEmbedded = sqrtEigenvalues.size();
generator.Resize(numberOfCellsEmbedded);
} }
void generateFineSample(const SampleID &id, void generateFineSample(const SampleID &id,
...@@ -168,6 +170,8 @@ public: ...@@ -168,6 +170,8 @@ public:
class MultilevelCirculantEmbedding : public SampleGenerator { class MultilevelCirculantEmbedding : public SampleGenerator {
private: private:
PlotMap plotMap;
std::map<int, CirculantEmbedding *> circulantEmbeddings; std::map<int, CirculantEmbedding *> circulantEmbeddings;
Vector *fineSample = nullptr; Vector *fineSample = nullptr;
...@@ -179,15 +183,16 @@ protected: ...@@ -179,15 +183,16 @@ protected:
if (!id.coarse) { if (!id.coarse) {
circulantEmbeddings[id.level.fine]->generateFineSample(id, fineSample, circulantEmbeddings[id.level.fine]->generateFineSample(id, fineSample,
coarseSample); coarseSample);
// if (plotting) SampleSolution _fineSample(*fineSample, id); // Todo
// plotter->PlotVector("kappa", *fineSample, 1, _fineSample.U = 1.0;
// id.level.coarse, "CellData"); _fineSample.name = "Kappa";
plotMap.VtkPlot(_fineSample);
} else { } else {
circulantEmbeddings[id.level.fine]->generateCoarseSample(id, fineSample, circulantEmbeddings[id.level.fine]->generateCoarseSample(id, fineSample,
coarseSample); coarseSample);
// if (plotting) SampleSolution _coarseSample(*coarseSample, id); // Todo
// plotter->PlotVector("kappa", *coarseSample, 1, _coarseSample.name = "Kappa";
// id.level.fine, "CellData"); plotMap.VtkPlot(_coarseSample);
} }
} }
...@@ -216,6 +221,4 @@ public: ...@@ -216,6 +221,4 @@ public:
string Name() const override { return "Circulant Embedding"; } string Name() const override { return "Circulant Embedding"; }
}; };
#endif //M_CIRCULANTEMBEDDING_H #endif //M_CIRCULANTEMBEDDING_H
\ No newline at end of file
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