MixedEllipticAssemble.cpp 14.1 KB
Newer Older
1
#include "MixedEllipticAssemble.hpp"
2
3


4
5
6
7
const char *MixedEllipticAssemble::Name() const {
    return "MixedEllipticAssemble";
}

8
9
10
void MixedEllipticAssemble::Initialize(Vector &u) const {
    u.ClearDirichletFlags();
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
11
12
        BFParts bnd(u.GetMesh(), c);
        if (!bnd.onBnd()) continue;
13
        for (int face = 0; face < c.Faces(); ++face) {
niklas.baumgarten's avatar
niklas.baumgarten committed
14
15
            RTFaceElementT<> faceElem(*disc, u, c, face);
            RowBndFaceValues uf_c(face, u, c, faceElem);
16
            if (bnd[face] == 2) {
niklas.baumgarten's avatar
niklas.baumgarten committed
17
                for (int j = 0; j < faceElem.size(); j++) {
18
                    double area = faceElem.Area();
niklas.baumgarten's avatar
niklas.baumgarten committed
19
                    VectorField F = problem->Flux(faceElem[0]());
20
                    VectorField N = faceElem.Normal();
niklas.baumgarten's avatar
niklas.baumgarten committed
21
22
                    uf_c(j) = area * (F * N);
                    uf_c.D(j) = true;
23
                }
24
            } else if (bnd[face] == 0) {
niklas.baumgarten's avatar
niklas.baumgarten committed
25
26
27
                for (int j = 0; j < faceElem.size(); j++) {
                    uf_c( j) = 0.0;
                    uf_c.D( j) = true;
28
                }
29
30
31
            }
        }
    }
32
    u.DirichletConsistent();
33
34
}

35
36
37
double MixedEllipticAssemble::Energy(const Vector &u) const {
    double energy = 0.0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
38
        RTLagrangeElementT<> elem(*disc, u, c);
39
40
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
41
            Tensor IK = Invert(problem->Permeability(c));
42
43
44
45
46
47
48
            VectorField Q = elem.VelocityField(q, u);
            energy += w * IK * Q * Q;
        }
    }
    return sqrt(PPM->Sum(energy));
}

49
void MixedEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
50
    RTLagrangeElementT<> elem(*disc, u, c);
51
    MixedRowBndValues r_c(r, c);
52
53
    for (int q = 0; q < elem.nQ(); ++q) {
        double w = elem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
54
        Scalar f = problem->Load(elem.QPoint(q));
55
        VectorField Q = elem.VelocityField(q, u);
niklas.baumgarten's avatar
niklas.baumgarten committed
56
        Tensor IK = Invert(problem->Permeability(c));
57
        Scalar divQ = elem.VelocityFieldDivergence(q, u);
58
        Scalar p = elem.PressureValue(q, u);
59
60
61
62
63
64
        for (int i = 0; i < elem.VelocitySize(); ++i) {
            for (int k = 0; k < elem.VelocityMaxk(i); ++k) {
                Velocity Q_i = elem.VelocityField(q, i, k);
                Scalar divQ_i = elem.VelocityFieldDivergence(q, i, k);
                r_c(elem.VelocityIndex(), i, k) += w * (((IK * Q) * Q_i) + divQ_i * p);
            }
65
        }
66
67
68
        for (unsigned int j = 0; j < elem.PressureSize(); ++j) {
            Scalar p_j = elem.PressureValue(q, j);
            r_c(elem.PressureIndex(), j, 0) += w * (divQ - f) * p_j;
69
70
71
        }
    }
    if (!r_c.onBnd()) return;
72
73
    for (int face = 0; face < c.Faces(); ++face) {
        if (r_c.bc(face) == 1) {
niklas.baumgarten's avatar
niklas.baumgarten committed
74
75
            RTFaceElementT<> faceElem(*disc, u, c, face);
            RowBndFaceValues rf_c(face, r, c, faceElem);
76
77
78
            for (int q = 0; q < faceElem.nQ(); ++q) {
                double w = faceElem.QWeight(q);
                Scalar P = problem->Solution(faceElem.QPoint(q));
niklas.baumgarten's avatar
niklas.baumgarten committed
79
80
81
82
                for (int i = 0; i < faceElem.size(); i++) {
                    Velocity face_Q_i = faceElem.VectorValue(q, i);
                    double uin = face_Q_i * faceElem.QNormal(q);
                    rf_c(0, i) -= w * uin * P;
83
                }
84
85
86
87
88
89
            }
        }
    }
}

void MixedEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const {
niklas.baumgarten's avatar
niklas.baumgarten committed
90
    RTLagrangeElementT<> elem(*disc, u, c);
91
92
93
    MixedRowEntries A_c(A, c, elem);
    for (int q = 0; q < elem.nQ(); ++q) {
        double w = elem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
94
        Tensor IK = Invert(problem->Permeability(c));
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
        for (int i = 0; i < elem.VelocitySize(); ++i) {
            for (int k = 0; k < elem.VelocityMaxk(i); ++k) {
                VectorField Q_i = elem.VelocityField(q, i, k);
                VectorField IKQ_i = IK * Q_i;
                Scalar divQ_i = elem.VelocityFieldDivergence(q, i, k);
                for (int j = 0; j < elem.VelocitySize(); ++j) {
                    for (int l = 0; l < elem.VelocityMaxk(j); ++l) {
                        Velocity Q_j = elem.VelocityField(q, j, l);
                        A_c(elem.VelocityIndex(), elem.VelocityIndex(), i, j, k, l) +=
                            w * (IKQ_i * Q_j);
                    }
                }
                for (unsigned int j = 0; j < elem.PressureSize(); ++j) {
                    Scalar p_j = elem.PressureValue(q, j);
                    A_c(elem.PressureIndex(), elem.VelocityIndex(), j, i, 0, k) +=
                        w * divQ_i * p_j;
                    A_c(elem.VelocityIndex(), elem.PressureIndex(), i, j, k, 0) +=
                        w * divQ_i * p_j;
                }
114
115
116
117
118
            }
        }
    }
}

119
120
121
double MixedEllipticAssemble::EnergyError(const Vector &u) const {
    double err = 0.0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
122
        RTLagrangeElementT<> elem(*disc, u, c);
123
124
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
125
            Tensor IK = Invert(problem->Permeability(c));
126
            VectorField Q = elem.VelocityField(q, u);
niklas.baumgarten's avatar
niklas.baumgarten committed
127
            VectorField F = problem->Flux(elem.QPoint(q));
128
129
            VectorField Diff = Q - F;
            err += w * IK * Diff * Diff;
130
131
132
133
134
        }
    }
    return sqrt(PPM->Sum(err));
}

135
136
137
double MixedEllipticAssemble::L2(const Vector &u) const {
    double l2 = 0.0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
138
        RTLagrangeElementT<> elem(*disc, u, c);
139
140
141
142
143
144
145
146
147
148
149
150
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
            Scalar p = elem.PressureValue(q, u);
            l2 += w * p * p;
        }
    }
    return sqrt(PPM->Sum(l2));
}

double MixedEllipticAssemble::H1(const Vector &u) const {
    double l2 = 0.0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
151
        RTLagrangeElementT<> elem(*disc, u, c);
152
153
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
154
            Tensor IK = Invert(problem->Permeability(c));
155
156
157
158
159
160
161
162
163
164
165
166
            VectorField Q = elem.VelocityField(q, u);
            Scalar p = elem.PressureValue(q, u);
            VectorField IK_Q = IK * Q;
            l2 += w * IK_Q * IK_Q + w * p * p;
        }
    }
    return sqrt(PPM->Sum(l2));
}

double MixedEllipticAssemble::L2Error(const Vector &u) const {
    double err = 0.0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
167
        RTLagrangeElementT<> elem(*disc, u, c);
168
169
170
171
172
173
174
175
176
177
178
179
180
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
            Scalar p = elem.PressureValue(q, u);
            Scalar Sol = problem->Solution(elem.QPoint(q));
            err += w * (p - Sol) * (p - Sol);
        }
    }
    return sqrt(PPM->Sum(err));
}

double MixedEllipticAssemble::L2CellAverageError(const Vector &u) const {
    double err = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
181
        RTLagrangeElementT<> elem(*disc, u, c);
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
        double w = 0;
        Scalar p = 0;
        Scalar Sol = 0;
        for (int q = 0; q < elem.nQ(); ++q) {
            w += elem.QWeight(q);
            p += elem.PressureValue(q, u);
            Sol += problem->Solution(elem.QPoint(q));
        }
        err += w * (p - Sol) * (p - Sol);
    }
    return sqrt(PPM->Sum(err));
}

double MixedEllipticAssemble::MaxError(const Vector &u) const {
    double err = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
198
        RTLagrangeElementT<> elem(*disc, u, c);
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
        for (int q = 0; q < elem.nQ(); ++q) {
            Scalar p = elem.PressureValue(q, u);
            err = max(err, abs(p - problem->Solution(elem.QPoint(q))));
        }
    }
    return PPM->Max(err);
}

double MixedEllipticAssemble::FluxError(const Vector &u) const {
    double flux_error = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        BFParts bnd(u.GetMesh(), c);
        if (!bnd.onBnd()) continue;
        for (int face = 0; face < c.Faces(); ++face) {
            if (bnd[face] == 2) {
niklas.baumgarten's avatar
niklas.baumgarten committed
214
                RTFaceElementT<> faceElem(*disc, u, c, face);
215
216
                for (int q = 0; q < faceElem.nQ(); ++q) {
                    double w = faceElem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
217
                    VectorField F = problem->Flux(faceElem.QPoint(q));
218
219
220
221
                    Scalar Diff = (faceElem.VectorValue(q, u) - F) * faceElem.QNormal(q);
                    flux_error += w * Diff * Diff;
                }
            } else if (bnd[face] == 0) {
niklas.baumgarten's avatar
niklas.baumgarten committed
222
                RTFaceElementT<> faceElem(*disc, u, c, face);
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
                for (int q = 0; q < faceElem.nQ(); ++q) {
                    double w = faceElem.QWeight(q);
                    Scalar FN = faceElem.VectorValue(q, u) * faceElem.QNormal(q);
                    flux_error += w * FN * FN;
                }
            }
        }
    }
    return sqrt(PPM->Sum(flux_error));
}

double MixedEllipticAssemble::FaceError(const Vector &u) const {
    double face_error = 0;  // For P0 without any pressure values on face
    return sqrt(PPM->Sum(face_error));
}

FluxPair MixedEllipticAssemble::InflowOutflow(const Vector &u) const {
    double inflow = 0;
    double outflow = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        BFParts bnd(u.GetMesh(), c);
        if (!bnd.onBnd()) continue;
        for (int face = 0; face < c.Faces(); ++face) {
            if (bnd[face] < 0) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
247
            RTFaceElementT<> faceElem(*disc, u, c, face);
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
            for (int q = 0; q < faceElem.nQ(); ++q) {
                double w = faceElem.QWeight(q);
                Scalar FN = faceElem.VectorValue(q, u) * faceElem.QNormal(q);
                if (FN > 0) { outflow += w * FN; }
                else { inflow += w * FN; }
            }
        }
    }
    return {PPM->Sum(inflow), PPM->Sum(outflow)};
}

FluxPair MixedEllipticAssemble::PrescribedInflowOutflow(const Vector &u) const {
    double inflow = 0;
    double outflow = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        BFParts bnd(u.GetMesh(), c);
        if (!bnd.onBnd()) continue;
        for (int face = 0; face < c.Faces(); ++face) {
            if (bnd[face] < 2) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
267
            RTFaceElementT<> faceElem(*disc, u, c, face);
268
269
            for (int q = 0; q < faceElem.nQ(); ++q) {
                double w = faceElem.QWeight(q);
niklas.baumgarten's avatar
niklas.baumgarten committed
270
                Scalar FN = problem->Flux(faceElem.QPoint(q)) *
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
                    faceElem.QNormal(q);
                if (FN > 0) { outflow += w * FN; }
                else { inflow += w * FN; }
            }
        }
    }
    return {PPM->Sum(inflow), PPM->Sum(outflow)};
}

FluxPair MixedEllipticAssemble::OutflowLeftRight(const Vector &u) const {
    double outflowLeft = 0;
    double outflowRight = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        BFParts bnd(u.GetMesh(), c);
        if (!bnd.onBnd()) continue;
        for (int face = 0; face < c.Faces(); ++face) {
            if (bnd[face] != 1) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
288
            RTFaceElementT<> faceElem(*disc, u, c, face);
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
            for (int q = 0; q < faceElem.nQ(); ++q) {
                double w = faceElem.QWeight(q);
                VectorField Q = faceElem.VectorValue(q, u);
                Scalar F = Q * faceElem.QNormal(q);
                if (F > 0 && c()[0] <= 1) { outflowLeft += w * F; }
                else if (F > 0 && c()[0] >= 2) outflowRight += w * F;
            }
        }
    }
    return {PPM->Sum(outflowLeft), PPM->Sum(outflowRight)};
}

double MixedEllipticAssemble::GoalFunctional(const Vector &u) const {
    double goal = 0;
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
        BFParts bnd(u.GetMesh(), c);
        if (!bnd.onBnd()) continue;
        for (int face = 0; face < c.Faces(); ++face) {
            if (bnd[face] != 1) continue;
niklas.baumgarten's avatar
niklas.baumgarten committed
308
            RTFaceElementT<> faceElem(*disc, u, c, face);
309
310
311
312
313
314
315
316
317
318
319
320
321
322
            for (int q = 0; q < faceElem.nQ(); ++q) {
                Point z = faceElem.QPoint(q);
                if (z[0] < 0.25) continue;
                if (z[0] > 0.5) continue;
                double w = faceElem.QWeight(q);
                VectorField Q = faceElem.VectorValue(q, u);
                Scalar F = Q * faceElem.QNormal(q);
                goal += w * F;
            }
        }
    }
    return PPM->Sum(goal);
}

323
324
void MixedEllipticAssemble::SetExactSolution(Vector &u) const {
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
325
        RTLagrangeElementT<> elem(*disc, u, c);
326
327
328
329
330
331
332
        vector<Point> nodalPoints;
        u.dof::NodalPoints(elem.PressureIndex(), c, nodalPoints);
        MixedRowValues u_c(u, c, elem);
        for (int i = 0; i < elem.PressureSize(); i++) {
            u_c(elem.PressureIndex(), i) = problem->Solution(nodalPoints[i]);
        }
        for (int face = 0; face < c.Faces(); ++face) {
niklas.baumgarten's avatar
niklas.baumgarten committed
333
334
            RTFaceElementT<> faceElem(*disc, u, c, face);
            RowFaceValues uf_c(face, u, c, faceElem);
335
336
            VectorField N = faceElem.Normal();
            vector<Point> nodalPointsOnFace;
niklas.baumgarten's avatar
niklas.baumgarten committed
337
338
339
340
            u.NodalPointsOnFace(0, c, face, nodalPointsOnFace);
            for (int j = 0; j < faceElem.size(); j++) {
                VectorField F = problem->Flux(c());
                uf_c( j) = faceElem.Area() * F * N;
341
342
343
344
345
346
347
            }
        }
    }
}

void MixedEllipticAssemble::SetFlux(const Vector &u, Vector &flux) {
    for (cell c = u.cells(); c != u.cells_end(); ++c) {
niklas.baumgarten's avatar
niklas.baumgarten committed
348
        RTLagrangeElementT<> elem(*disc, u, c);
349
350
351
352
353
354
355
356
357
        VectorField F = zero;
        double area = 0;
        for (int q = 0; q < elem.nQ(); ++q) {
            double w = elem.QWeight(q);
            VectorField Q = elem.VelocityField(q, u);
            F += w * Q;
            area += w;
        }
        F *= (1 / area);
niklas.baumgarten's avatar
niklas.baumgarten committed
358
        for (int d = 0; d < SpaceDimension; ++d) {
359
            flux(c(), d) = F[d];
360
361
362
363
        }
    }
}

364
void MixedEllipticAssemble::AssembleTransfer(TransferMatrix &TM) const {
365
    // Only works for rt0_p0
366
367
368
369
370
371
372
373
374
375
376
377
378
379
    TM = 0;
    const matrixgraph &cg = TM.CoarseMatrixGraph();
    for (cell c = cg.cells(); c != cg.cells_end(); ++c) {
        for (int cf = 0; cf < c.Faces(); ++cf) {
            for (int k = 0; k < 2; ++k)
                TM(c.Face(cf), 0.5 * (c.Face(cf) + c.FaceCorner(cf, k)))[0] = 0.5;
            if (c.Faces() == 3) {
                for (int f = cf + 1; f < c.Faces(); ++f)
                    TM(c.Face(cf), 0.5 * (c.Face(cf) + c.Face(f)))[0] = 0.5;
            } else
                TM(c.Face(cf), 0.5 * (c.Face(cf) + c()))[0] = 0.5;
        }
    }
}