Commit 63d1b81d authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

refactored assemble classes

parent 0024e799
set(assemble
assemble/LagrangeEllipticAssemble.cpp
assemble/MixedEllipticAssemble.cpp
assemble/DGEllipticAssemble.cpp
assemble/HybridEllipticAssemble.cpp
assemble/DGEllipticAssemble.cpp
assemble/DGTransportAssemble.cpp
assemble/DGReactionAssemble.cpp
assemble/PGReactionAssemble.cpp
assemble/DGReactionAssemble.cpp
)
\ No newline at end of file
This diff is collapsed.
#ifndef _DGLAPLACE_H_
#define _DGLAPLACE_H_
#include "LagrangeEllipticAssemble.hpp"
#include "discretization/DGDiscretization.hpp"
#include "discretization/DGUtils.h"
#include "elements/Elements.h"
#include "IEllipticAssemble.hpp"
class DGEllipticAssemble : public LagrangeEllipticAssemble, public DGUtils {
class DGEllipticAssemble : public IEllipticAssemble {
public:
DGDiscretization *disc;
double penalty = 1.0;
double penalty = 20.0;
int sign = 1;
DGEllipticAssemble(IDiscretization *disc, StochasticEllipticProblem *problem)
: LagrangeEllipticAssemble(disc, problem) {
this->disc = dynamic_cast<DGDiscretization *>(disc);
: IEllipticAssemble(problem), disc(dynamic_cast<DGDiscretization *> (disc)) {
this->disc = dynamic_cast<DGDiscretizationT<> *>(disc);
config.get("penalty", penalty);
config.get("sign", sign);
}
const char *Name() const override { return "Discontinuous-Galerkin"; }
Element *createElement(const cell &c, const Vector &u) const override {
return new DGElement(*disc, u, c);
void PrintInfo() const override {
mout.PrintInfo("Assemble", verbose,
PrintInfoEntry("Assemble Name", Name()),
PrintInfoEntry("penalty", penalty),
PrintInfoEntry("sign", sign),
PrintInfoEntry("Problem", problem->Name()));
}
FaceElement *createFaceElement(const cell &c, const Vector &u,
int face) const override {
return new DGFaceElement(*disc, u, c, face);
}
const char *Name() const override;
Scalar getValue(const cell &c, const Vector &u,
Element *elem, int q) const override {
auto dgEllipticElem = dynamic_cast<DGElement *>(elem);
return dgEllipticElem->Value(q, u);
}
void Initialize(Vector &u) const override;
Scalar getFaceValue(const cell &c, const Vector &u,
FaceElement *faceElem, int q, int face) const override {
auto dgEllipticFaceElem = dynamic_cast<DGFaceElement *>(faceElem);
return dgEllipticFaceElem->Value(q, u);
}
double Energy(const Vector &u) const override;
VectorField getFlux(const cell &c, const Vector &u,
Element *elem, int q) const override {
auto dgEllipticElem = dynamic_cast<DGElement *>(elem);
Tensor K = problem->Permeability(c);
return K * dgEllipticElem->Derivative(q, u);
}
void Residual(const cell &c, const Vector &u, Vector &r) const override;
VectorField getFaceFlux(const cell &c, const Vector &u,
FaceElement *faceElem, int q,
int face) const override {
auto dgEllipticFaceElem = dynamic_cast<DGFaceElement *>(faceElem);
Tensor K = problem->Permeability(c);
return K * dgEllipticFaceElem->Derivative(q, u);
}
void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
VectorField getGradient(const cell &c, const Vector &u,
Element *elem, int q) const override {
auto dgEllipticElem = dynamic_cast<DGElement *>(elem);
return dgEllipticElem->Derivative(q, u);
}
virtual double EnergyError(const Vector &u) const override;
void Residual(const cell &c, const Vector &u, Vector &r) const override;
virtual double L2(const Vector &u) const override;
void Jacobi(const cell &c, const Vector &u, Matrix &A) const override;
virtual double H1(const Vector &u) const override;
void SetExactSolution(Vector &u) const override;
virtual double L2Error(const Vector &u) const override;
void AssembleTransfer(TransferMatrix &TM) const override;
virtual double L2CellAverageError(const Vector &u) const override;
virtual double MaxError(const Vector &u) const override;
virtual double FluxError(const Vector &u) const override;
virtual double FaceError(const Vector &u) const override;
FluxPair InflowOutflow(const Vector &u) const override;
FluxPair PrescribedInflowOutflow(const Vector &u) const override;
void plotU(const Vector &u, Meshes &meshes);
FluxPair OutflowLeftRight(const Vector &u) const override;
virtual double GoalFunctional(const Vector &u) const override;
virtual void SetExactSolution(Vector &uEx) const override;
virtual void SetFlux(const Vector &u, Vector &flux) override;
void AssembleTransfer(TransferMatrix &TM) const override;
void plotUex(const Vector &u, Meshes &meshes);
};
#endif
//
// Created by niklas on 10.08.20.
//
#include "DGReactionAssemble.hpp"
void DGReactionAssemble::Dirichlet(double t, Vector &u) const {
// Todo has to be checked
u.ClearDirichletFlags();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
RowBndValues u_c(u, c);
if (!u_c.onBnd()) continue;
int sd = c.Subdomain();
DGElement E(*disc, u, c);
for (int i = 0; i < c.Faces(); ++i) {
if (u_c.bc(i) == -1) continue;
DGFaceElement FE(*disc, u, c, i);
Scalar F = problem->FaceConvection(sd, c, i, FE.QNormal(0),
FE.QPoint(0));
if (F > -1e-6) continue;
for (int j = 0; j < disc->NodalPointsOnFace(c, i); ++j) {
int k = disc->NodalPointOnFace(c, i, j);
u_c(k) = problem->Concentration(t, E.NodalPoint(k));
u_c.D(k) = true;
}
}
}
DirichletConsistent(u);
}
double DGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const {
r = 0;
const Vector &u_old = U_old();
for (cell c = u.cells(); c != u.cells_end(); ++c) {
int sd = c.Subdomain();
RowBndValues r_c(r, c);
DGElement elem(*disc, u, c);
double kappa = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem->Diffusion(sd, z);
kappa = max(kappa, norm(K));
Scalar f = problem->Load(sd, z);
Scalar U = elem.Value(q, u);
VectorField DU = elem.Derivative(q, u);
VectorField K_DU = K * DU;
Scalar Uold = elem.Value(q, u_old);
Scalar R = problem->Reaction(z, U);
VectorField B = problem->CellConvection(sd, c, z);
for (int i = 0; i < elem.dimension(); ++i) {
Scalar U_i = elem.Value(q, i);
VectorField DU_i = elem.Derivative(q, i);
r(c(), i) += w * (K_DU * DU_i
+ (B * DU) * U_i
+ (1.0 / dt_) * (U - Uold) * U_i
- (R + f) * U_i);
}
}
for (int f = 0; f < c.Faces(); ++f) {
Scalar s = penalty * c.FaceArea(f);
DGFaceElement felem(*disc, u, c, f);
if (r.GetMesh().onBndDG(c, f)) {
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar BN = problem->FaceConvection(sd, c, f, N, z);
Tensor K = problem->Diffusion(c.Subdomain(), z);
Scalar U_diff = felem.Value(q, u)
- problem->Concentration(t, z);
VectorField K_DU = K * felem.Derivative(q, u);
if (BN < 0) {
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
+ BN * U_diff * phi_i
+ U_diff *
(NDphi_i - s * phi_i));
}
} else {
Scalar F = problem->Flux(2, z) * N;
for (int i = 0; i < felem.dimension(); ++i) {
Scalar phi_i = felem.Value(q, i);
r(c(), i) -= w * F * phi_i;
}
}
}
} else {
cell cf = r.GetMesh().find_neighbour_cell(c, f);
if (cf() < c()) continue;
int f1 = r.GetMesh().find_neighbour_face_id( c.Face(f), cf);
DGFaceElement felem_1(*disc, r, cf, f1);
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar BN = problem->FaceConvection(sd, c, f, N, z);
Tensor K = problem->Diffusion(c.Subdomain(), z);
Scalar U = felem.Value(q, u);
VectorField K_DU = K * felem.Derivative(q, u);
Point Qf_c = felem.QPoint(q);
int q1 = felem.findQPointID(felem_1, Qf_c);
Tensor K_1 = problem->Diffusion(cf.Subdomain(), z);
Scalar U_1 = felem_1.Value(q1, u);
VectorField K_DU_1 = K_1 * felem_1.Derivative(q1, u);
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 * ((s * U - 0.5 * (K_DU * N)) * phi_i
- 0.5 * U * NDphi_i);
r(c(), i) +=
w * ((-s * U_1 - 0.5 * (K_DU_1 * N)) * phi_i
+ 0.5 * U_1 * NDphi_i);
if (BN < 0) {
r(c(), i) -= w * BN * U * phi_i;
r(c(), i) += w * BN * U_1 * phi_i;
}
}
for (int j = 0; j < felem_1.dimension(); ++j) {
Scalar phi_j = felem_1.Value(q1, j);
Scalar NDphi_j = (K_1 * felem_1.Derivative(q1, j)) * N;
r(cf(), j) +=
w * ((0.5 * K_DU_1 * N + s * U_1) * phi_j
+ 0.5 * U_1 * NDphi_j);
r(cf(), j) += w * ((0.5 * K_DU * N - s * U) * phi_j
- 0.5 * U * NDphi_j);
if (BN > 0) {
r(cf(), j) += w * BN * U_1 * phi_j;
r(cf(), j) -= w * BN * U * phi_j;
}
}
}
}
}
}
Collect(r);
return r.norm();
}
void DGReactionAssemble::Jacobi(double t, const Vector &u, Matrix &A) const {
A = 0;
for (cell c = A.cells(); c != A.cells_end(); ++c) {
int sd = c.Subdomain();
DGElement elem(*disc, A, c);
DGRowEntries A_c(A, c, c);
double kappa = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem->Diffusion(sd, z);
kappa = max(kappa, norm(K));
Scalar U = elem.Value(q, u);
Scalar DR = problem->DerivativeReaction(z, U);
VectorField B = problem->CellConvection(sd, c, z);
for (int i = 0; i < elem.dimension(); ++i) {
Scalar U_i = elem.Value(q, i);
VectorField DU_i = elem.Derivative(q, i);
VectorField K_DU_i = K * DU_i;
for (int j = 0; j < elem.dimension(); ++j) {
Scalar U_j = elem.Value(q, j);
VectorField DU_j = elem.Derivative(q, j);
A_c(i, j) += w * (K_DU_i * DU_j
+ (B * DU_j) * U_i
+ (1.0 / dt_) * U_i * U_j
- DR * U_i * U_j);
}
}
}
BFParts bnd(u.GetMesh(), c);
for (int f = 0; f < c.Faces(); ++f) {
Scalar s = penalty * c.FaceArea(f);
DGFaceElement felem(*disc, u, c, f);
if (u.GetMesh().onBndDG(c, f)) {
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar BN = problem->FaceConvection(sd, c, f, N, z);
if (BN < 0) {
Tensor K = problem->Diffusion(c.Subdomain(), z);
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(i, j) -= w * (NDphi_i * phi_j
+ BN * phi_j * phi_i
+ phi_i *
(NDphi_j - s * phi_j));
}
}
}
}
} else {
cell cf = u.GetMesh().find_neighbour_cell(c, f);
if (cf() < c()) continue;
DGRowEntries A_cf(A, c, cf);
DGRowEntries A_fc(A, cf, c);
DGRowEntries A_ff(A, cf, cf);
int f1 = u.GetMesh().find_neighbour_face_id(c.Face(f), cf);
DGFaceElement felem_1(*disc, u, cf, f1);
for (int q = 0; q < felem.nQ(); ++q) {
double w = felem.QWeight(q);
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar BN = problem->FaceConvection(sd, c, f, N, z);
Tensor K = problem->Diffusion(c.Subdomain(), z);
Point Qf_c = felem.QPoint(q);
int q1 = felem.findQPointID( felem_1, Qf_c);
Tensor K_1 = problem->Diffusion(cf.Subdomain(), z);
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(i, j) += w * (-0.5 * NDphi_i * phi_j
- 0.5 * phi_i * NDphi_j
+ s * phi_i * phi_j);
if (BN < 0)
A_c(i, j) -= w * BN * phi_j * phi_i;
}
for (int j = 0; j < felem_1.dimension(); ++j) {
Scalar phi_j = felem_1.Value(q1, j);
Scalar NDphi_j =
(K_1 * felem_1.Derivative(q1, j)) * N;
A_cf(i, j) += w * (0.5 * NDphi_i * phi_j
- phi_i * 0.5 * NDphi_j
- s * phi_i * phi_j);
A_fc(j, i) += w * (0.5 * NDphi_i * phi_j
- phi_i * 0.5 * NDphi_j
- s * phi_i * phi_j);
if (BN < 0)
A_cf(i, j) += w * BN * phi_j * phi_i;
else
A_fc(j, i) -= w * BN * phi_j * phi_i;
}
}
for (int i = 0; i < felem_1.dimension(); ++i) {
Scalar phi_i = felem_1.Value(q1, i);
Scalar NDphi_i =
(K_1 * felem_1.Derivative(q1, i)) * N;
for (int j = 0; j < felem_1.dimension(); ++j) {
Scalar phi_j = felem_1.Value(q1, j);
Scalar NDphi_j =
(K_1 * felem_1.Derivative(q1, j)) * N;
A_ff(i, j) += w * (0.5 * NDphi_i * phi_j
+ phi_i * 0.5 * NDphi_j
+ s * phi_i * phi_j);
if (BN > 0) A_ff(i, j) += w * BN * phi_j * phi_i;
}
}
}
}
}
}
}
std::pair<double, double> DGReactionAssemble::InFlowOutFlowRate(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;
int sd = c.Subdomain();
for (int f = 0; f < c.Faces(); ++f) {
if (bnd[f] == -1) continue;
DGFaceElement FE(*disc, u, c, f);
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Scalar F = problem->FaceConvection(sd, c, f, FE.QNormal(q),
FE.QPoint(q));
Scalar U = FE.Value(q, u);
if (F > 0) outflow += w * F * U;
else inflow += w * F * U;
}
}
}
return {PPM->Sum(inflow), PPM->Sum(outflow)};
}
void DGReactionAssemble::SetInitialValue(Vector &u) const {
u = 0;
Point z[MaxNodalPoints];
for (cell c = u.cells(); c != u.cells_end(); ++c) {
row r = u.find_row(c());
DGElement Elem(*disc, u, c);
Elem.NodalPoints(z);
double lambda = 1e-9;
for (int j = 0; j < Elem.dimension(); ++j)
u(r, j) = problem->Concentration(0, lambda * c() + (1 - lambda) * z[j]);
}
Accumulate(u);
}
double DGReactionAssemble::Mass(const Vector &u) const {
Scalar e = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
e += w * U;
}
}
return PPM->Sum(e);
}
void DGReactionAssemble::VtkPlotting_cell(double t, const Vector &u, char *filename) const {
Vector u_tmp(u);
for (cell c = u_tmp.cells(); c != u_tmp.cells_end(); ++c) {
DGElement elem(*disc, u_tmp, c);
Scalar U = 0.0;
double a = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
U += w * elem.Value(q, u);
a += w;
}
U *= (1 / a);
row r = u.find_row(c());
u_tmp(r)[0] = U;
}
plot->celldata(u_tmp, 1);
NumberName("U", filename, step);
plot->vtk_celldata(filename, 0);
plot->gnu_celldata(filename, 0);
}
......@@ -6,7 +6,7 @@
#include "IReactionAssemble.hpp"
class DGReactionAssemble : public IReactionAssemble, public DGUtils {
class DGReactionAssemble : public IReactionAssemble {
public:
DGDiscretization *disc;
ReactionProblem *problem;
......
......@@ -75,7 +75,7 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
}
for (int f = 0; f < c.Faces(); ++f) {
DGFaceElement faceElem(*disc, systemMatrix, c, f);
if (checkIfBND(systemMatrix, c, f)) {
if (systemMatrix.GetMesh().onBndDG(c, f)) {
for (int q = 0; q < faceElem.nQ(); ++q) {
const Point &Qf_c = faceElem.QPoint(q);
VectorField Nq = faceElem.QNormal(q);
......@@ -91,8 +91,8 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
}
}
} else {
cell cf = findNBCell(systemMatrix, c, f);
int f1 = findNBFaceID(systemMatrix, c.Face(f), cf);
cell cf = systemMatrix.GetMesh().find_neighbour_cell(c, f);
int f1 = systemMatrix.GetMesh().find_neighbour_face_id(c.Face(f), cf);
DGFaceElement felem_1(*disc, systemMatrix, cf, f1);
DGRowEntries A_cf(systemMatrix, c, cf);
for (int q = 0; q < faceElem.nQ(); ++q) {
......@@ -100,7 +100,7 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
VectorField Nq = faceElem.QNormal(q);
Scalar BN = problem->faceB(c, f, Nq, Qf_c);
if (BN > 0) continue;
int q1 = findQPointID(faceElem, felem_1, Qf_c);
int q1 = felem_1.findQPointID(faceElem, Qf_c);
double w = faceElem.QWeight(q);
for (int i = 0; i < faceElem.dimension(); ++i) {
Scalar phi_i = faceElem.Value(q, i);
......@@ -122,7 +122,7 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
const Point &z = faceElem.QPoint(q);
const Point &N = faceElem.QNormal(q);
const Point &Qf_c = faceElem.QPoint(q);
int q1 = findQPointID(faceElem, felem_1, Qf_c);
int q1 = felem_1.findQPointID(faceElem, Qf_c);
double s = 1;
for (int i = 0; i < faceElem.dimension(); ++i) {
Scalar phi_i = faceElem.Value(q, i);
......@@ -173,7 +173,7 @@ void DGTransportAssemble::RHS(double t, Vector &rhs) const {
for (cell c = rhs.cells(); c != rhs.cells_end(); ++c) {
row r = rhs.find_row(c());
for (int f = 0; f < c.Faces(); ++f) {
if (!checkIfBND(rhs, c, f)) continue;
if (!rhs.GetMesh().onBndDG(c, f)) continue;
DGFaceElement felem(*disc, rhs, c, f);
for (int q = 0; q < felem.nQ(); ++q) {
const Point &z = felem.QPoint(q);
......@@ -213,7 +213,6 @@ void DGTransportAssemble::Energy(const cell &c, const Vector &u, double &energy)
}
}
double DGTransportAssemble::Error(double t, const Vector &u) const {
double err = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
......
......@@ -15,7 +15,7 @@ public:
DGTransportAssemble(IDiscretization *disc, StochasticTransportProblem *problem,
Plot *plot)
: DGTAssemble(dynamic_cast<DGDiscretization *>(disc),
: DGTAssemble(dynamic_cast<DGDiscretizationT<> *>(disc),
plot), problem(problem) {
config.get("flux_alpha", flux_alpha);
}
......@@ -25,7 +25,7 @@ public:
void PrintInfo(const Vector &u) const override {
std::pair<double, double> rate = InFlowOutFlowRate(u);
mout << "n=" << std::setw(3) << std::left << step
<< " t=" << std::setw(9) << std::left << t_
<< " t=" << std::setw(7) << std::left << t_
<< " Energy=" << std::setw(13) << std::left << Energy(u)
<< " Mass=" << std::setw(13) << std::left << Mass(u);
if (abs(rate.first) >= 10e-10)
......@@ -75,12 +75,12 @@ public:
}
void VtkPlotting(double t, const Vector &u) const override {
// const Mesh &Mp = plot->GetMesh();
// const Mesh &Mu = u.GetMesh();