StochasticFields.C 4.92 KB
Newer Older
1
2
#include "StochasticFields.h"

niklas.baumgarten's avatar
niklas.baumgarten committed
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
62
63
64
65
66
67
68
69
class StochasticFieldPair1D : public StochasticFieldPair {
    int N;
    CirculantEmbedding1D *circEmbedding1D;
    vector<double> fieldf;
    vector<double> fieldc;
public:
    explicit StochasticFieldPair1D(int l) : StochasticFieldPair() {
        // Todo rmv Hack!
        N = (int) pow(2, l);
        circEmbedding1D = new CirculantEmbedding1D(N);
    }

    ~StochasticFieldPair1D() { delete circEmbedding1D; }

    void generate() override {
        if (PPM->master()) {
            PPM->BroadcastInt(N);
            fieldf = circEmbedding1D->getFieldf();
            for (double &cell_value : fieldf) {
                PPM->BroadcastDouble(cell_value);
            }
            fieldc = circEmbedding1D->generateInterpolField(fieldf);
            for (double &cell_value : fieldc) {
                PPM->BroadcastDouble(cell_value);
            }
        } else {
            N = PPM->BroadcastInt();
            fieldf.resize(N);
            for (double &cell_value : fieldf)
                cell_value = PPM->BroadcastDouble();
            fieldc.resize(N / 2);
            for (double &cell_value : fieldc)
                cell_value = PPM->BroadcastDouble();
        }
    }

    double getFineValue(const Point &z) const override {
        int i = floor(z[0] * N);
        return fieldf[i];
    }

    double getCoarseValue(const Point &z) const override {
        int i = floor(z[0] * N / 2);
        return fieldc[i];
    }

    pair<double, double> getFineMinMaxValues() const override {
        double max = 0.0;
        double min = infty;
        for (int i = 0; i < N; i++) {
            if (fieldf[i] > max) max = fieldf[i];
            if (fieldf[i] < min) min = fieldf[i];
        }
        return {min, max};
    }

    pair<double, double> getCoarseMinMaxValues() const override {
        double max = 0.0;
        double min = infty;
        for (int i = 0; i < N; i++) {
            if (fieldc[i] > max) max = fieldc[i];
            if (fieldc[i] < min) min = fieldc[i];
        }
        return {min, max};
    }
};

70
71
72
73
class StochasticFieldPair2D : public StochasticFieldPair {
    int N1;
    int N2;
    CirculantEmbedding2D *circEmbedding2D;
niklas.baumgarten's avatar
niklas.baumgarten committed
74
75
    std::vector<std::vector<double>> fieldf;
    std::vector<std::vector<double>> fieldc;
76
77
78
79
80
81
82
83
84
85
86
87
88
89
public:
    explicit StochasticFieldPair2D(int l) : StochasticFieldPair() {
        // Todo rmv Hack!
        N1 = (int) pow(2, l);
        N2 = N1;
        circEmbedding2D = new CirculantEmbedding2D(N1, N2);
    }

    ~StochasticFieldPair2D() { delete circEmbedding2D; }

    void generate() override {
        if (PPM->master()) {
            PPM->BroadcastInt(N1);
            PPM->BroadcastInt(N2);
niklas.baumgarten's avatar
niklas.baumgarten committed
90
            fieldf = circEmbedding2D->getFieldf();
91
92
            for (int i = 0; i < N2; i++)
                for (int j = 0; j < N1; j++)
niklas.baumgarten's avatar
niklas.baumgarten committed
93
94
                    PPM->BroadcastDouble(fieldf[i][j]);
            fieldc = circEmbedding2D->generateInterpolField(fieldf);
95
96
            for (int i = 0; i < N2 / 2; i++)
                for (int j = 0; j < N1 / 2; j++)
niklas.baumgarten's avatar
niklas.baumgarten committed
97
                    PPM->BroadcastDouble(fieldc[i][j]);
98
99
100
        } else {
            N1 = PPM->BroadcastInt();
            N2 = PPM->BroadcastInt();
niklas.baumgarten's avatar
niklas.baumgarten committed
101
            fieldf.resize(N2);
102
            for (int i = 0; i < N2; i++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
103
                fieldf[i].resize(N1);
104
                for (int j = 0; j < N1; j++)
niklas.baumgarten's avatar
niklas.baumgarten committed
105
                    fieldf[i][j] = PPM->BroadcastDouble();
106
            }
niklas.baumgarten's avatar
niklas.baumgarten committed
107
            fieldc.resize(N2 / 2);
108
            for (int i = 0; i < N2 / 2; i++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
109
                fieldc[i].resize(N1 / 2);
110
                for (int j = 0; j < N1 / 2; j++)
niklas.baumgarten's avatar
niklas.baumgarten committed
111
                    fieldc[i][j] = PPM->BroadcastDouble();
112
113
114
115
116
117
118
            }
        }
    }

    double getFineValue(const Point &z) const override {
        int i = floor(z[0] * N2);
        int j = floor(z[1] * N1);
niklas.baumgarten's avatar
niklas.baumgarten committed
119
        return fieldf[i][j];
120
121
122
123
124
    }

    double getCoarseValue(const Point &z) const override {
        int i = floor(z[0] * N2 / 2);
        int j = floor(z[1] * N1 / 2);
niklas.baumgarten's avatar
niklas.baumgarten committed
125
        return fieldc[i][j];
126
127
128
129
130
131
132
    }

    pair<double, double> getFineMinMaxValues() const override {
        double max = 0.0;
        double min = infty;
        for (int i = 0; i < N2; i++) {
            for (int j = 0; j < N1; j++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
133
134
                if (fieldf[i][j] > max) max = fieldf[i][j];
                if (fieldf[i][j] < min) min = fieldf[i][j];
135
136
137
138
139
140
141
142
143
144
            }
        }
        return {min, max};
    }

    pair<double, double> getCoarseMinMaxValues() const override {
        double max = 0.0;
        double min = infty;
        for (int i = 0; i < N2 / 2; i++) {
            for (int j = 0; j < N1 / 2; j++) {
niklas.baumgarten's avatar
niklas.baumgarten committed
145
146
                if (fieldc[i][j] > max) max = fieldc[i][j];
                if (fieldc[i][j] < min) min = fieldc[i][j];
147
148
149
150
151
152
            }
        }
        return {min, max};
    }
};

niklas.baumgarten's avatar
niklas.baumgarten committed
153
154
155
156
StochasticFieldPair *getStochasticFieldPair(int l, int dim, string fieldName) {
    if (dim == 1)
        return new StochasticFieldPair1D(l);
    else if (dim == 2)
157
        return new StochasticFieldPair2D(l);
niklas.baumgarten's avatar
niklas.baumgarten committed
158
159
    else
        return nullptr;
160
}