Skip to content

GitLab

Commit 145af950 authored by jonathan.froehlich's avatar jonathan.froehlich
Browse files

[426-make-egdisc-a-mixeddisc] Started adaption of assemble

parent 1a335d97
Hide whitespace changes
Inline Side-by-side
tutorial/src/assemble/EGEllipticAssemble.cpp
View file @ 145af950
#include "EGEllipticAssemble.hpp"
#include "MixedEGElement.hpp"
const char *EGEllipticAssemble::Name() const {
return "EGEllipticAssemble";
......@@ -9,21 +10,21 @@ void EGEllipticAssemble::Initialize(Vector &u) const {
ExchangeBuffer EB;
Point z = FindDirichlet(u);
if (z != infty)
for (int q = 0; q < PPM->size(); ++q)
for (int q = 0; q < PPM->size(); ++q)
EB.Send(q) << z;
EB.Communicate();
z = infty;
for (int q = 0; q < PPM->size(); ++q) {
while (EB.Receive (q).size() < EB.Receive (q).Size()) {
EB.Receive (q) >> z;
while (EB.Receive(q).size() < EB.Receive(q).Size()) {
EB.Receive(q) >> z;
break;
}
if (z != infty) break;
}
row r = u.find_row(z);
if (r != u.rows_end()) {
u(r,0) = problem.Solution(z);
u.D(r,0) = true;
u(r, 0) = problem.Solution(z);
u.D(r, 0) = true;
}
u.DirichletConsistent();
}
......@@ -31,7 +32,7 @@ void EGEllipticAssemble::Initialize(Vector &u) const {
double EGEllipticAssemble::Energy(const Vector &u) const {
double energy = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
EGElement elem(u, c);
MixedEGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
VectorField DU = elem.Derivative(q, u);
......@@ -44,22 +45,23 @@ double EGEllipticAssemble::Energy(const Vector &u) const {
}
void EGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) const {
EGElement elem(u, c);
RowBndValues r_c(r, c);
MixedEGElement elem(u, c);
MixedRowBndValues r_c(r, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem.Permeability(z);
Scalar f = problem.Load(z);
VectorField K_DU = K * elem.Derivative(q, u);
for (int i = 0; i < elem.NodalPoints(); ++i) {
for (int i = 0; i < elem.ValueSize(); ++i) {
Scalar U_i = elem.Value(q, i);
VectorField DU_i = elem.Derivative(q, i);
r_c(i) += w * (K_DU * DU_i - f * U_i);
r_c(elem.ValueIndex(), i) += w * (K_DU * DU_i - f * U_i);
}
r_c(elem.PenaltyIndex(), 0) += w * elem.PenaltyValue(q);
}
for (int f = 0; f < c.Faces(); ++f) {
EGFaceElement felem(u, c, f);
MixedEGFaceElement felem(u, c, f);
if (r.GetMesh().onBndDG(c, f)) {
int bnd = r_c.bc(f);
if (bnd == 2) {
......@@ -68,9 +70,9 @@ void EGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) con
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar F = problem.Flux(z) * N;
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
r_c(i) -= w * F * phi_i;
for (int i = 0; i < felem.ValueSize(); ++i) {
Scalar phi_i = felem.Value(q, i);
r_c(elem.ValueIndex, i) -= w * F * phi_i;
}
}
} else if (bnd == 1) {
......@@ -82,95 +84,93 @@ void EGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) con
Scalar U_diff = felem.Value(q, u) - problem.Solution(z);
VectorField K_DU = K * felem.Derivative(q, u);
Scalar F = K_DU * N;
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q,i)) * N;
r_c(i) -= w * (F * phi_i
- U_diff * (sign * NDphi_i
+ (penalty / c.FaceArea(f)) * (K * N) * N * phi_i));
}
}
for (int i = 0; i < felem.ValueSize(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
r_c(felem.ValueIndex(), i) -= w * (F * phi_i
- U_diff * (sign * NDphi_i
+
(penalty / c.FaceArea(f)) * (K * N) * N *
phi_i));
}
}
}
} else {
cell cf = r.GetMesh().find_neighbour_cell(c, f);
if (cf() < c()) continue;
EGElement elem_1(r, cf);
int nV = elem.NodalPoints() - 1;
int nV_1 = elem_1.NodalPoints() - 1;
if (cf() < c()) continue;
MixedEGElement elem_1(r, cf);
int f1 = r.GetMesh().find_neighbour_face_id(c.Face(f), cf);
EGFaceElement felem_1(r, cf, f1);
RowBndValues r_cf(r, cf);
MixedEGFaceElement felem_1(r, cf, f1);
MixedRowBndValues r_cf(r, cf);
Tensor K = problem.Permeability(c());
Tensor K_1 = problem.Permeability(cf());
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight (q);
Point N = felem.QNormal (q);
double K_e = 2 * (K * N) * N * (K_1 * N) * N / ((K * N) * N + (K_1 * N) * N);
Scalar K_DU = K_e * ((felem.Derivative (q, u) + felem_1.Derivative (q, u)) * N);
r_c(nV) -= w * 0.5 * K_DU;
r_cf(nV) += w * 0.5 * K_DU;
Scalar U = felem.Value (q, u) - felem_1.Value (q, u);
Scalar s = K_e * penalty / c.FaceArea(f);
r_c(nV) += w * s * U;
r_cf(nV) -= w * s * U;
for (int i = 0; i < nV; ++i) {
Scalar NDphi_i = K_e * (felem.Derivative (q, i) * N);
r_c(i) -= w * 0.5 * sign * U * NDphi_i;
}
for (int j = 0; j < nV; ++j) {
Scalar NDphi_1_j= K_e * (felem_1.Derivative (q, j) * N);
r_cf(j) -= w * 0.5 * sign * U * NDphi_1_j;
}
double w = felem.QWeight(q);
Point N = felem.QNormal(q);
double K_e = 2 * (K * N) * N * (K_1 * N) * N / ((K * N) * N + (K_1 * N) * N);
Scalar K_DU = K_e * ((felem.Derivative(q, u) + felem_1.Derivative(q, u)) * N);
r_c(elem.PenaltyIndex(), 0) -= w * 0.5 * K_DU;
r_cf(elem.PenaltyIndex(), 0) += w * 0.5 * K_DU;
Scalar U = felem.Value(q, u) - felem_1.Value(q, u);
Scalar s = K_e * penalty / c.FaceArea(f);
r_c(elem.PenaltyIndex(), 0) += w * s * U;
r_cf(elem.PenaltyIndex(), 0) -= w * s * U;
for (int i = 0; i < elem.ValueSize(); ++i) {
Scalar NDphi_i = K_e * (felem.Derivative(q, i) * N);
r_c(elem.ValueIndex(), i) -= w * 0.5 * sign * U * NDphi_i;
}
for (int j = 0; j < elem.ValueSize(); ++j) {
Scalar NDphi_1_j = K_e * (felem_1.Derivative(q, j) * N);
r_cf(elem.ValueIndex(), j) -= w * 0.5 * sign * U * NDphi_1_j;
}
}
}
}
}
void EGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const {
EGElement elem(A, c);
EGRowEntries A_c(A, c);
MixedEGElement elem(A, c);
MixedRowEntries A_c(A, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem.Permeability(z);
for (int i = 0; i < elem.NodalPoints(); ++i) {
for (int i = 0; i < elem.ValueSize(); ++i) {
VectorField K_DU_i = K * elem.Derivative(q, i);
for (int j = 0; j < elem.NodalPoints(); ++j) {
for (int j = 0; j < elem.ValueSize(); ++j) {
VectorField DU_j = elem.Derivative(q, j);
A_c(i, j) += w * K_DU_i * DU_j;
A_c(elem.ValueIndex(), elem.ValueIndex(), i, j) += w * K_DU_i * DU_j;
}
}
}
BFParts bnd(u.GetMesh(), c);
for (int f = 0; f < c.Faces(); ++f) {
EGFaceElement felem(u, c, f);
MixedEGFaceElement felem(u, c, f);
if (elem.Bnd(f) == 1) {
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
const Point &z = felem.QPoint(q);
const Point &N = felem.QNormal(q);
Tensor K = problem.Permeability(c());
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
for (int j = 0; j < felem.dimension(); ++j) {
Scalar phi_j = felem.Value(q, j);
Scalar NDphi_j = (K * felem.Derivative(q, j)) * N;
A_c(j, i) -= w * (NDphi_i * phi_j
- phi_i * (sign * NDphi_j
+ (penalty / c.FaceArea(f)) * (K * N) * N * phi_j));
}
}
double w = felem.QWeight(q);
const Point &z = felem.QPoint(q);
const Point &N = felem.QNormal(q);
Tensor K = problem.Permeability(c());
for (int i = 0; i < felem.ValueSize(); ++i) {
Scalar phi_i = felem.Value(q, i);
Scalar NDphi_i = (K * felem.Derivative(q, i)) * N;
for (int j = 0; j < felem.ValueSize(); ++j) {
Scalar phi_j = felem.Value(q, j);
Scalar NDphi_j = (K * felem.Derivative(q, j)) * N;
A_c(elem.ValueIndex(), elem.ValueIndex(),j, i) -= w * (NDphi_i * phi_j
- phi_i * (sign * NDphi_j
+ (penalty / c.FaceArea(f)) * (K * N) * N * phi_j));
}
}
}
} else if (elem.Bnd(f) == -1) {
cell cf = u.GetMesh().find_neighbour_cell(c, f);
if (cf() < c()) continue;
EGElement elem_1(A, cf);
EGRowEntriesNB A_cf(A, c, cf);
EGRowEntriesNB A_fc(A, cf, c);
EGRowEntries A_f(A, cf);
int nV = elem.NodalPoints() - 1;
int nV_1 = elem_1.NodalPoints() - 1;
MixedEGElement elem_1(A, cf);
MixedRowEntries A_cf(A, c, cf);
MixedRowEntries A_fc(A, cf, c);
MixedRowEntries A_f(A, cf);
int f1 = u.GetMesh().find_neighbour_face_id(c.Face(f), cf);
EGFaceElement felem_1(u, cf, f1);
Tensor K = problem.Permeability(c());
......@@ -178,26 +178,26 @@ void EGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
Point N = felem.QNormal(q);
double K_e = 2 * (K * N) * N * (K_1 * N) * N / ((K * N) * N + (K_1 * N) * N);
Scalar s = K_e * penalty / c.FaceArea(f);
for (int i = 0; i < nV; ++i) {
Scalar NDphi_i = K_e * (felem.Derivative (q, i) * N);
A_c(nV,i) -= 0.5 * w * NDphi_i;
A_fc(i) += 0.5 * w * NDphi_i;
A_c(i,nV) -= 0.5 * w * sign * NDphi_i;
A_cf(i) += 0.5 * w * sign * NDphi_i;
}
for (int i = 0; i < nV_1; ++i) {
Scalar NDphi_1_i = K_e * (felem_1.Derivative (q, i) * N);
A_cf(i) -= 0.5 * w * NDphi_1_i;
A_f(nV,i) += 0.5 * w * NDphi_1_i;
A_c(i,nV) += 0.5 * w * sign * NDphi_1_i;
A_cf(i) -= 0.5 * w * sign * NDphi_1_i;
}
A_c(nV,nV) += w * s;
A_cf(nV) -= w * s;
A_fc(nV) -= w * s;
A_f(nV,nV) += w * s;
double K_e = 2 * (K * N) * N * (K_1 * N) * N / ((K * N) * N + (K_1 * N) * N);
Scalar s = K_e * penalty / c.FaceArea(f);
for (int i = 0; i < nV; ++i) {
Scalar NDphi_i = K_e * (felem.Derivative(q, i) * N);
A_c(nV, i) -= 0.5 * w * NDphi_i;
A_fc(i) += 0.5 * w * NDphi_i;
A_c(i, nV) -= 0.5 * w * sign * NDphi_i;
A_cf(i) += 0.5 * w * sign * NDphi_i;
}
for (int i = 0; i < nV_1; ++i) {
Scalar NDphi_1_i = K_e * (felem_1.Derivative(q, i) * N);
A_cf(i) -= 0.5 * w * NDphi_1_i;
A_f(nV, i) += 0.5 * w * NDphi_1_i;
A_c(i, nV) += 0.5 * w * sign * NDphi_1_i;
A_cf(i) -= 0.5 * w * sign * NDphi_1_i;
}
A_c(nV, nV) += w * s;
A_cf(nV) -= w * s;
A_fc(nV) -= w * s;
A_f(nV, nV) += w * s;
}
}
}
......
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

Imprint