Skip to content

GitLab

Commit fc85251a authored by daniele.corallo's avatar daniele.corallo
Browse files

[402-move-stelements-to-stproject] rmv ST-Elements

parent a8590fb7
Hide whitespace changes
Inline Side-by-side
src/lib4_fem/CMakeLists.txt
View file @ fc85251a
......@@ -51,8 +51,6 @@ endif()
if(USE_SPACETIME)
set (fem_src ${fem_src}
elements/STDGDGViscoAcousticElement.cpp
elements/STDGDGViscoAcousticFaceElement.cpp
discretization/SpaceTimeDiscretization.cpp
discretization/SpaceTimeDiscretizations.cpp
discretization/SpaceTimeDiscretization_PGDG.cpp
......
src/lib4_fem/elements/STDGDGViscoAcousticElement.cpp deleted 100644 → 0
View file @ a8590fb7
#include "STDGDGViscoAcousticElement.hpp"
const SpaceTimeShape & GetSpaceTimeShape(const VectorMatrixBase &base, cell c){
auto &disc = dynamic_cast<const STDiscretization &>(base.GetDisc());
int sd = base.GetDoFs().get_space_deg(c);
int td = base.GetDoFs().get_time_deg(c);
return disc.GetSpaceTimeShape(sd, td);
}
const SpaceTimeShape & GetSpaceTimeShapeHighOrder(const VectorMatrixBase &base, cell c){
auto &disc = dynamic_cast<const STDiscretization &>(base.GetDisc());
return disc.GetSpaceTimeShape(6, 6);
}
const SpaceTimeQuadrature & GetSpaceTimeQuadrature(const VectorMatrixBase &base, cell c){
auto &disc = dynamic_cast<const STDiscretization &>(base.GetDisc());
int sd = base.GetDoFs().get_space_deg(c);
int td = base.GetDoFs().get_time_deg(c);
return disc.GetSpaceTimeQuad(sd, td);
}
const SpaceTimeQuadrature & GetSpaceTimeQuadratureHighOrder(const VectorMatrixBase &base){
auto &disc = dynamic_cast<const STDiscretization &>(base.GetDisc());
return disc.GetSpaceTimeQuad(6, 6);
}
STDGDGViscoAcousticElement::STDGDGViscoAcousticElement
(const VectorMatrixBase &g, const cell &tc, int nL, bool max_quad_order)
:
shape(GetSpaceTimeShape(g, tc)),
quad(GetSpaceTimeQuadrature(g, tc)),
TC(tc), r(g.find_row(tc())), dim(tc.dim()),
gradients(quad.size(), shape.size()),
numL(nL),
space_quad(max_quad_order ? GetSpaceTimeQuadratureHighOrder(g).GetSpaceQuad() : quad.GetSpaceQuad()){
qWeight.resize(nQ());
qPoint.resize(nQ());
velocityST.resize(nQ() * j_dimension());
dtVelocityST.resize(nQ() * j_dimension());
divVelocityST.resize(nQ() * j_dimension());
pressureST.resize(nQ() * j_dimension());
dtPressureST.resize(nQ() * j_dimension());
gradPressureST.resize(nQ() * j_dimension());
const TransformationT<double, SpaceDimension, 1> &T = TC.GetSpaceTimeTransformation();
det_time = T.DetTime();
det_space = T.DetSpace();
for (int q = 0; q < quad.size(); q++) {
qWeight[q] = T.Det() * quad.Weight(q);
qPoint[q] = TC.LocalToGlobal(quad.QPoint(q));
for (int i = 0; i < shape.size(); i++) {
VectorFieldT<double, SpaceDimension, 1> STGrad = shape.LocalSpaceTimeGradient(q, i);
gradients[q][i] = T * STGrad;
}
}
for (int tq = 0; tq < quad.GetTimeQuad().size(); ++tq) {
for (int sq = 0; sq < quad.GetSpaceQuad().size(); ++sq) {
int quad_i = tq * quad.GetSpaceQuad().size() + sq;
for (int ti = 0; ti < shape.GetTimeShape().size(); ++ti) {
for (int si = 0; si < shape.GetSpaceShape().size(); ++si) {
int shape_i = ti * shape.GetSpaceShape().size() + si;
for (int k = 0; k < dim; ++k) {
int index = (tq * quad.GetSpaceQuad().size() + sq) * j_dimension()
+ ti * shape.GetSpaceShape().size() * (1 + numL + dim)
+ k * shape.GetSpaceShape().size() + si;
velocityST[index][k] = shape(quad_i, shape_i);
dtVelocityST[index][k] = gradients[quad_i][shape_i].t();
divVelocityST[index] = gradients[quad_i][shape_i][k];
}
for (int k = 0; k <= numL; ++k) {
int index = (tq * quad.GetSpaceQuad().size() + sq) * j_dimension()
+ ti * shape.GetSpaceShape().size() * (1 + numL + dim)
+ (dim + k) * shape.GetSpaceShape().size()
+ si;
pressureST[index] = shape.GetSpaceShape()(sq, si) * shape.GetTimeShape()(tq, ti);
dtPressureST[index] = gradients[quad_i][shape_i].t();
gradPressureST[index] = gradients[quad_i][shape_i].sliceTime();
}
}
}
}
}
}
src/lib4_fem/elements/STDGDGViscoAcousticElement.hpp deleted 100644 → 0
View file @ a8590fb7
#ifndef SPACETIME_STDGDGVISCOACOUSTICELEMENT_HPP
#define SPACETIME_STDGDGVISCOACOUSTICELEMENT_HPP
#include "SpaceTimeDiscretization.hpp"
#include "Element.hpp"
#include "SpaceTimeQuadrature.hpp"
#include "SpaceTimeShape.hpp"
class STDGDGViscoAcousticElement {
public:
double det_space;
double det_time;
vector<double> space_qWeight;
vector<double> time_qWeight;
vector<double> qWeight;
vector<Point> qPoint;
int dim;
const SpaceTimeShape &shape;
const SpaceTimeQuadrature &quad;
const Quadrature &space_quad;
const cell &TC;
const row r;
ShapeValues<VectorFieldT<double, SpaceDimension, 1>> gradients;
vector<VectorField> velocityST;
vector<VectorField> dtVelocityST;
vector<double> divVelocityST;
vector<double> pressureST;
vector<double> dtPressureST;
vector<VectorField> gradPressureST;
int numL;
STDGDGViscoAcousticElement(const VectorMatrixBase &,
const cell &,
int nL = 0,
bool max_quad_order = false);
int nQ() const { return quad.size(); }
int nSpaceQ() const { return space_quad.size(); }
int nTimeQ() const { return quad.GetTimeQuad().size(); }
double QWeight(int q) const { return qWeight[q]; }
double SpaceQWeight(int qq) {
if (space_qWeight.empty()) {
space_qWeight.resize(space_quad.size());
for (int q = 0; q < space_qWeight.size(); q++) {
space_qWeight[q] = det_space * space_quad.Weight(q);
}
}
return space_qWeight[qq];
}
Point SpaceQPoint(int q) const {
return space_quad.QPoint(q);
}
double TimeQWeight(int qq) {
if (time_qWeight.empty()) {
time_qWeight.resize(quad.GetTimeQuad().size());
for (int q = 0; q < time_qWeight.size(); q++) {
time_qWeight[q] = det_time * quad.GetTimeQuad().Weight(q);
}
}
return time_qWeight[qq];
}
const Point &QPoint(int q) const { return qPoint[q]; }
Point LocalToGlobal(const Point &local) const;
Point GlobalToLocal(const Point &) const;
double Area() const {
double a = 0;
for (int q = 0; q < nQ(); ++q) a += qWeight[q];
return a;
}
int dimension() const {Exit("Not implemented"); }
int i_dimension() const { return shape.size() * (1 + dim + numL); }
int j_dimension() const { return shape.size() * (1 + dim + numL); }
int variable(int j) {
int ss_size = shape.GetSpaceShape().size();
int tmp = j % (ss_size * (dim + 1 + numL));
if (tmp < dim * ss_size) return 0;
for (int nL = 0; nL <= numL; nL++)
if (tmp - ss_size * dim < (nL + 1) * ss_size) return 1 + nL;
Exit("BLUPP");
}
VectorField Velocity(int nq, int j) {
return velocityST[nq * j_dimension() + j];
}
VectorField VelocityLocal(const Point localP, int j) {
int shape_idx = j % shape.size();
int component = j / shape.size();
if (component == 0) {
return VectorField(shape(localP, shape_idx), 0.0);
} else if (component == 1) {
return VectorField(0.0, shape(localP, shape_idx));
} else {
return VectorField();
}
}
VectorField Velocity(const Point p, int j) {
for (int tq = 0; tq < quad.GetTimeQuad().size(); ++tq) {
for (int sq = 0; sq < quad.GetSpaceQuad().size(); ++sq) {
for (int ti = 0; ti < shape.GetTimeShape().size(); ++ti) {
for (int si = 0; si < shape.GetSpaceShape().size(); ++si) {
int shape_i = ti * shape.GetSpaceShape().size() + si;
for (int k = 0; k < dim; ++k) {
int index = (tq * quad.GetSpaceQuad().size() + sq) * j_dimension()
+ ti * shape.GetSpaceShape().size() * (1 + numL + dim)
+ k * shape.GetSpaceShape().size() + si;
if (index == j) {
double value = shape(p, shape_i);
if (k == 0)
return VectorField(value, 0.0);
else
return VectorField(0.0, value);
}
}
}
}
}
}
return zero;
// int shape_idx = j % shape.size();
// int component = j / shape.size();
// if (component == 0) {
// return VectorField(shape(p, shape_idx), 0.0);
// } else if (component == 1) {
// return VectorField(0.0, shape(p, shape_idx));
// } else {
// return zero;
// }
}
VectorField Velocity(Point P, const double *u) {
VectorField V = zero;
for (int j = 0; j < j_dimension(); ++j)
V += u[j] * Velocity(P, j);
return V;
}
VectorField VelocityGlobal(const Point &globalPoint, const Vector &u) {
Point localPoint = TC.GlobalToLocal(globalPoint);
VectorField V = zero;
for (int j = 0; j < j_dimension(); ++j)
V += u(r, j) * Velocity(localPoint, j);
return V;
}
[[deprecated]] VectorField Velocity(const Point &globalPoint, const Vector &u) {
Point localPoint = TC.GlobalToLocal(globalPoint);
VectorField V = zero;
for (int j = 0; j < j_dimension(); ++j)
V += u(r, j) * Velocity(localPoint, j);
return V;
}
VectorField Velocity(int nq, const double *u) {
VectorField V = zero;
for (int j = 0; j < j_dimension(); ++j)
V += u[j] * Velocity(nq, j);
return V;
}
VectorField DtVelocity(int nq, int j) {
return dtVelocityST[nq * j_dimension() + j];
}
double DivVelocity(int nq, int j) {
return divVelocityST[nq * j_dimension() + j];
}
double Pressure(int nq, int j) {
return pressureST[nq * j_dimension() + j];
}
double Pressure(const Point p, int j) {
for (int tq = 0; tq < quad.GetTimeQuad().size(); ++tq) {
for (int sq = 0; sq < quad.GetSpaceQuad().size(); ++sq) {
for (int ti = 0; ti < shape.GetTimeShape().size(); ++ti) {
for (int si = 0; si < shape.GetSpaceShape().size(); ++si) {
int shape_i = ti * shape.GetSpaceShape().size() + si;
int index = (tq * quad.GetSpaceQuad().size() + sq) * j_dimension()
+ ti * shape.GetSpaceShape().size() * (1 + numL + dim)
+ dim * shape.GetSpaceShape().size() + si;
if (index == j)
return shape(p, shape_i);
}
}
}
}
return 0.0;
// int shape_idx = j / (1 + dim + numL);
// int component = j % (1 + dim + numL);
// if (component > 1) {
// return shape(p, shape_idx);
// }
// return 0.0;
}
double Pressure(const Point p, const double *u) {
double P = 0;
for (int j = 0; j < j_dimension(); ++j)
P += u[j] * Pressure(p, j);
return P;
}
double PressureGlobal(const Point &globalPoint, const Vector &u) {
Point localPoint = TC.GlobalToLocal(globalPoint);
double P = 0;
for (int j = 0; j < j_dimension(); ++j)
P += u(r, j) * Pressure(localPoint, j);
return P;
}
[[deprecated]] double Pressure(const Point &globalPoint, const Vector &u) {
Point localPoint = TC.GlobalToLocal(globalPoint);
double P = 0;
for (int j = 0; j < j_dimension(); ++j)
P += u(r, j) * Pressure(localPoint, j);
return P;
}
double Pressure(int nq, const double *u) {
double P = 0;
for (int j = 0; j < j_dimension(); ++j)
P += u[j] * Pressure(nq, j);
return P;
}
double DtPressure(int nq, int j) {
return dtPressureST[nq * j_dimension() + j];
}
VectorField GradPressure(int nq, int j) {
return gradPressureST[nq * j_dimension() + j];
}
};
using SpaceTimeViscoAcousticDGTElement = STDGDGViscoAcousticElement;
#endif //SPACETIME_STDGDGVISCOACOUSTICELEMENT_HPP
src/lib4_fem/elements/STDGDGViscoAcousticFaceElement.cpp deleted 100644 → 0
View file @ a8590fb7
#include "STDGDGViscoAcousticFaceElement.hpp"
const SpaceTimeShape &GetSpaceTimeShape(const VectorMatrixBase &base, cell c);
SpaceTimeViscoAcousticDGTFaceElement::SpaceTimeViscoAcousticDGTFaceElement
(const STDiscretization &D, const VectorMatrixBase &g, const cell &tc, int f_id, int nL)
:
shape(GetSpaceTimeShape(g, tc)),
faceqpoints(D.FaceQPoints()),
TC(tc), f(g.GetMesh().find_face(tc.Face(f_id))),
fid(f_id), dim(tc.dim()),
SpaceQ(D.GetCellFaceQuad(f_id, g.GetDoFs().get_space_deg(tc))),
TimeQ(D.GetTimeFaceQuad(f_id, g.GetDoFs().get_time_deg(tc))),
numL(nL) {
double ww = TC.SpaceCell().LocalFaceArea(fid);
double a = TC.min();
double b = TC.max();
det_time = b - a;
qWeight.resize(nQ());
qPoint.resize(nQ());
qLocal.resize(nQ());
qNormal.resize(nQ());
velocityST.resize(nQ() * j_dimension());
pressureST.resize(nQ() * j_dimension());
const Transformation &T = tc.SpaceCell().GetTransformation(faceqpoints[fid][0]);
for (int tq = 0; tq < TimeQ.size(); ++tq) {
for (int sq = 0; sq < SpaceQ.size(); ++sq) {
Point localQPoint = faceqpoints[fid][sq].CopyWithT(TimeQ.QPoint(tq)[0]);
qPoint[tq * SpaceQ.size() + sq] = TC.LocalToGlobal(localQPoint);
qNormal[tq * SpaceQ.size() + sq] = T * TC.SpaceCell().LocalFaceNormal(fid);
double w = norm(qNormal[tq * SpaceQ.size() + sq]);
qWeight[tq * SpaceQ.size() + sq] = T.Det()
* w * ww * SpaceQ.Weight(sq) * det_time * TimeQ.Weight(tq);
qNormal[tq * SpaceQ.size() + sq] /= w;
}
}
const Shape &SS = shape.GetSpaceShape();
ShapeValues<double> space_values(SpaceQ.size(), SS.size());
for (int sq = 0; sq < SpaceQ.size(); ++sq) {
for (int i = 0; i < SS.size(); ++i) {
space_values[sq][i] = SS(fid, sq, i);
}
}
const Shape &TS = shape.GetTimeShape();
ShapeValues<double> time_values(TimeQ.size(), TS.size());
for (int tq = 0; tq < TimeQ.size(); ++tq) {
for (int i = 0; i < TS.size(); ++i) {
time_values[tq][i] = TS(TimeQ.QPoint(tq), i);
}
}
for (int tq = 0; tq < nTimeQ(); ++tq) {
for (int sq = 0; sq < nSpaceQ(); ++sq) {
for (int si = 0; si < SS.size(); ++si) {
for (int ti = 0; ti < TS.size(); ++ti) {
for (int k = 0; k < dim; ++k) {
int index = (tq * nSpaceQ() + sq) * j_dimension()
+ ti * SS.size() * (1 + numL + dim) + k * SS.size() + si;
velocityST[index][k] = space_values[sq][si] * time_values[tq][ti];
}
for (int k = 0; k <= numL; ++k) {
int index = (tq * nSpaceQ() + sq) * j_dimension()
+ ti * SS.size() * (dim + 1 + numL) + (dim + k) * SS.size()
+ si;
pressureST[index] = space_values[sq][si] * time_values[tq][ti];
}
}
}
}
}
j_zero.resize(nQ() * j_dimension());
i_zero.resize(nQ() * i_dimension());
for (int q = 0; q < nQ(); ++q) {
for (int j = 0; j < j_dimension(); ++j)
j_zero[q * j_dimension() + j] =
((std::abs(Pressure(q, j)) == 0) && (norm(Velocity(q, j)) == 0));
for (int i = 0; i < i_dimension(); ++i)
i_zero[q * i_dimension() + i] =
((std::abs(Pressure(q, i)) == 0) && (norm(Velocity(q, i)) == 0));
}
}
SpaceTimeViscoAcousticDGTFaceElement::SpaceTimeViscoAcousticDGTFaceElement
(const STDiscretization &D, const VectorMatrixBase &g, const cell &tc,
int f_id, int nL, string dgdg)
:
shape(GetSpaceTimeShape(g, tc)),
faceqpoints(D.FaceQPoints()),
TC(tc), f(g.GetMesh().find_face(tc.Face(f_id))),
fid(f_id), dim(tc.dim()),
SpaceQ(D.GetCellQuad(g.GetDoFs().get_space_deg(tc))),
TimeQ(GetQuadrature("Qint1")),
numL(nL) {
qWeight.resize(nQ());
qPoint.resize(nQ());
qLocal.resize(nQ());
velocityST.resize(nQ() * j_dimension());
pressureST.resize(nQ() * j_dimension());
const Transformation &T = tc.SpaceCell().GetTransformation(SpaceQ.QPoint(0));
for (int tq = 0; tq < TimeQ.size(); ++tq) {
for (int sq = 0; sq < SpaceQ.size(); ++sq) {
qPoint[tq * SpaceQ.size() + sq] = TC.SpaceCell().
LocalToGlobal(SpaceQ.QPoint(sq)).
CopyWithT(f().t());
qWeight[tq * SpaceQ.size() + sq] = T.Det() * SpaceQ.Weight(sq);
}
}
const Shape &SS = shape.GetSpaceShape();
ShapeValues<double> space_values(SpaceQ.size(), SS.size());
for (int sq = 0; sq < SpaceQ.size(); ++sq) {
for (int i = 0; i < SS.size(); ++i) space_values[sq][i] = SS(sq, i);
}
const Shape &TS = shape.GetTimeShape();
for (int tq = 0; tq < nTimeQ(); ++tq) {
for (int sq = 0; sq < nSpaceQ(); ++sq) {
for (int si = 0; si < SS.size(); ++si) {
for (int ti = 0; ti < TS.size(); ++ti) {
double timeValue = 0.0;
if (fid == TC.Faces() - 2) timeValue = TS(Point(0.0), ti);
if (fid == TC.Faces() - 1) timeValue = TS(Point(1.0), ti);
for (int k = 0; k < dim; ++k) {
int index = (tq * nSpaceQ() + sq) * j_dimension()
+ ti * SS.size() * (1 + numL + dim) + k * SS.size() + si;
velocityST[index][k] = space_values[sq][si] * timeValue;
}
for (int k = 0; k <= numL; ++k) {
int index = (tq * nSpaceQ() + sq) * j_dimension()
+ ti * SS.size() * (dim + 1 + numL) + (dim + k) * SS.size()
+ si;
pressureST[index] = space_values[sq][si] * timeValue;
}
}
}
}
}
j_zero.resize(nQ() * j_dimension());
i_zero.resize(nQ() * i_dimension());
for (int q = 0; q < nQ(); ++q) {
for (int j = 0; j < j_dimension(); ++j)
j_zero[q * j_dimension() + j] = ((std::abs(Pressure(q, j)) == 0) &&
(norm(Velocity(q, j)) == 0));
for (int i = 0; i < i_dimension(); ++i)
i_zero[q * i_dimension() + i] = ((std::abs(Pressure(q, i)) == 0)
&& (norm(Velocity(q, i)) == 0));
}
}
SpaceTimeViscoAcousticDGTFaceElement::SpaceTimeViscoAcousticDGTFaceElement
(const STDiscretization &D, const VectorMatrixBase &g, const cell &tc,
int f_id, int nL, string dgdg, const cell &tcf)
:
shape(GetSpaceTimeShape(g, tc)),
faceqpoints(D.FaceQPoints()),
TC(tc), f(g.GetMesh().find_face(tc.Face(f_id))),
fid(f_id), dim(tc.dim()),
SpaceQ(D.GetCellQuad(g.GetDoFs().get_space_deg(tcf))),
TimeQ(GetQuadrature("Qint1")),
numL(nL) {
qWeight.resize(nQ());
qPoint.resize(nQ());
qLocal.resize(nQ());
velocityST.resize(nQ() * j_dimension());
pressureST.resize(nQ() * j_dimension());
const Transformation &T = tc.SpaceCell().GetTransformation(SpaceQ.QPoint(0));
for (int tq = 0; tq < TimeQ.size(); ++tq) {
for (int sq = 0; sq < SpaceQ.size(); ++sq) {
qPoint[tq * SpaceQ.size() + sq] = TC.SpaceCell().
LocalToGlobal(SpaceQ.QPoint(sq)).
CopyWithT(f().t());
qWeight[tq * SpaceQ.size() + sq] = T.Det() * SpaceQ.Weight(sq);
}
}
const Shape &SS = shape.GetSpaceShape();
ShapeValues<double> space_values(SpaceQ.size(), SS.size());
for (int sq = 0; sq < SpaceQ.size(); ++sq) {
Point qpLoc = TC.GlobalToLocal(qPoint[sq]);
for (int i = 0; i < SS.size(); ++i) space_values[sq][i] = SS(qpLoc, i);
}
const Shape &TS = shape.GetTimeShape();
for (int tq = 0; tq < nTimeQ(); ++tq) {
for (int sq = 0; sq < nSpaceQ(); ++sq) {
for (int si = 0; si < SS.size(); ++si) {
for (int ti = 0; ti < TS.size(); ++ti) {
double timeValue = 0.0;
if (fid == TC.Faces() - 2) timeValue = TS(Point(0.0), ti);
if (fid == TC.Faces() - 1) timeValue = TS(Point(1.0), ti);
for (int k = 0; k < dim; ++k) {