Commit 289e9914 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

rmaining files commited

parent defee528
......@@ -3,8 +3,8 @@ Model = LinearFEM
#Model = MixedFEM
#Model = HybridFEM
#Model = FVFEM # Don't forget to swith on Overlap, sign and penalty
#Model = LinearDGFEM # Don't forget to swith on Overlaü, sign and penalty
#Model = QuadraticDGFEM # Don't forget to swith on Overlaü, sign and penalty
#Model = LinearDGFEM # Don't forget to swith on Overlap, sign and penalty
#Model = QuadraticDGFEM # Don't forget to swith on Overlap, sign and penalty
#Overlap = dG1
#sign = 1
......@@ -18,23 +18,23 @@ Problem = StochasticLaplace2D
StochasticField = LogNormal
#StochasticField = Gaussian
#Experiment = ConvergenceTest
Experiment = ConvergenceTest
#Experiment = MLMCExperiment
Experiment = MLMCOverEpsilon
#Experiment = MLMCOverEpsilon
#size_init = 7
#l_init = 3 4 5 6 7 8 9
#M_init = 500 500 500 500 500 500 500
size_init = 4
l_init = 3 4 5 6
M_init = 3 3 3 3
size_init = 3
l_init = 3 4 5
M_init = 12 6 3
M_init = 45 15 5
epsilon = 0.01
size_lst = 4
epsilon_lst = 0.01 0.005 0.001 0.0005
size_lst = 2
epsilon_lst = 0.05 0.01
Lmax = 11
Lmax = 9
plevel = 2
functional = L2
#functional = Energy
......@@ -45,7 +45,7 @@ functional = L2
MCVerbose = 1
MLMCVerbose = 2
MCPlotting = 3
MLMCPlotting = 0
MLMCPlotting = 2
enablePython = 1
......@@ -55,7 +55,5 @@ sigma = 1.0
norm_p = 2
lambda1 = 0.15
lambda2 = 0.15
alpha = 1.0
evtol = 1e-10
logfile = ref
alpha = 1.8
evtol = 1e-10
\ No newline at end of file
#include "DGDiscretization.h"
int DGDoF::get_cell_deg(const Point &z, int ud) const {
return cell_std_deg;
}
void DGDoF::set_cell_deg(const Point &z, int ud, short deg) {
vector<short> *c_deg = nullptr;
if (cell_deg.find(z) != cell_deg.end())
c_deg = &cell_deg[z];
else {
c_deg = &cell_deg[z];
c_deg->resize(1, -1);
}
(*c_deg)[ud] = deg;
}
int DGDoF::CellDoFs(const cell &c) const {
switch (c.dim()) {
case 1:
return get_cell_deg(c(), 0) + 1;
case 2:
switch (c.Corners()) {
case 3: {
int result = 0;
for (int i = get_cell_deg(c(), 0) + 1; i != 0; i--)
result += i;
return result;
}
case 4:
return (int) power(get_cell_deg(c(), 0) + 1, 2);
}
case 3:
switch (c.Corners()) {
case 4:
switch (get_cell_deg(c(), 0)) {
case 0:
return 1;
case 1:
return 4;
case 2:
return 10;
case 3:
return 20;
}
case 8:
return (int) power(get_cell_deg(c(), 0) + 1, 3);
}
}
Exit("not implemented");
}
void DGDoF::NodalDoFs(const cell &c, vector<short> &z) const {
z.resize(1);
z[0] = CellDoFs(c);
}
void DGDoF::NodalPoints(const cell &c, vector <Point> &z) const {
z.resize(1);
z[0] = c();
}
int DGDoF::NodalPoints(const cell &c) const { return 1; }
DGDiscretization::DGDiscretization(DGDoF *dof) : Discretization(dof, "DG"), D(*dof) {
S_int[0] = Get_Shape("P0int", "Qint1", "Qint1");
S_int[1] = Get_Shape("P1int", "Qint2", "Qint1");
S_int[2] = Get_Shape("P2int", "Qint3", "Qint1");
S_int[3] = Get_Shape("P3int", "Qint4", "Qint1");
S_int[4] = Get_Shape("P4int", "Qint5", "Qint1");
S_tri[0] = Get_Shape("P0tri", "Qtri1", "Qint1");
S_tri[1] = Get_Shape("P1tri", "Qtri7", "Qint2");
S_tri[2] = Get_Shape("P2tri", "Qtri7", "Qint3");
S_tri[3] = Get_Shape("P3tri", "Qtri16", "Qint4");
S_tri[4] = Get_Shape("P4tri", "Qtri25", "Qint5");
T_tri[1] = Get_Shape("RT0tri");
S_quad[0] = Get_Shape("P0quad", "Qquad1", "Qint1");
S_quad[1] = Get_Shape("P1quad", "Qquad4", "Qint2");
S_quad[2] = Get_Shape("P2quad", "Qquad9", "Qint3");
S_quad[3] = Get_Shape("P3quad", "Qquad16", "Qint4");
S_quad[4] = Get_Shape("P4quad", "Qquad25", "Qint5");
S_quad[5] = Get_Shape("P5quad", "Qquad25", "Qint5");
S_tet[1] = Get_Shape("P1tet");
S_tet[2] = Get_Shape("P2tet");
S_hex[1] = Get_Shape("P1hex");
S_hex[2] = Get_Shape("P2hex");
S_hex[3] = Get_Shape("P3hex");
T_quad[1] = Get_Shape("RT0quad");
T_tet[1] = Get_Shape("RT0tet");
T_hex[1] = Get_Shape("RT0hex");
L_P[0] = Get_Shape("P0int", "Qint1", "Qint1");
L_P[1] = Get_Shape("P1int", "Qint2", "Qint1");
L_P[2] = Get_Shape("P2int", "Qint3", "Qint1");
L_P[3] = Get_Shape("P3int", "Qint4", "Qint1");
L_P[4] = Get_Shape("P4int", "Qint5", "Qint1");
}
int DGDiscretization::getCellDegree(const Point &z, int component) const {
return D.get_cell_deg(z, component);
}
const Shape &DGDiscretization::getCellShape(const cell &c, int p, int component) const {
if (p > 5) Exit("GetShape: Not implemented");
switch (c.Corners()) {
case 2:
return *S_int[p]; // *L_P[p];
case 3:
return *S_tri[p];
case 4:
switch (c.dim()) {
case 2:
return *S_quad[p];
case 3:
return *S_tet[p];
}
case 8:
return *S_hex[p];
}
Exit("GetShape: Not implemented!");
}
const Quadrature &DGDiscretization::GetQuad(const cell &c, int p, int component) const {
const Shape &S = getCellShape(c, p, component);
return S.GetQuad();
}
const Quadrature &
DGDiscretization::GetFaceQuad(const cell &c, int p, int component) const {
const Shape &S = getCellShape(c, p, component);
return S.GetFaceQuad();
}
int DGDiscretization::CellDoFs(const cell &c) const { return D.CellDoFs(c); }
string DGDiscretization::DiscName() const { return "DGDiscretization"; }
Point DGElement::GlobalToLocal(const Point &z) const {
Tensor T;
T[0][0] = C[1][0] - C[0][0];
T[0][1] = C[2][0] - C[0][0];
T[1][0] = C[1][1] - C[0][1];
T[1][1] = C[2][1] - C[0][1];
if (dim == 2) {
T[0][2] = 0;
T[1][2] = 0;
T[2][0] = 0;
T[2][1] = 0;
T[2][2] = 1;
} else {
T[0][2] = C[3][0] - C[0][0];
T[1][2] = C[3][1] - C[0][1];
T[2][0] = C[1][2] - C[0][2];
T[2][1] = C[2][2] - C[0][2];
T[2][2] = C[3][2] - C[0][2];
}
Tensor IT = Invert(T);
Point d = z - C[0];
Point v;
v[0] = real(IT[0][0] * d[0] + IT[0][1] * d[1] + IT[0][2] * d[2]);
v[1] = real(IT[1][0] * d[0] + IT[1][1] * d[1] + IT[1][2] * d[2]);
v[2] = real(IT[2][0] * d[0] + IT[2][1] * d[1] + IT[2][2] * d[2]);
return v;
}
int DGElement::dimension() const { return S.size(); }
int DGElement::get_deg() const { return deg; }
DGElement::DGElement(const DGDiscretization &disc, const matrixgraph &g, const cell &c)
: Element(g, c, disc.GetQuad(c, disc.getCellDegree(c()))),
deg(disc.getCellDegree(c())),
dim(c.dim()),
S(disc.getCellShape(c, deg)),
value(S.values()), C(c) {
for (int q = 0; q < Q.size(); ++q) {
T[q] = c.GetTransformation(Q.QPoint(q));
qWeight[q] = T[q].Det() * Q.Weight(q);
qPoint[q] = c[Q.QPoint(q)];
for (int i = 0; i < S.size(); ++i)
gradient[q][i] = T[q] * S.LocalGradient(q, i);
}
}
void DGElement::NodalPoints(Point *z) const { S.NodalPoints(C, z); }
int DGElement::nQ() const { return Q.size(); }
double DGElement::QWeight(int q) const { return qWeight[q]; }
const Point &DGElement::QPoint(int q) const { return qPoint[q]; }
double DGElement::Area() const {
double a = 0;
for (int q = 0; q < nQ(); ++q)
a += qWeight[q];
return a;
}
double DGElement::Value(int q, int i) const { return value[q][i]; }
double DGElement::Value(int q, const Vector &u) const {
double p = 0;
for (int i = 0; i < S.size(); ++i)
p += u(r(0), i) * value[q][i];
return p;
}
VectorField DGElement::Derivative(int q, int i) const { return gradient[q][i]; }
VectorField DGElement::Derivative(int q, const Vector &u) const {
VectorField Du = zero;
for (int i = 0; i < S.size(); ++i)
Du += u(r(0), i) * gradient[q][i];
return Du;
}
const Point &DGFaceElement::LocalFaceCorner(int i) const {
return c.LocalCorner(c.facecorner(fid, i));
}
Point DGFaceElement::Local(const Point &local) const {
switch (c.FaceCorners(fid)) {
case 1:
return LocalFaceCorner(0);
case 2:
return (1 - local[0]) * LocalFaceCorner(0)
+ local[0] * LocalFaceCorner(1);
case 3:
return (1 - local[0] - local[1]) * LocalFaceCorner(0)
+ local[0] * LocalFaceCorner(1)
+ local[1] * LocalFaceCorner(2);
case 4:
return (1 - local[0]) * (1 - local[1]) * LocalFaceCorner(0)
+ local[0] * (1 - local[1]) * LocalFaceCorner(1)
+ local[0] * local[1] * LocalFaceCorner(2)
+ (1 - local[0]) * local[1] * LocalFaceCorner(3);
default: Exit("Not implemented!")
return local;
}
}
DGFaceElement::DGFaceElement(const DGDiscretization &D, const matrixgraph &g,
const cell &C, int f_id)
: FaceElement(g, C, f_id, D.GetFaceQuad(C, D.getCellDegree(C()))),
S(D.getCellShape(C, D.getCellDegree(C()))),
c(C), f(g.GetMesh().find_face(C.Face(f_id))),
fid(f_id), dim(C.dim()), Q(S.GetFaceQuad()) {
d = dim * S.size();
ww = 0;
for (int q = 0; q < Q.size(); ++q) ww += Q.Weight(q);
ww = c.LocalFaceArea(fid) / ww;
for (int q = 0; q < Q.size(); ++q) {
const Point &x = Q.QPoint(q);
Point y = Local(x);
qLocal[q] = y;
qPoint[q] = c[y];
T[q] = c.GetTransformation(y);
qNormal[q] = T[q] * c.LocalFaceNormal(fid);
double w = norm(qNormal[q]);
qWeight[q] = T[q].Det() * w * ww * Q.Weight(q);
qNormal[q] *= (1 / w);
}
for (int q = 0; q < Q.size(); ++q) {
for (int i = 0; i < S.size(); ++i) {
value[q][i] = S(QLocal(q), i);
gradient[q][i] = T[q] * S.LocalGradient(QLocal(q), i);
}
}
}
int DGFaceElement::dimension() const { return S.size(); }
double DGFaceElement::QWeight(const Point &y) const {
Transformation TR = c.GetTransformation(y);
Point N = TR * c.LocalFaceNormal(fid);
return norm(N);
}
double DGFaceElement::Value(int q, const Vector &u) const {
double p = 0;
for (int i = 0; i < S.size(); ++i)
p += u(c(), i) * value[q][i];
return p;
}
VectorField DGFaceElement::Derivative(int q, int i) const { return gradient[q][i]; }
double DGFaceElement::Value(int q, int i) const { return value[q][i]; }
int DGFaceElement::nQ() const { return Q.size(); }
double DGFaceElement::QWeight(int q) const { return qWeight[q]; }
const Point &DGFaceElement::QLocal(int q) const { return qLocal[q]; }
const Point &DGFaceElement::QPoint(int q) const { return qPoint[q]; }
const Point &DGFaceElement::QNormal(int q) const { return qNormal[q]; }
VectorField DGFaceElement::Derivative(int q, const Vector &u) const {
VectorField Du = zero;
for (int i = 0; i < S.size(); ++i)
Du += u(c(), i) * gradient[q][i];
return Du;
}
#ifndef TUTORIAL_DGDISCRETIZATION_H
#define TUTORIAL_DGDISCRETIZATION_H
#include "m++.h"
#include "DGMatrixGraph.h"
class DGDoF : public DoF {
int cell_std_deg;
hash_map<Point, vector < short>, Hash>
cell_deg;
public:
explicit DGDoF(int cellDegree) : DoF(0, false), cell_std_deg(cellDegree) {}
int get_cell_deg(const Point &z, int ud = 0) const;
void set_cell_deg(const Point &z, int ud, short deg);
int CellDoFs(const cell &c) const;
void NodalDoFs(const cell &c, vector<short> &z) const override;
int NodalPoints(const cell &c) const override;
void NodalPoints(const cell &c, vector <Point> &z) const override;
string Name() const override { return "DGDoF"; }
int TypeDoF(int) const override { return CELL; }
};
class DGDiscretization : public Discretization {
const Shape *L_P[5]{};
const Shape *S_int[5]{};
const Shape *S_tri[5]{};
const Shape *S_quad[6]{};
const Shape *S_tet[4]{};
const Shape *S_hex[4]{};
const Shape *T_tri[5]{};
const Shape *T_quad[5]{};
const Shape *T_tet[2]{};
const Shape *T_hex[2]{};
DGDoF &D;
public:
explicit DGDiscretization(DGDoF *dof);
int getCellDegree(const Point &z, int component = 0) const;
const Shape &getCellShape(const cell &c, int p, int component = 0) const;
const Quadrature &GetQuad(const cell &c, int p, int component = 0) const;
const Quadrature &GetFaceQuad(const cell &c, int p, int component = 0) const;
int CellDoFs(const cell &c) const;
string DiscName() const override;
};
class DGElement : public Element {
int deg;
int dim;
const Shape &S;
const ShapeValue &value;
ShapeGradient gradient;
const cell &C;
public:
Point GlobalToLocal(const Point &z) const;
int dimension() const;
int get_deg() const;
DGElement(const DGDiscretization &disc,
const matrixgraph &g, const cell &c);
void NodalPoints(Point *z) const;
int nQ() const override;
double QWeight(int q) const override;
const Point &QPoint(int q) const override;
double Area() const override;
double Value(int q, int i) const;
double Value(int q, const Vector &u) const;
VectorField Derivative(int q, int i) const;
VectorField Derivative(int q, const Vector &u) const;
};
class DGFaceElement : public FaceElement {
int d;
int dim;
const Shape &S;
const cell &c;
face f;
int fid;
const Quadrature &Q;
Point qPoint[MaxQuadraturePoints];
Point qLocal[MaxQuadraturePoints];
Point qNormal[MaxQuadraturePoints];
double qWeight[MaxQuadraturePoints];
Transformation T[MaxQuadraturePoints];
ShapeValue value;
ShapeGradient gradient;
double ww;
public:
const Point &LocalFaceCorner(int i) const;
Point Local(const Point &local) const;
DGFaceElement(const DGDiscretization &D, const matrixgraph &g,
const cell &C, int f_id);
int dimension() const;
double QWeight(const Point &y) const;
const Point &QLocal(int q) const override;
const Point &QPoint(int q) const override;
const Point &QNormal(int q) const override;
double QWeight(int q) const override;
int nQ() const override;
double Value(int q, int i) const;
double Value(int q, const Vector &u) const;
VectorField Derivative(int q, int i) const;
VectorField Derivative(int q, const Vector &u) const;
};
#endif //TUTORIAL_DGDISCRETIZATION_H
#include "DGMatrixGraph.h"
DGRowEntries::DGRowEntries(Matrix &A, const rows &R, const rows &Rf) {
n = R[0].n();
nf = Rf[0].n();
a = A(R[0], Rf[0]);
}
DGRowEntries::DGRowEntries(Matrix &A, const cell &c, const cell &cf) {
row r = A.find_row(c());
row rf = A.find_row(cf());
n = r.n();
nf = rf.n();
a = A(r, rf);
}
DGRowEntries::DGRowEntries(Matrix &A, const cell &c) {
row r = A.find_row(c());
n = nf = r.n();
a = A(r, r);
}
Scalar &DGRowEntries::operator()(int k, int l) {
if (k * nf + l >= n * nf) {
mout << OUT (n) << OUT (nf) << OUT (k) << OUT (l) << endl;
Exit ("too large");
}
return a[k * nf + l];
}
void EGRowEntries::fill(Matrix &A, const rows &R, const rows &R2) {
for (unsigned int i = 0; i < R.size(); ++i) {
n[i] = R[i].n();
for (unsigned int j = 0; j < R2.size(); ++j) {
if ((i == 3) || (j == 3))
a[i][j] = A(R[i], R2[j]);
}
}
}
EGRowEntries::EGRowEntries(Matrix &A, const rows &R, const rows &R2)
: a{}, n{} { fill(A, R, R2); }
EGRowEntries::EGRowEntries(Matrix &A, const cell &c, const cell &c2)
: a{}, n{} {
rows R(A, c);
rows R2(A, c2);
fill(A, R, R2);
}
Scalar &EGRowEntries::operator()(int i, int j, int k, int l) {
return a[i][j][k * n[j] + l];
}
int CellMatrixGraph::CellDoFs(const Discretization &disc, const cell &c) {
vector<short> z;
disc.NodalDoFs(c, z);
return z[0];
}
void CellMatrixGraph::SeedCells(const Discretization &disc) {
for (cell c = M.cells(); c != M.cells_end(); ++c)
Rows::Insert(c(), CellDoFs(disc, c));