Skip to content

GitLab

Commit 75d5cb0c authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

adapted assemble

parent a1b20a85
Hide whitespace changes
Inline Side-by-side
mlmc/src/assemble/DGEllipticAssemble.cpp
View file @ 75d5cb0c
......@@ -47,7 +47,7 @@ void DGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) con
double w = elem.QWeight(q);
const Point &z = elem.QPoint(q);
Tensor K = problem->Permeability(c);
Scalar Load = problem->Load(c.Subdomain(), z);
Scalar Load = problem->Load(z);
VectorField K_DU = K * elem.Derivative(q, u);
for (int i = 0; i < elem.dimension(); ++i) {
Scalar U_i = elem.Value(q, i);
......@@ -64,7 +64,7 @@ void DGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) con
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
const Point &z = faceElem.QPoint(q);
Scalar F = problem->Flux(bnd, z) * faceElem.QNormal(q);
Scalar F = problem->Flux(z) * faceElem.QNormal(q);
for (int i = 0; i < faceElem.dimension(); ++i) {
Scalar U_i = faceElem.Value(q, i);
r(c(), i) -= w * U_i * F;
......@@ -224,7 +224,7 @@ double DGEllipticAssemble::EnergyError(const Vector &u) const {
Tensor K = problem->Permeability(c);
Tensor IK = Invert(K);
VectorField diff = (K * elem.Derivative(q, u) -
problem->Flux(c.Subdomain(), z));
problem->Flux(z));
err += w * (IK * diff) * diff;
}
}
......@@ -311,7 +311,7 @@ double DGEllipticAssemble::FluxError(const Vector &u) const {
DGFaceElement faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
VectorField G = problem->Flux(bnd[face], faceElem.QPoint(q));
VectorField G = problem->Flux(faceElem.QPoint(q));
Tensor K = problem->Permeability(c);
VectorField DU = faceElem.Derivative(q, u);
Scalar F = (K * DU - G) * faceElem.QNormal(q);
......@@ -381,7 +381,7 @@ FluxPair DGEllipticAssemble::PrescribedInflowOutflow(const Vector &u) const {
DGFaceElement faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar FN = problem->Flux(bnd[face], faceElem.QPoint(q)) *
Scalar FN = problem->Flux(faceElem.QPoint(q)) *
faceElem.QNormal(q);
if (FN > 0) { outflow += w * FN; }
else { inflow += w * FN; }
......
mlmc/src/assemble/DGEllipticAssemble.hpp
View file @ 75d5cb0c
......@@ -12,7 +12,7 @@ public:
double penalty = 20.0;
int sign = 1;
DGEllipticAssemble(IDiscretization *disc, StochasticEllipticProblem *problem)
DGEllipticAssemble(IDiscretization *disc, IStochasticEllipticProblem *problem)
: IEllipticAssemble(problem), disc(dynamic_cast<DGDiscretization *> (disc)) {
this->disc = dynamic_cast<DGDiscretizationT<> *>(disc);
config.get("penalty", penalty);
......
mlmc/src/assemble/DGReactionAssemble.cpp
View file @ 75d5cb0c
......@@ -13,7 +13,7 @@ void DGReactionAssemble::Dirichlet(double t, Vector &u) const {
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),
Scalar F = problem->FaceConvection(c, i, FE.QNormal(0),
FE.QPoint(0));
if (F > -1e-6) continue;
for (int j = 0; j < disc->NodalPointsOnFace(c, i); ++j) {
......@@ -37,15 +37,15 @@ double DGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem->Diffusion(sd, z);
Tensor K = problem->Diffusion(z);
kappa = max(kappa, norm(K));
Scalar f = problem->Load(sd, z);
Scalar f = problem->Load(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);
VectorField B = problem->CellConvection(c, z);
for (int i = 0; i < elem.dimension(); ++i) {
Scalar U_i = elem.Value(q, i);
VectorField DU_i = elem.Derivative(q, i);
......@@ -63,8 +63,8 @@ double DGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const
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 BN = problem->FaceConvection(c, f, N, z);
Tensor K = problem->Diffusion(z);
Scalar U_diff = felem.Value(q, u)
- problem->Concentration(t, z);
VectorField K_DU = K * felem.Derivative(q, u);
......@@ -96,13 +96,13 @@ double DGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const
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 BN = problem->FaceConvection(c, f, N, z);
Tensor K = problem->Diffusion(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);
Tensor K_1 = problem->Diffusion(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) {
......@@ -149,11 +149,11 @@ void DGReactionAssemble::Jacobi(double t, const Vector &u, Matrix &A) const {
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
Tensor K = problem->Diffusion(sd, z);
Tensor K = problem->Diffusion(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);
VectorField B = problem->CellConvection(c, z);
for (int i = 0; i < elem.dimension(); ++i) {
Scalar U_i = elem.Value(q, i);
VectorField DU_i = elem.Derivative(q, i);
......@@ -177,9 +177,9 @@ void DGReactionAssemble::Jacobi(double t, const Vector &u, Matrix &A) const {
double w = felem.QWeight(q);
Point z = felem.QPoint(q);
Point N = felem.QNormal(q);
Scalar BN = problem->FaceConvection(sd, c, f, N, z);
Scalar BN = problem->FaceConvection(c, f, N, z);
if (BN < 0) {
Tensor K = problem->Diffusion(c.Subdomain(), z);
Tensor K = problem->Diffusion(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;
......@@ -206,11 +206,11 @@ void DGReactionAssemble::Jacobi(double t, const Vector &u, Matrix &A) const {
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 BN = problem->FaceConvection(c, f, N, z);
Tensor K = problem->Diffusion(z);
Point Qf_c = felem.QPoint(q);
int q1 = felem.findQPointID( felem_1, Qf_c);
Tensor K_1 = problem->Diffusion(cf.Subdomain(), z);
Tensor K_1 = problem->Diffusion(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;
......@@ -275,8 +275,8 @@ std::pair<double, double> DGReactionAssemble::InFlowOutFlowRate(const Vector &u)
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 F = problem->FaceConvection(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;
......
mlmc/src/assemble/DGReactionAssemble.hpp
View file @ 75d5cb0c
......@@ -9,14 +9,14 @@
class DGReactionAssemble : public IReactionAssemble {
public:
DGDiscretization *disc;
StochasticReactionProblem *problem;
IStochasticReactionProblem *problem;
Plot *plot;
double penalty = 1.0;
int sign = 1;
double flux_alpha = 1.0;
DGReactionAssemble(IDiscretization *disc, StochasticReactionProblem *problem,
DGReactionAssemble(IDiscretization *disc, IStochasticReactionProblem *problem,
Plot *plot)
: disc(dynamic_cast<DGDiscretization *> (disc)),
plot(plot), problem(problem) {
......
mlmc/src/assemble/DGTransportAssemble.cpp
View file @ 75d5cb0c
......@@ -64,7 +64,7 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
DGRowEntries A_c(systemMatrix, c, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
VectorField B = problem->cellB(c, elem.QPoint(q));
VectorField B = problem->CellFlux(c, elem.QPoint(q));
for (int i = 0; i < elem.dimension(); ++i) {
Scalar Phi_i = elem.Value(q, i);
for (int j = 0; j < elem.dimension(); ++j) {
......@@ -79,7 +79,7 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
for (int q = 0; q < faceElem.nQ(); ++q) {
const Point &Qf_c = faceElem.QPoint(q);
VectorField Nq = faceElem.QNormal(q);
Scalar BN = problem->faceB(c, f, Nq, Qf_c);
Scalar BN = problem->FaceNormalFlux(c, f, Nq, Qf_c);
if (BN > 0) continue;
double w = faceElem.QWeight(q);
for (int i = 0; i < faceElem.dimension(); ++i) {
......@@ -98,7 +98,7 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
for (int q = 0; q < faceElem.nQ(); ++q) {
const Point &Qf_c = faceElem.QPoint(q);
VectorField Nq = faceElem.QNormal(q);
Scalar BN = problem->faceB(c, f, Nq, Qf_c);
Scalar BN = problem->FaceNormalFlux(c, f, Nq, Qf_c);
if (BN > 0) continue;
int q1 = felem_1.findQPointID(faceElem, Qf_c);
double w = faceElem.QWeight(q);
......@@ -169,7 +169,7 @@ void DGTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
void DGTransportAssemble::RHS(double t, Vector &rhs) const {
rhs = 0;
if (!problem->RHS()) return;
// if (!problem->RHS()) return;
for (cell c = rhs.cells(); c != rhs.cells_end(); ++c) {
row r = rhs.find_row(c());
for (int f = 0; f < c.Faces(); ++f) {
......@@ -179,9 +179,9 @@ void DGTransportAssemble::RHS(double t, Vector &rhs) const {
const Point &z = felem.QPoint(q);
double w = felem.QWeight(q);
VectorField N = felem.QNormal(q);
Scalar BN = problem->faceB(c, f, N, z);
Scalar BN = problem->FaceNormalFlux(c, f, N, z);
if (BN > 0) continue;
Scalar U = problem->ut(t, z);
Scalar U = problem->Solution(t, z);
for (int i = 0; i < felem.dimension(); ++i) {
Scalar Phi_i = felem.Value(q, i);
rhs(r, i) -= w * U * BN * Phi_i;
......@@ -220,7 +220,7 @@ double DGTransportAssemble::Error(double t, const Vector &u) const {
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
Scalar Sol = problem->ut(t, elem.QPoint(q));
Scalar Sol = problem->Solution(t, elem.QPoint(q));
err += w * (U - Sol) * (U - Sol);
}
}
......@@ -248,7 +248,7 @@ std::pair<double, double> DGTransportAssemble::InFlowOutFlowRate(const Vector &u
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
VectorField Nq = faceElem.QNormal(q);
Scalar BN = problem->faceB(c, f, Nq, elem.QPoint(q));
Scalar BN = problem->FaceNormalFlux(c, f, Nq, elem.QPoint(q));
if (BN > 0) outflow += w * BN * U;
else inflow += w * BN * U;
}
......@@ -271,7 +271,7 @@ void DGTransportAssemble::SetExactSolution(double t, Vector &u_ex) const {
DGElement Elem(*disc, u_ex, c);
Elem.NodalPoints(z);
for (int j = 0; j < Elem.dimension(); ++j)
u_ex(r, j) = problem->ut(t, z[j]);
u_ex(r, j) = problem->Solution(t, z[j]);
}
Accumulate(u_ex);
}
......
mlmc/src/assemble/DGTransportAssemble.hpp
View file @ 75d5cb0c
......@@ -4,6 +4,8 @@
#include "timestepping/DGTAssemble.hpp"
#include "problem/StochasticTransportProblem.hpp"
#include "TimeSeries.h"
#include "utility/Logging.hpp"
#include "elements/DGFaceElement.h"
class DGTransportAssemble : public DGTAssemble {
......@@ -23,22 +25,22 @@ public:
const char *Name() const override { return "DGTransportAssemble"; }
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(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)
mout << " IFR=" << std::setw(13) << std::left << rate.first;
if (abs(rate.second) >= 10e-10)
mout << " OFR=" << std::setw(13) << std::left << rate.second;
mout << endl;
// std::pair<double, double> rate = InFlowOutFlowRate(u);
// mout << "n=" << std::setw(3) << std::left << step
// << " 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)
// mout << " IFR=" << std::setw(13) << std::left << rate.first;
// if (abs(rate.second) >= 10e-10)
// mout << " OFR=" << std::setw(13) << std::left << rate.second;
// mout << endl;
}
void PrintInfo() const {
mout.PrintInfo("Assemble", verbose,
PrintInfoEntry("Name", Name()),
PrintInfoEntry("Problem", problem->Name()),
// PrintInfoEntry("Problem", problem->Name()),
PrintInfoEntry("Discretization", disc->DiscName()));
}
......@@ -69,18 +71,18 @@ public:
Elem.NodalPoints(z);
double lambda = 1e-9;
for (int j = 0; j < Elem.dimension(); ++j)
u(r, j) = problem->ut(0, lambda * c() + (1 - lambda) * z[j]);
u(r, j) = problem->Solution(0, lambda * c() + (1 - lambda) * z[j]);
}
Accumulate(u);
}
void VtkPlotting(double t, const Vector &u) const override {
const Mesh &Mp = plot->GetMesh();
const Mesh &Mu = u.GetMesh();
if (Mp.CellCount() != Mu.CellCount()) return;
vout (2) << " => VtkPlotting result of step " << step << ".\n";
char filename[128];
VtkPlotting_cell(t, u, filename);
// const Mesh &Mp = plot->GetMesh();
// const Mesh &Mu = u.GetMesh();
// if (Mp.CellCount() != Mu.CellCount()) return;
// vout (2) << " => VtkPlotting result of step " << step << ".\n";
// char filename[128];
// VtkPlotting_cell(t, u, filename);
}
void PrintMatrixInfo(Matrix &A, int diagonal) const override {
......
mlmc/src/assemble/HybridEllipticAssemble.cpp
View file @ 75d5cb0c
......@@ -58,7 +58,7 @@ void HybridEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r)
RT0_LagrangeFaceElementT faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar F = problem->Flux(bnd, faceElem.QPoint(q)) * faceElem.QNormal(q);
Scalar F = problem->Flux(faceElem.QPoint(q)) * faceElem.QNormal(q);
r_c(face) -= w * F;
}
}
......@@ -156,9 +156,9 @@ FluxPair HybridEllipticAssemble::OutflowLeftRight(const Vector &p) const {
void HybridEllipticAssemble::SetFlux(const Vector &u, Vector &flux) {
for (cell c = u.cells(); c != u.cells_end(); ++c) {
int N = c.Faces();
SmallMatrix A(N);
SmallVector B(N), R(N);
int Faces = c.Faces();
SmallMatrix A(Faces);
SmallVector B(Faces), R(Faces);
LocalProblem(c, u, A, B, R);
InverseSmallMatrix IA(A);
SmallVector IAB(IA, B);
......@@ -166,10 +166,10 @@ void HybridEllipticAssemble::SetFlux(const Vector &u, Vector &flux) {
SmallVector RU(R);
SmallVector UU(R);
RT0_LagrangeElementT<> elem(*disc, u, c);
for (int i = 0; i < c.Faces(); ++i)
UU[i] = u(elem[i], 0);
for (int i = 0; i < c.Faces(); ++i)
RU[i] *= UU[i];
for (int face = 0; face < c.Faces(); ++face)
UU[face] = u(elem[face], 0);
for (int face = 0; face < c.Faces(); ++face)
RU[face] *= UU[face];
SmallVector IARU(IA, RU);
Scalar BIARU = B * IARU;
Scalar P = BIARU / BIAB;
......@@ -364,7 +364,7 @@ double HybridEllipticAssemble::EnergyError(const Vector &u) const {
Point z = elem.QPoint(q);
Tensor IK = Invert(problem->Permeability(c));
VectorField diff = (Flux(c, u, elem, q) -
problem->Flux(c.Subdomain(), z));
problem->Flux(z));
err += w * (IK * diff) * diff;
}
}
......@@ -452,7 +452,7 @@ double HybridEllipticAssemble::FluxError(const Vector &u) const {
RT0_LagrangeFaceElementT<> faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
VectorField G = problem->Flux(bnd[face], faceElem.QPoint(q));
VectorField G = problem->Flux(faceElem.QPoint(q));
Scalar F = (FaceFlux(c, u, faceElem, q, face) - G) *
faceElem.QNormal(q);
flux_error += w * F * F;
......@@ -498,7 +498,7 @@ FluxPair HybridEllipticAssemble::PrescribedInflowOutflow(const Vector &u) const
RT0_LagrangeFaceElementT<> faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar F = problem->Flux(bnd[face], faceElem.QPoint(q)) *
Scalar F = problem->Flux(faceElem.QPoint(q)) *
faceElem.QNormal(q);
if (F > 0) { outflow += w * F; }
else { inflow += w * F; }
......
mlmc/src/assemble/LagrangeEllipticAssemble.cpp
View file @ 75d5cb0c
......@@ -46,7 +46,7 @@ void LagrangeEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &
double w = elem.QWeight(q);
const Point &z = elem.QPoint(q);
Tensor K = problem->Permeability(c);
Scalar f = problem->Load(c.Subdomain(), z);
Scalar f = problem->Load(z);
VectorField DU = elem.Derivative(q, u);
VectorField K_DU = K * DU;
for (int i = 0; i < elem.size(); ++i) {
......@@ -63,7 +63,7 @@ void LagrangeEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &
RowValues r_f(r, faceElem);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar F = problem->Flux(bnd, faceElem.QPoint(q)) * faceElem.QNormal(q);
Scalar F = problem->Flux(faceElem.QPoint(q)) * faceElem.QNormal(q);
for (int i = 0; i < faceElem.size(); ++i) {
Scalar U_i = faceElem.Value(q, i);
r_f(i) -= w * U_i * F;
......@@ -98,7 +98,7 @@ double LagrangeEllipticAssemble::EnergyError(const Vector &u) const {
Tensor K = problem->Permeability(c);
Tensor IK = Invert(K);
VectorField diff = (K * elem.Derivative(q, u) -
problem->Flux(c.Subdomain(), z));
problem->Flux(z));
err += w * (IK * diff) * diff;
}
}
......@@ -185,7 +185,7 @@ double LagrangeEllipticAssemble::FluxError(const Vector &u) const {
ScalarFaceElement faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
VectorField G = problem->Flux(bnd[face], faceElem.QPoint(q));
VectorField G = problem->Flux(faceElem.QPoint(q));
Tensor K = problem->Permeability(c);
VectorField DU = faceElem.Derivative(q, u);
Scalar F = (K * DU - G) * faceElem.QNormal(q);
......@@ -255,7 +255,7 @@ FluxPair LagrangeEllipticAssemble::PrescribedInflowOutflow(const Vector &u) cons
ScalarFaceElement faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar F = problem->Flux(bnd[face], faceElem.QPoint(q)) *
Scalar F = problem->Flux(faceElem.QPoint(q)) *
faceElem.QNormal(q);
if (F > 0) { outflow += w * F; }
else { inflow += w * F; }
......
mlmc/src/assemble/MixedEllipticAssemble.cpp
View file @ 75d5cb0c
......@@ -16,7 +16,7 @@ void MixedEllipticAssemble::Initialize(Vector &u) const {
if (bnd[face] == 2) {
for (int j = 0; j < faceElem.VelocitySizeOnFace(); j++) {
double area = faceElem.Area();
VectorField F = problem->Flux(c.Subdomain(), faceElem[0]());
VectorField F = problem->Flux(faceElem[0]());
VectorField N = faceElem.Normal();
uf_c(faceElem.VelocityIndex(), j) = area * (F * N);
uf_c.D(faceElem.VelocityIndex(), j) = true;
......@@ -51,7 +51,7 @@ void MixedEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r)
MixedRowBndValues r_c(r, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar f = problem->Load(c.Subdomain(), elem.QPoint(q));
Scalar f = problem->Load(elem.QPoint(q));
VectorField Q = elem.VelocityField(q, u);
Tensor IK = Invert(problem->Permeability(c));
Scalar divQ = elem.VelocityFieldDivergence(q, u);
......@@ -125,7 +125,7 @@ double MixedEllipticAssemble::EnergyError(const Vector &u) const {
Point z = elem.QPoint(q);
Tensor IK = Invert(problem->Permeability(c));
VectorField Q = elem.VelocityField(q, u);
VectorField F = problem->Flux(c.Subdomain(), z);
VectorField F = problem->Flux(z);
VectorField Diff = Q - F;
err += w * IK * Diff * Diff;
}
......@@ -215,7 +215,7 @@ double MixedEllipticAssemble::FluxError(const Vector &u) const {
RT0_LagrangeFaceElementT<> faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
VectorField F = problem->Flux(bnd[face], faceElem.QPoint(q));
VectorField F = problem->Flux(faceElem.QPoint(q));
Scalar Diff = (faceElem.VectorValue(q, u) - F) * faceElem.QNormal(q);
flux_error += w * Diff * Diff;
}
......@@ -268,7 +268,7 @@ FluxPair MixedEllipticAssemble::PrescribedInflowOutflow(const Vector &u) const {
RT0_LagrangeFaceElementT<> faceElem(*disc, u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar FN = problem->Flux(bnd[face], faceElem.QPoint(q)) *
Scalar FN = problem->Flux(faceElem.QPoint(q)) *
faceElem.QNormal(q);
if (FN > 0) { outflow += w * FN; }
else { inflow += w * FN; }
......@@ -337,7 +337,7 @@ void MixedEllipticAssemble::SetExactSolution(Vector &u) const {
vector<Point> nodalPointsOnFace;
u.NodalPointsOnFace(faceElem.VelocityIndex(), c, face, nodalPointsOnFace);
for (int j = 0; j < faceElem.VelocitySizeOnFace(); j++) {
VectorField F = problem->Flux(c.Subdomain(), nodalPointsOnFace[j]);
VectorField F = problem->Flux(nodalPointsOnFace[j]);
uf_c(faceElem.VelocityIndex(), j) = faceElem.Area() * F * N;
}
}
......
mlmc/src/assemble/PGReactionAssemble.cpp
View file @ 75d5cb0c
......@@ -11,7 +11,7 @@ void PGReactionAssemble::Dirichlet(double t, Vector &u) const {
for (int i = 0; i < c.Faces(); ++i) {
if (u_c.bc(i) == -1) continue;
ScalarFaceElement FE(*disc, u, c, i);
Scalar F = problem->FaceConvection(sd, c, i, FE.QNormal(0),
Scalar F = problem->FaceConvection(c, i, FE.QNormal(0),
FE.QPoint(0));
if (F > -1e-6) continue;
for (int j = 0; j < disc->NodalPointsOnFace(c, i); ++j) {
......@@ -35,10 +35,10 @@ double PGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const
double w = E.QWeight(q);
Scalar U = E.Value(q, u);
Scalar Uold = E.Value(q, u_old);
Tensor K = problem->Diffusion(c.Subdomain(), E.QPoint(q));
Scalar f = problem->Load(c.Subdomain(), E.QPoint(q));
VectorField C = problem->CellConvection(c.Subdomain(),
c, E.QPoint(q));
Tensor K = problem->Diffusion(E.QPoint(q));
Scalar f = problem->Load(E.QPoint(q));
VectorField C = problem->CellConvection(
c, E.QPoint(q));
Scalar R = problem->Reaction(E.QPoint(q), U);
VectorField DU = E.Derivative(q, u);
for (unsigned int i = 0; i < E.size(); ++i) {
......@@ -60,7 +60,7 @@ double PGReactionAssemble::Residual(double t, const Vector &u, Vector &r) const
RowValues r_f(r, FE);
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Scalar F = problem->FaceConvection(sd, c, f, FE.QNormal(q),
Scalar F = problem->FaceConvection(c, f, FE.QNormal(q),
FE.QPoint(q));
if (F < 1e-6) continue;
Scalar U = 0;
......@@ -90,9 +90,9 @@ void PGReactionAssemble::Jacobi(double t, const Vector &u, Matrix &A) const {
for (int q = 0; q < E.nQ(); ++q) {
double w = E.QWeight(q);
Scalar U = E.Value(q, u);
Tensor K = problem->Diffusion(c.Subdomain(), E.QPoint(q));
VectorField C = problem->CellConvection(c.Subdomain(),
c, E.QPoint(q));
Tensor K = problem->Diffusion(E.QPoint(q));
VectorField C = problem->CellConvection(
c, E.QPoint(q));
Scalar DR = problem->DerivativeReaction(E.QPoint(q), U);
for (unsigned int i = 0; i < E.size(); ++i) {
VectorField DU_i = E.Derivative(q, i);
......@@ -117,7 +117,7 @@ void PGReactionAssemble::Jacobi(double t, const Vector &u, Matrix &A) const {
RowEntries A_f(A, FE);
for (int q = 0; q < FE.nQ(); ++q) {
double w = FE.QWeight(q);
Scalar F = problem->FaceConvection(sd, c, f, FE.QNormal(q),
Scalar F = problem->FaceConvection(c, f, FE.QNormal(q),
FE.QPoint(q));
if (F < 1e-6) continue;
for (unsigned int i = 0; i < FE.size(); ++i) {
......@@ -146,7 +146,7 @@ std::pair<double, double> PGReactionAssemble::InFlowOutFlowRate(const Vector &u)
ScalarFaceElement 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),
Scalar F = problem->FaceConvection(c, f, FE.QNormal(q),
FE.QPoint(q));
Scalar U = FE.Value(q, u);
if (F > 0) outflow += w * F * U;
......@@ -189,8 +189,8 @@ void PGReactionAssemble::VtkPlotting_cell(double t, const Vector &u, char *filen
double PGReactionAssemble::Delta(const cell &c) const {
double h = diam(c);
Tensor K = problem->Diffusion(c.Subdomain(), c());
VectorField C = problem->CellConvection(c.Subdomain(), c, c());
Tensor K = problem->Diffusion(c());
VectorField C = problem->CellConvection(c, c());
double Pe = 0.5 * h * norm(C) / norm(K);
if (Pe > 1) return h * delta;
return h * h * delta / norm(K);
......
mlmc/src/assemble/PGReactionAssemble.hpp
View file @ 75d5cb0c
......@@ -12,12 +12,12 @@
class PGReactionAssemble : public IReactionAssemble {
public:
IDiscretization *disc;
StochasticReactionProblem *problem;
IStochasticReactionProblem *problem;