Commit 3af3ab6b authored by niklas.baumgarten's avatar niklas.baumgarten
Browse files

[Add-GetCell-To-Elements] added GetCell to elements

parent a8590fb7
......@@ -47,6 +47,8 @@ public:
int Bnd(int i) const { return bnd[i]; }
const cell &GetCell() const { return c; }
virtual int nQ() const { return Q.size(); }
virtual TT QWeight(int q) const { return qWeight[q]; }
......
......@@ -4,6 +4,7 @@
#include "IDiscretization.hpp"
#include "Algebra.hpp"
#ifdef BUILD_IA
#include "IAInterval.hpp"
......@@ -17,7 +18,8 @@ class BasicFaceElementT : public rows {
for (int q = 0; q < Q.size(); ++q)
refArea += Q.Weight(q);
if constexpr (!std::is_same_v<TT, double>) {
WarningOnMaster("LocalFaceArea only used in double precision: -> implement for IAInterval")
WarningOnMaster(
"LocalFaceArea only used in double precision: -> implement for IAInterval")
}
refArea = c.LocalFaceArea(faceID) / refArea;
......@@ -43,7 +45,7 @@ class BasicFaceElementT : public rows {
void computeQPointCellLocal(int q) {
const PointT<TT, sDim, tDim> &x = Q.QPoint(q);
switch (FaceCellType(c.Type())) {
switch(FaceCellType(c.Type())) {
case POINT: {
qPointCellLocal[q] = c.LocalCorner(c.facecorner(faceID, 0));
break;
......@@ -93,47 +95,47 @@ class BasicFaceElementT : public rows {
if ((y0 < y1) && (y0 < y2) && (y0 < y3)) {
if (y1 < y3)
qPointCellLocal[q] = (1 - x[0]) * (1 - x[1]) * z0
+ x[0] * (1 - x[1]) * z1
+ x[0] * x[1] * z2
+ (1 - x[0]) * x[1] * z3;
+ x[0] * (1 - x[1]) * z1
+ x[0] * x[1] * z2
+ (1 - x[0]) * x[1] * z3;
else
qPointCellLocal[q] = (1 - x[0]) * (1 - x[1]) * z0
+ x[0] * (1 - x[1]) * z3
+ x[0] * x[1] * z2
+ (1 - x[0]) * x[1] * z1;
+ x[0] * (1 - x[1]) * z3
+ x[0] * x[1] * z2
+ (1 - x[0]) * x[1] * z1;
} else if ((y1 < y0) && (y1 < y2) && (y1 < y3)) {
if (y2 < y0)
qPointCellLocal[q] = (1 - x[0]) * (1 - x[1]) * z1
+ x[0] * (1 - x[1]) * z2
+ x[0] * x[1] * z3
+ (1 - x[0]) * x[1] * z0;
+ x[0] * (1 - x[1]) * z2
+ x[0] * x[1] * z3
+ (1 - x[0]) * x[1] * z0;
else
qPointCellLocal[q] = (1 - x[0]) * (1 - x[1]) * z1
+ x[0] * (1 - x[1]) * z0
+ x[0] * x[1] * z3
+ (1 - x[0]) * x[1] * z2;
+ x[0] * (1 - x[1]) * z0
+ x[0] * x[1] * z3
+ (1 - x[0]) * x[1] * z2;
} else if ((y2 < y0) && (y2 < y1) && (y2 < y3)) {
if (y3 < y1)
qPointCellLocal[q] = (1 - x[0]) * (1 - x[1]) * z2
+ x[0] * (1 - x[1]) * z3
+ x[0] * x[1] * z0
+ (1 - x[0]) * x[1] * z1;
+ x[0] * (1 - x[1]) * z3
+ x[0] * x[1] * z0
+ (1 - x[0]) * x[1] * z1;
else
qPointCellLocal[q] = (1 - x[0]) * (1 - x[1]) * z2
+ x[0] * (1 - x[1]) * z1
+ x[0] * x[1] * z0
+ (1 - x[0]) * x[1] * z3;
+ x[0] * (1 - x[1]) * z1
+ x[0] * x[1] * z0
+ (1 - x[0]) * x[1] * z3;
} else if ((y3 < y0) && (y3 < y2) && (y3 < y1)) {
if (y0 < y2)
qPointCellLocal[q] = (1 - x[0]) * (1 - x[1]) * z3
+ x[0] * (1 - x[1]) * z0
+ x[0] * x[1] * z1
+ (1 - x[0]) * x[1] * z2;
+ x[0] * (1 - x[1]) * z0
+ x[0] * x[1] * z1
+ (1 - x[0]) * x[1] * z2;
else
qPointCellLocal[q] = (1 - x[0]) * (1 - x[1]) * z3
+ x[0] * (1 - x[1]) * z2
+ x[0] * x[1] * z1
+ (1 - x[0]) * x[1] * z0;
+ x[0] * (1 - x[1]) * z2
+ x[0] * x[1] * z1
+ (1 - x[0]) * x[1] * z0;
}
break;
}
......@@ -171,10 +173,10 @@ protected:
public:
BasicFaceElementT(const VectorMatrixBase &base, const cell &c, int face, int n = 0)
: faceID(face), rows(base.GetMatrixGraph(), c, face, n),
Q(base.GetDiscT<TT, sDim, tDim>().GetFaceQuad(c)), c(c), qWeight(this->nQ()),
qPoint(this->nQ()), qNormal(this->nQ()), qPointCellLocal(this->nQ()),
transformation(1) {
: faceID(face), rows(base.GetMatrixGraph(), c, face, n),
Q(base.GetDiscT<TT, sDim, tDim>().GetFaceQuad(c)), c(c), qWeight(this->nQ()),
qPoint(this->nQ()), qNormal(this->nQ()), qPointCellLocal(this->nQ()),
transformation(1) {
init();
}
......@@ -182,6 +184,8 @@ public:
return area;
}
const cell &GetCell() const { return c; }
virtual int nQ() const { return Q.size(); }
virtual int NodalPoints() const { return this->size(); }
......@@ -234,8 +238,8 @@ protected:
const ShapeT<TT, sDim, tDim> &shape;
public:
FaceElementT(const VectorMatrixBase &base, const cell &c, int face, int n = 0)
: BasicFaceElementT<TT, sDim, tDim>(base, c, face, n),
shape(base.GetDiscT<TT, sDim, tDim>().GetShape(c, n)) {}
: BasicFaceElementT<TT, sDim, tDim>(base, c, face, n),
shape(base.GetDiscT<TT, sDim, tDim>().GetShape(c, n)) {}
};
typedef FaceElementT<> FaceElement;
......@@ -248,10 +252,11 @@ protected:
vector<const ShapeT<TT, sDim, tDim> *> shapes;
public:
MixedFaceElementT(const VectorMatrixBase &base, const cell &c, int face, int n = 0)
: BasicFaceElementT<TT, sDim, tDim>(base, c, face, n),
nodalPointData(base.GetMatrixGraph().GetShapeDoFs().NodalPointData(c)),
nodalPointFaceData_face(base.GetMatrixGraph().GetShapeDoFs().NodalPointFaceData(c, face)),
shapes(base.GetMatrixGraph().GetShapeDoFs().NumberOfDoFs()) {
: BasicFaceElementT<TT, sDim, tDim>(base, c, face, n),
nodalPointData(base.GetMatrixGraph().GetShapeDoFs().NodalPointData(c)),
nodalPointFaceData_face(base.GetMatrixGraph().GetShapeDoFs()
.NodalPointFaceData(c, face)),
shapes(base.GetMatrixGraph().GetShapeDoFs().NumberOfDoFs()) {
for (int m = 0; m < shapes.size(); ++m)
shapes[m] = &base.GetDiscT<TT, sDim, tDim>().GetShape(c, m);
}
......
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