Commit b9cca575 authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

removed *disc in elements

parent e9dd726c
......@@ -2,53 +2,53 @@
#include "pdesolver/PDESolverCreator.hpp"
#include "pdesolver/assembling/elliptic/HybridEllipticAssemble.hpp"
void HybridFaceNormalFluxGenerator::createPDESolver() {
string problemName;
if (meshes.dim() == 1)
problemName = "StochasticLaplace1D";
if (meshes.dim() == 2)
problemName = "StochasticLaplace2D";
pdeSolver = PDESolverCreator().
WithModel("HybridElliptic").
WithProblem(problemName).
WithQuantity("L2"). // Todo flux error and check?
Create(meshes);
string problemName;
if (meshes.dim() == 1)
problemName = "StochasticLaplace1D";
if (meshes.dim() == 2)
problemName = "StochasticLaplace2D";
pdeSolver = PDESolverCreator().
WithModel("HybridElliptic").
WithProblem(problemName).
WithQuantity("L2"). // Todo flux error and check?
Create(meshes);
}
// Todo delete sample solutions
void HybridFaceNormalFluxGenerator::drawSample(const SampleID &id) {
pdeSolver->DrawSample(id);
solutionFaceValues = new SampleSolution(pdeSolver->GetDisc(), id);
solutionFaceFlux = new SampleSolution(pdeSolver->GetDisc(), id);
pdeSolver->Run(*solutionFaceValues);
auto assemble = dynamic_cast<HybridEllipticAssemble *>(pdeSolver->GetAssemble());
assemble->SetNormalFlux(solutionFaceValues->U, solutionFaceFlux->U);
pdeSolver->DrawSample(id);
solutionFaceValues = new SampleSolution(pdeSolver->GetDisc(), id);
solutionFaceFlux = new SampleSolution(pdeSolver->GetDisc(), id);
pdeSolver->Run(*solutionFaceValues);
auto assemble = dynamic_cast<HybridEllipticAssemble *>(pdeSolver->GetAssemble());
assemble->SetNormalFlux(solutionFaceValues->U, solutionFaceFlux->U);
}
Scalar HybridFaceNormalFluxGenerator::EvalSample(int face, const cell &c) {
RTLagrangeElement elem(*pdeSolver->GetDisc(), solutionFaceFlux->U, c);
RTFaceElement faceElem(*pdeSolver->GetDisc(), solutionFaceFlux->U, c, face);
return (solutionFaceFlux->U)(elem[face], 0) * elem.Sign(face) / faceElem.Area();
RTLagrangeElement elem(solutionFaceFlux->U, c);
RTFaceElement faceElem(solutionFaceFlux->U, c, face);
return (solutionFaceFlux->U)(elem[face], 0) * elem.Sign(face) / faceElem.Area();
}
HybridFaceNormalFluxGenerator::~HybridFaceNormalFluxGenerator() {
if (!pdeSolver) delete pdeSolver;
if (!solutionFaceValues) delete solutionFaceValues;
if (!solutionFaceFlux) delete solutionFaceFlux;
if (!pdeSolver) delete pdeSolver;
if (!solutionFaceValues) delete solutionFaceValues;
if (!solutionFaceFlux) delete solutionFaceFlux;
}
VectorField HybridCellFluxGenerator::EvalSample(const cell &c) {
RTLagrangeElement elem(*generator.pdeSolver->GetDisc(),
generator.solutionFaceFlux->U, c);
VectorField F = zero;
double area = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
area += w;
F += w * elem.VelocityField(q, generator.solutionFaceFlux->U);
}
F *= (1 / area);
return F;
RTLagrangeElement elem(generator.solutionFaceFlux->U, c);
VectorField F = zero;
double area = 0;
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
area += w;
F += w * elem.VelocityField(q, generator.solutionFaceFlux->U);
}
F *= (1 / area);
return F;
}
......@@ -5,7 +5,7 @@
void DGLinearTransportAssemble::MassMatrix(Matrix &massMatrix) const {
massMatrix = 0;
for (cell c = massMatrix.cells(); c != massMatrix.cells_end(); ++c) {
DGElement elem(*disc, massMatrix, c);
DGElement elem(massMatrix, c);
DGRowEntries M_c(massMatrix, c, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
......@@ -23,7 +23,7 @@ void DGLinearTransportAssemble::MassMatrix(Matrix &massMatrix) const {
void DGLinearTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
systemMatrix = 0;
for (cell c = systemMatrix.cells(); c != systemMatrix.cells_end(); ++c) {
DGElement elem(*disc, systemMatrix, c);
DGElement elem(systemMatrix, c);
DGRowEntries A_c(systemMatrix, c, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
......@@ -37,7 +37,7 @@ void DGLinearTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
}
}
for (int f = 0; f < c.Faces(); ++f) {
DGFaceElement faceElem(*disc, systemMatrix, c, f);
DGFaceElement faceElem(systemMatrix, c, f);
if (systemMatrix.GetMesh().onBndDG(c, f)) {
for (int q = 0; q < faceElem.nQ(); ++q) {
const Point &Qf_c = faceElem.QPoint(q);
......@@ -56,7 +56,7 @@ void DGLinearTransportAssemble::SystemMatrix(Matrix &systemMatrix) const {
} else {
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);
DGFaceElement felem_1(systemMatrix, cf, f1);
DGRowEntries A_cf(systemMatrix, c, cf);
for (int q = 0; q < faceElem.nQ(); ++q) {
const Point &Qf_c = faceElem.QPoint(q);
......@@ -137,7 +137,7 @@ void DGLinearTransportAssemble::RHS(double t, Vector &rhs) const {
row r = rhs.find_row(c());
for (int f = 0; f < c.Faces(); ++f) {
if (!rhs.GetMesh().onBndDG(c, f)) continue;
DGFaceElement felem(*disc, rhs, c, f);
DGFaceElement felem(rhs, c, f);
for (int q = 0; q < felem.nQ(); ++q) {
const Point &z = felem.QPoint(q);
double w = felem.QWeight(q);
......@@ -157,7 +157,7 @@ void DGLinearTransportAssemble::RHS(double t, Vector &rhs) const {
double DGLinearTransportAssemble::Energy(const Vector &u) const {
double energy = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
......@@ -170,7 +170,7 @@ double DGLinearTransportAssemble::Energy(const Vector &u) const {
void DGNonLinearTransportAssemble::Energy(const cell &c,
const Vector &u,
double &energy) const {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
......@@ -181,7 +181,7 @@ void DGNonLinearTransportAssemble::Energy(const cell &c,
double DGLinearTransportAssemble::Error(double t, const Vector &u) const {
double err = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
......@@ -198,7 +198,7 @@ std::pair<double, double> DGLinearTransportAssemble::InflowOutflow(const Vector
for (cell c = u.cells(); c != u.cells_end(); ++c) {
BFParts bnd(u.GetMesh(), c);
if (!bnd.onBnd()) continue;
DGElement elem(*disc, u, c);
DGElement elem(u, c);
double U = 0;
double area = 0;
for (int q = 0; q < elem.nQ(); ++q) {
......@@ -209,7 +209,7 @@ std::pair<double, double> DGLinearTransportAssemble::InflowOutflow(const Vector
U *= (1 / area);
for (int f = 0; f < c.Faces(); ++f) {
if (bnd[f] == -1) continue;
DGFaceElement faceElem(*disc, u, c, f);
DGFaceElement faceElem(u, c, f);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
VectorField Nq = faceElem.QNormal(q);
......@@ -232,7 +232,7 @@ void DGLinearTransportAssemble::SetExactSolution(double t, Vector &u_ex) const {
u_ex = 0;
for (cell c = u_ex.cells(); c != u_ex.cells_end(); ++c) {
row r = u_ex.find_row(c());
DGElement Elem(*disc, u_ex, c);
DGElement Elem(u_ex, c);
for (int j = 0; j < Elem.NodalPoints(); ++j)
u_ex(r, j) = problem->Solution(t, Elem.NodalPoint(j));
}
......@@ -242,7 +242,7 @@ void DGLinearTransportAssemble::SetExactSolution(double t, Vector &u_ex) const {
double DGLinearTransportAssemble::Mass(const Vector &u) const {
Scalar e = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
......@@ -256,7 +256,7 @@ void DGLinearTransportAssemble::SetInitialValue(Vector &u) const {
u = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
row r = u.find_row(c());
DGElement Elem(*disc, u, c);
DGElement Elem(u, c);
double lambda = 1e-9;
for (int j = 0; j < Elem.NodalPoints(); ++j)
u(r, j) = problem->Solution(0, lambda * c() + (1 - lambda) * Elem.NodalPoint(j));
......@@ -266,7 +266,7 @@ void DGLinearTransportAssemble::SetInitialValue(Vector &u) const {
void DGLinearTransportAssemble::PrintMatrixInfo(Matrix &A, int diagonal) const {
for (cell c = A.cells(); c != A.cells_end(); ++c) {
DGElement elem(*disc, A, c);
DGElement elem(A, c);
DGRowEntries A_c(A, c, c);
for (int i = 0; i < elem.NodalPoints(); ++i) {
for (int j = 0; j < elem.NodalPoints(); ++j)
......@@ -299,7 +299,7 @@ void DGLinearTransportAssemble::VtkPlotting_cell(double t,
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);
DGElement elem(u_tmp, c);
Scalar U = 0.0;
double a = 0;
for (int q = 0; q < elem.nQ(); ++q) {
......
......@@ -11,7 +11,7 @@ void DGEllipticAssemble::Initialize(Vector &u) const {
// for (cell c = u.cells(); c != u.cells_end(); ++c) {
// RowBndValues u_c(u, c);
// if (!u_c.onBnd()) continue;
// DGElement elem(*disc, u, c);
// DGElement elem(u, c);
// for (int face = 0; face < c.Faces(); ++face) {
// if (u_c.bc(face) == 1) {
// for (int j = 0; j < u.NodalPointsOnFace(c, face); ++j) {
......@@ -28,7 +28,7 @@ void DGEllipticAssemble::Initialize(Vector &u) const {
double DGEllipticAssemble::Energy(const Vector &u) const {
double energy = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
VectorField DU = elem.Derivative(q, u);
......@@ -41,7 +41,7 @@ double DGEllipticAssemble::Energy(const Vector &u) const {
}
void DGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) const {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
RowBndValues r_c(r, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
......@@ -57,7 +57,7 @@ void DGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) con
}
for (int f = 0; f < c.Faces(); ++f) {
Scalar scaledPenalty = penalty / c.FaceArea(f);
DGFaceElement faceElem(*disc, u, c, f);
DGFaceElement faceElem(u, c, f);
if (r.GetMesh().onBndDG(c, f)) {
int bnd = r_c.bc(f);
if (bnd == 2) {
......@@ -91,7 +91,7 @@ void DGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) con
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 otherFaceElem(*disc, r, cf, f1);
DGFaceElement otherFaceElem(r, cf, f1);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
const Point &N = faceElem.QNormal(q);
......@@ -125,7 +125,7 @@ void DGEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r) con
}
void DGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const {
DGElement elem(*disc, A, c);
DGElement elem(A, c);
DGRowEntries A_c(A, c, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
......@@ -142,7 +142,7 @@ void DGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const
BFParts bnd(u.GetMesh(), c);
for (int f = 0; f < c.Faces(); ++f) {
Scalar s = penalty / c.FaceArea(f);
DGFaceElement faceElem(*disc, u, c, f);
DGFaceElement faceElem(u, c, f);
if (u.GetMesh().onBndDG(c, f)) {
if (bnd[f] == 1) {
for (int q = 0; q < faceElem.nQ(); ++q) {
......@@ -169,7 +169,7 @@ void DGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const
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 otherFaceElem(*disc, u, cf, f1);
DGFaceElement otherFaceElem(u, cf, f1);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
const Point &N = faceElem.QNormal(q);
......@@ -217,7 +217,7 @@ void DGEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &A) const
double DGEllipticAssemble::EnergyError(const Vector &u) const {
double err = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Point z = elem.QPoint(q);
......@@ -234,7 +234,7 @@ double DGEllipticAssemble::EnergyError(const Vector &u) const {
double DGEllipticAssemble::L2(const Vector &u) const {
double l2 = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
......@@ -247,7 +247,7 @@ double DGEllipticAssemble::L2(const Vector &u) const {
double DGEllipticAssemble::H1(const Vector &u) const {
double err = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
......@@ -261,7 +261,7 @@ double DGEllipticAssemble::H1(const Vector &u) const {
double DGEllipticAssemble::L2Error(const Vector &u) const {
double err = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Scalar U = elem.Value(q, u);
......@@ -275,7 +275,7 @@ double DGEllipticAssemble::L2Error(const Vector &u) const {
double DGEllipticAssemble::L2CellAverageError(const Vector &u) const {
double err = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
double w = 0;
Scalar U = 0;
Scalar Sol = 0;
......@@ -292,7 +292,7 @@ double DGEllipticAssemble::L2CellAverageError(const Vector &u) const {
double DGEllipticAssemble::MaxError(const Vector &u) const {
double err = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
Scalar p = elem.Value(q, u);
err = max(err, abs(p - problem->Solution(elem.QPoint(q))));
......@@ -308,7 +308,7 @@ double DGEllipticAssemble::FluxError(const Vector &u) const {
if (!bnd.onBnd()) continue;
for (int face = 0; face < c.Faces(); ++face) {
if (bnd[face] == 2) {
DGFaceElement faceElem(*disc, u, c, face);
DGFaceElement faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
VectorField G = problem->Flux(faceElem.QPoint(q));
......@@ -318,7 +318,7 @@ double DGEllipticAssemble::FluxError(const Vector &u) const {
flux_error += w * F * F;
}
} else if (bnd[face] == 0) {
DGFaceElement faceElem(*disc, u, c, face);
DGFaceElement faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Tensor K = problem->Permeability(c);
......@@ -336,7 +336,7 @@ double DGEllipticAssemble::FaceError(const Vector &u) const {
double face_error = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
for (int face = 0; face < c.Faces(); ++face) {
DGFaceElement faceElem(*disc, u, c, face);
DGFaceElement faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar U = faceElem.Value(q, u);
......@@ -356,7 +356,7 @@ FluxPair DGEllipticAssemble::InflowOutflow(const Vector &u) const {
if (!bnd.onBnd()) continue;
for (int face = 0; face < c.Faces(); ++face) {
if (bnd[face] < 0) continue;
DGFaceElement faceElem(*disc, u, c, face);
DGFaceElement faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Tensor K = problem->Permeability(c);
......@@ -379,7 +379,7 @@ FluxPair DGEllipticAssemble::PrescribedInflowOutflow(const Vector &u) const {
if (!bnd.onBnd()) continue;
for (int face = 0; face < c.Faces(); ++face) {
if (bnd[face] < 2) continue;
DGFaceElement faceElem(*disc, u, c, face);
DGFaceElement faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar FN = problem->Flux(faceElem.QPoint(q)) *
......@@ -401,7 +401,7 @@ FluxPair DGEllipticAssemble::OutflowLeftRight(const Vector &u) const {
if (!bnd.onBnd()) continue;
for (int face = 0; face < c.Faces(); ++face) {
if (bnd[face] != 1) continue;
DGFaceElement faceElem(*disc, u, c, face);
DGFaceElement faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Tensor K = problem->Permeability(c);
......@@ -423,7 +423,7 @@ double DGEllipticAssemble::GoalFunctional(const Vector &u) const {
if (!bnd.onBnd()) continue;
for (int face = 0; face < c.Faces(); ++face) {
if (bnd[face] != 1) continue;
DGFaceElement faceElem(*disc, u, c, face);
DGFaceElement faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
Point z = faceElem.QPoint(q);
if (z[0] < 0.25) continue;
......@@ -443,7 +443,7 @@ void DGEllipticAssemble::SetExactSolution(Vector &u) const {
u = 0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
row r = u.find_row(c());
DGElement elem(*disc, u, c);
DGElement elem(u, c);
for (int j = 0; j < elem.NodalPoints(); ++j)
u(r, j) = problem->Solution(elem.NodalPoint(j));
}
......
......@@ -37,7 +37,7 @@ void HybridEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r)
RVector RU(R);
RVector UU(R);
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
for (int face = 0; face < Faces; ++face)
UU[face] = u(elem[face], 0);
for (int face = 0; face < Faces; ++face)
......@@ -56,7 +56,7 @@ void HybridEllipticAssemble::Residual(const cell &c, const Vector &u, Vector &r)
for (int face = 0; face < Faces; ++face) {
int bnd = r_c.bc(face);
if (bnd != 2) continue;
RTFaceElementT<> faceElem(*disc, u, c, face);
RTFaceElementT<> faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar F = problem->Flux(faceElem.QPoint(q)) * faceElem.QNormal(q);
......@@ -74,7 +74,7 @@ void HybridEllipticAssemble::Jacobi(const cell &c, const Vector &u, Matrix &J) c
RVector IAB = IA * B;
Scalar BIAB = B * IAB;
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
RowEntries J_c(J, elem);
for (int face_i = 0; face_i < Faces; ++face_i)
for (int face_j = 0; face_j < Faces; ++face_j)
......@@ -86,7 +86,7 @@ void HybridEllipticAssemble::LocalProblem(const cell &c, const Vector &u,
RMatrix &A,
RVector &B, RVector &R) const {
A = 0, B = 0, R = 0;
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Tensor IK = Invert(problem->Permeability(c));
......@@ -103,7 +103,7 @@ void HybridEllipticAssemble::LocalProblem(const cell &c, const Vector &u,
}
}
for (int face = 0; face < c.Faces(); ++face) {
RTFaceElementT<> faceElem(*disc, u, c, face);
RTFaceElementT<> faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
for (int i = 0; i < faceElem.size(); i++) {
......@@ -122,7 +122,7 @@ std::pair<double, double> HybridEllipticAssemble::InflowOutflow(const Vector &u)
if (!bnd.onBnd()) continue;
for (int face = 0; face < c.Faces(); ++face) {
if (bnd[face] < 0) continue;
RTFaceElementT<> faceElem(*disc, u, c, face);
RTFaceElementT<> faceElem(u, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
double w = faceElem.QWeight(q);
Scalar F = FaceFlux(c, u, faceElem, q, face) *
......@@ -144,7 +144,7 @@ FluxPair HybridEllipticAssemble::OutflowLeftRight(const Vector &p) const {
if (!bnd.onBnd()) continue;
for (int face = 0; face < c.Faces(); ++face) {
if (bnd[face] < 0) continue;
RTFaceElementT<> faceElem(*disc, p, c, face);
RTFaceElementT<> faceElem(p, c, face);
for (int q = 0; q < faceElem.nQ(); ++q) {
Point z = faceElem.QPoint(q);
double w = faceElem.QWeight(q);
......@@ -171,7 +171,7 @@ void HybridEllipticAssemble::SetPressure(const Vector &u, Vector &p) const {
Scalar BIAB = B * IAB;
RVector RU(R);
RVector UU(R);
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
for (int i = 0; i < c.Faces(); ++i)
UU[i] = u(elem[i], 0);
for (int i = 0; i < c.Faces(); ++i)
......@@ -193,7 +193,7 @@ void HybridEllipticAssemble::SetFlux(const Vector &u, Vector &flux) {
Scalar BIAB = B * IAB;
RVector RU(R);
RVector UU(R);
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
for (int i = 0; i < c.Faces(); ++i)
UU[i] = u(elem[i], 0);
for (int i = 0; i < c.Faces(); ++i)
......@@ -239,7 +239,7 @@ void HybridEllipticAssemble::SetNormalFlux(const Vector &u, Vector &flux) {
Scalar BIAB = B * IAB;
RVector RU(R);
RVector UU(R);
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
for (int face = 0; face < c.Faces(); ++face)
UU[face] = u(elem[face], 0);
for (int face = 0; face < c.Faces(); ++face)
......@@ -259,7 +259,7 @@ VectorField HybridEllipticAssemble::EvaluateCellFlux(const Vector &flux,
const cell &c) const {
const Mesh &M = flux.GetMesh();
if (M.find_cell(c()) == M.cells_end()) return zero;
RTLagrangeElement elem(*disc, flux, c);
RTLagrangeElement elem(flux, c);
VectorField F = zero;
double area = 0;
for (int q = 0; q < elem.nQ(); ++q) {
......@@ -275,14 +275,14 @@ double HybridEllipticAssemble::EvaluateNormalFlux(const Vector &flux, const cell
int i) const {
const Mesh &M = flux.GetMesh();
if (M.find_cell(c()) == M.cells_end()) return 0;
RTLagrangeElement elem(*disc, flux, c);
RTFaceElementT<> faceElem(*disc, flux, c, i);
RTLagrangeElement elem(flux, c);
RTFaceElementT<> faceElem(flux, c, i);
return flux(elem[i], 0) * elem.Sign(i) / faceElem.Area();
}
VectorField HybridEllipticAssemble::EvaluateNormal(const Vector &flux,
const cell &c, int i) const {
RTLagrangeElement elem(*disc, flux, c);
RTLagrangeElement elem(flux, c);
return elem.Normal(i);
}
......@@ -297,7 +297,7 @@ Scalar HybridEllipticAssemble::Value(const cell &c, const Vector &u,
Scalar BIAB = B * IAB;
RVector RU(R);
RVector UU(R);
RTLagrangeElement E(*disc, u, c);
RTLagrangeElement E(u, c);
for (int face = 0; face < c.Faces(); ++face)
UU[face] = u(E[face], 0);
for (int face = 0; face < c.Faces(); ++face)
......@@ -318,7 +318,7 @@ Scalar HybridEllipticAssemble::FaceValue(const cell &c, const Vector &u,
Scalar BIAB = B * IAB;
RVector RU(R);
RVector UU(R);
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
for (int face = 0; face < c.Faces(); ++face)
UU[face] = u(elem[face], 0);
for (int face = 0; face < Faces; ++face)
......@@ -360,7 +360,7 @@ VectorField HybridEllipticAssemble::Flux(const cell &c, const Vector &u,
VectorField HybridEllipticAssemble::FaceFlux(const cell &c, const Vector &u,
RTFaceElementT<> &faceElem,
int q, int face) const {
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
int Faces = c.Faces();
RMatrix A(Faces);
RVector B(Faces), R(Faces);
......@@ -386,7 +386,7 @@ VectorField HybridEllipticAssemble::FaceFlux(const cell &c, const Vector &u,
double HybridEllipticAssemble::EnergyError(const Vector &u) const {
double err = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);
for (int q = 0; q < elem.nQ(); ++q) {
double w = elem.QWeight(q);
Tensor IK = Invert(problem->Permeability(c));
......@@ -401,7 +401,7 @@ double HybridEllipticAssemble::EnergyError(const Vector &u) const {
double HybridEllipticAssemble::L2(const Vector &u) const {
double l2 = 0.0;
for (cell c = u.cells(); c != u.cells_end(); ++c) {
RTLagrangeElement elem(*disc, u, c);
RTLagrangeElement elem(u, c);