Skip to content

GitLab

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

[Hack-For-ST-In-VtuPlot] Merge branch 'feature' into Hack-For-ST-In-VtuPlot

parents 6d9a4c74 c6836617
Hide whitespace changes
Inline Side-by-side
src/lib2_mesh/cells/Cell.hpp
View file @ d639942d
......@@ -71,6 +71,10 @@ public:
return (*this)->second->GlobalToLocal(z);
}
bool PointInCell(const Point &z) const {
return (*this)->second->PointInCell(z);
}
Point operator[](const Point &z) const {
return (*this)->second->LocalToGlobal(z);
}
......
src/lib2_mesh/cells/ICell.hpp
View file @ d639942d
......@@ -166,8 +166,13 @@ public:
virtual Point FaceLocalToGlobal(int face, const Point &) const = 0;
virtual Point LocalToGlobal(const Point &) const = 0;
virtual Point GlobalToLocal(const Point &) const = 0;
virtual bool PointInCell(const Point &) const {
Exit("Not implemented")
}
virtual Transformation GetTransformation(const Point &q = Origin) const = 0;
#ifdef BUILD_IA
......@@ -228,19 +233,10 @@ public:
return Corner(edgecorner(edge, corner));
};
/**
* Returns the index of the specified corner.
*/
virtual short edgecorner(unsigned short edge, unsigned short corner) const = 0;
/**
* Returns the amount of (d-1) dimensional faces of the cell.
*/
virtual int Faces() const = 0;
/**
* Returns the center of the (d-1) dimensional Face associated with the given integer
*/
virtual Point Face(int) const = 0;
virtual short FarFace(int) const { return -1; }
......@@ -282,18 +278,15 @@ public:
virtual double LocalFaceArea(int) const = 0;
/// Returns the area of the given face
virtual double FaceArea(int face) const = 0;
virtual const vector<Rule> &Refine(vector<Point> &) const = 0;
/// Returns the dimension of the cell
virtual int dim() const = 0;
/// Returns true, if the cell is planar
virtual bool plane() const = 0;
/// Returns the amount of children this cell generates when refined
virtual int Children() const = 0;
/// Returns the center point of the selected child
......
src/lib2_mesh/cells/Intval.cpp
View file @ d639942d
......@@ -139,4 +139,9 @@ Point Intval::GlobalToLocal(const Point &z) const {
return Point (x,0.0,0.0);
}
bool Intval::PointInCell(const Point &z) const {
return (abs(z[0] - Center()[0]) < 0.5 * abs(Corner(0)[0] - Corner(1)[0]) + GeometricTolerance);
}
src/lib2_mesh/cells/Intval.hpp
View file @ d639942d
......@@ -83,6 +83,9 @@ public:
Point Child(int i) const override;
Point GlobalToLocal(const Point &point) const override;
virtual bool PointInCell(const Point &z) const override;
};
#endif //INTVAL_H
src/lib2_mesh/cells/Quadrilateral.cpp
View file @ d639942d
......@@ -4,11 +4,11 @@
Quadrilateral::Quadrilateral(const vector<Point> &z,
short sd,
bool checkOrientation)
: Cell(z, sd), midpoint(0.25 * (corners[0] + corners[1] + corners[2] + corners[3])) {
: Cell(z, sd), midPoint(0.25 * (corners[0] + corners[1] + corners[2] + corners[3])) {
if (checkOrientation) this->checkOrientation(corners);
}
Point Quadrilateral::Center() const { return midpoint; }
Point Quadrilateral::Center() const { return midPoint; }
CELLTYPE Quadrilateral::Type() const { return QUADRILATERAL; }
......@@ -32,8 +32,8 @@ Point Quadrilateral::FaceLocalToGlobal(int face, const Point &local) const {
}
Point Quadrilateral::LocalToGlobal(const Point &z) const {
Point v = corners[0] + z[0] * (corners[1]-corners[0]);
Point w = corners[3] + z[0] * (corners[2]-corners[3]);
Point v = corners[0] + z[0] * (corners[1] - corners[0]);
Point w = corners[3] + z[0] * (corners[2] - corners[3]);
return v + z[1] * (w - v);
//TODO: Compare Performance
......@@ -79,7 +79,8 @@ Point Quadrilateral::Edge(int i) const {
short Quadrilateral::edgecorner(unsigned short i, unsigned short j) const {
if (i < Edges() && j < 2)
return (i + j) % 4;
else THROW("Not Implemented")
else
THROW("Not Implemented")
}
int Quadrilateral::Faces() const { return 4; }
......@@ -103,7 +104,8 @@ short Quadrilateral::FaceCorners(int i) const { return 2; }
short Quadrilateral::facecorner(unsigned short i, unsigned short j) const {
if (i < Faces() && j < FaceCorners(i))
return (i + j) % 4;
else THROW("Not Implemented")
else
THROW("Not Implemented")
}
short Quadrilateral::FaceEdges(int i) const { return 1; }
......@@ -149,7 +151,7 @@ Point Quadrilateral::Child(int i) const {
const Quadrilateral *Quadrilateral::ReferenceQuadrilateral() {
static const Quadrilateral
*refQuad = new Quadrilateral(Points(4, Quadrilateral::LocalCornersQuad()), -1);
*refQuad = new Quadrilateral(Points(4, Quadrilateral::LocalCornersQuad()), -1);
return refQuad;
}
......@@ -178,6 +180,24 @@ const Point *Quadrilateral::LocalCenterQuad() {
return localCenter;
}
bool Quadrilateral::PointInCell(const Point &z) const {
Point diff = Corner(0) - Corner(1);
double edgelength = std::max(abs(diff[0]), abs(diff[1]));
double l1distx = abs(z[0] - midPoint[0]);
double l1radius = edgelength * 0.5;
if (l1distx < l1radius + GeometricTolerance) {
double l1disty = abs(z[1] - midPoint[1]);
if (l1disty < l1radius + GeometricTolerance) {
if (abs(z[2] - midPoint[2]) < l1radius + GeometricTolerance) return true;
else return false;
} else {
return false;
}
}
return false;
}
Point Quadrilateral2::LocalToGlobal(const Point &z) const {
return ((1 - z[0]) * (1 - z[1]) * (1 - 2 * z[0] - 2 * z[1])) * corners[0]
+ (-z[0] * (1 - z[1]) * (1 - 2 * z[0] + 2 * z[1])) * corners[1]
......
src/lib2_mesh/cells/Quadrilateral.hpp
View file @ d639942d
......@@ -34,7 +34,7 @@ class Quadrilateral : public Cell {
}
private:
Point midpoint;
Point midPoint;
public:
Quadrilateral(const vector<Point> &z,
......@@ -75,6 +75,8 @@ public:
virtual Point LocalToGlobal(const Point &z) const override;
bool PointInCell(const Point &z) const override;
virtual Transformation GetTransformation(const Point &x) const override;
int Edges() const override;
......
src/lib2_mesh/distribution/Distribution.cpp
View file @ d639942d
......@@ -4,6 +4,7 @@
#include <algorithm>
#include <list>
#include <map>
#include <functional>
void Distribution::DistributeMesh(Mesh &mesh) {
......@@ -69,7 +70,7 @@ void Distribution::DistributeMesh(Mesh &mesh) {
if (distribute) {
for(int i = startDist; i<cells.size(); ++i){
for (int i = startDist; i < cells.size(); ++i) {
distributedCells.emplace_back(cells[i]());
}
......@@ -85,8 +86,7 @@ void Distribution::DistributeMesh(Mesh &mesh) {
<< endl;
if (PPM->Min(cellCount, commSplit) == 0)
VerboseWarning(
"At least one processor has no cells!", 1);
VerboseWarning("At least one processor has no cells!", 1);
emptyData(clearDistributed);
mout.EndBlock(true);
......@@ -97,7 +97,7 @@ void Distribution::Communicate(Mesh &mesh) {
for (int q = 0; q < PPM->Size(commSplit); ++q) {
for (auto c = markedCells[q].begin(); c != markedCells[q].end(); ++c) {
for (int i = 0; i < c->Corners(); ++i) {
mesh.procSets.Add(c->Corner(i), q);
mesh.procSets.Add(c->Corner(i), q);
}
for (int i = 0; i < c->SpaceCell().Edges(); ++i)
mesh.procSets.Add(c->SpaceCell().Edge(i), q);
......@@ -327,6 +327,21 @@ bool TLess_t(const cell &c0, const cell &c1) {
return false;
}
std::function<bool(const cell &, const cell &)> GetLess(char var) {
switch (var) {
case 'x':
return Less_x;
case 'y':
return Less_y;
case 'z':
return Less_z;
case 't':
return Less_t;
default:
THROW("Dimension " + std::to_string(var) + " not known");
}
}
void Distribution::RCB_x(int rekDepth, int begin, int end) {
if (rekDepth == 0) return;
sort(cells.begin() + begin, cells.begin() + end, Less_x);
......@@ -794,19 +809,7 @@ void Distribution::timeOpt(int tProc, int level, int begin, int end) {
void Distribution::timercb(int level, int begin, int end, char d) {
int N = end - begin;
switch (d) {
case 'x':
sort(cells.begin() + begin, cells.begin() + end, TLess_x);
break;
case 'y':
sort(cells.begin() + begin, cells.begin() + end, TLess_y);
break;
case 't':
sort(cells.begin() + begin, cells.begin() + end, TLess_t);
break;
default:
break;
}
sort(cells.begin() + begin, cells.begin() + end, GetLess(d));
if (level == 0) {
return;
} else {
......@@ -865,23 +868,8 @@ void Distribution::timercb_weighted(int tlevel, int xlevel, int ylevel,
int default_weight) {
int N = end - begin;
int level = tlevel + xlevel + ylevel;
switch (d) {
case 'x':
sort(cells.begin() + begin, cells.begin() + end, TLess_x);
break;
case 'y':
sort(cells.begin() + begin, cells.begin() + end, TLess_y);
break;
case 't':
sort(cells.begin() + begin, cells.begin() + end, TLess_t);
break;
default:
break;
}
sort(cells.begin() + begin, cells.begin() + end, GetLess(d));
if (PPM->proc() == 0) {
int W_sum = 0;
for (int i = begin; i < end; ++i) {
if (W.find(cells[i]()) != W.end())
......@@ -919,7 +907,7 @@ void Distribution::timercb_weighted(int tlevel, int xlevel, int ylevel,
while (cells[mid]().t() == cells[mid - 1]().t() || a == b) {
--mid;
if (mid < begin) break;
if (mid < begin || mid == 0) break;
a = cells[mid]().CopyWithT(0.0);
b = cells[mid + 1]().CopyWithT(0.0);
}
......@@ -1046,8 +1034,13 @@ void Distribution::DistributeSTMesh(Mesh &mesh) {
} else if (_distName == "deformed_optimized") {
int S = mesh.steps() - 1;
int Slice_Maximum = int(log(double(S)) / log(2.0) + 1e-10);
int X_Maximum = int(log(double(sqrt(N / (S)))) / log(2.0) + 1e-10);
int X_Maximum = int(log(double(sqrt(N / S))) / log(2.0) + 1e-10);
int Y_Maximum = X_Maximum;
if (Slice_Maximum == 0 && X_Maximum == 0 && Y_Maximum == 0) {
THROW("Not enough cells to distribute. Level: " + mesh.Level().str())
}
int tL = 0;
int yL = 0;
int xL = 0;
......
src/lib2_mesh/distribution/Overlap.cpp
View file @ d639942d
......@@ -707,7 +707,7 @@ void Overlap::CellsOverlap_dG2(Mesh &mesh) const {
}
}
mesh.InsertOverlapCell(c->second);
//mesh.InsertOverlapCell(c->second);
mesh.procSets.Add(c(), PPM->proc());
if (!send_list.empty()) {
......@@ -752,4 +752,5 @@ void Overlap::CellsOverlap_dG2(Mesh &mesh) const {
delete C;
}
}
mesh.procSets.RemoveSingle();
}
src/lib2_mesh/mesh/CoarseGeometry.cpp
View file @ d639942d
......@@ -77,7 +77,7 @@ std::vector<double> CoarseGeometry::refineTimeSteps(const std::vector<double> &t
double dt = init_timesteps[S - 1][k + 2] - init_timesteps[S - 1][k + 1];
for (unsigned int l = 1; l < SF_inv; ++l)
init_timesteps[S][k * SF_inv + l + 1] =
init_timesteps[S - 1][k + 1] + dt * 1 / double(SF_inv) * l;
init_timesteps[S - 1][k + 1] + dt * 1 / double(SF_inv) * l;
}
}
return init_timesteps[refineTimeStepCount];
......@@ -136,7 +136,7 @@ void CoarseGeometry::CommunicateGeometry() {
ExchangeBuffer exBuffer;
if (PPM->Master()) {
for (int q = 1; q < PPM->Size(); ++q) {
exBuffer.Send(q) << timeSteps.size();
exBuffer.Send(q) << int(timeSteps.size());
for (double timeStep : timeSteps) {
exBuffer.Send(q) << double(timeStep);
}
......@@ -147,8 +147,9 @@ void CoarseGeometry::CommunicateGeometry() {
int m;
exBuffer.Receive(0) >> m;
timeSteps.resize(m);
for (double & timeStep : timeSteps) {
exBuffer.Receive(0) >> timeStep;
for (int i = 0; i < m; i++) {
exBuffer.Receive(0) >> timeSteps[i];
}
}
exBuffer.ClearBuffers();
......@@ -157,7 +158,7 @@ void CoarseGeometry::CommunicateGeometry() {
CoarseGeometry::CoarseGeometry(const MeshPoints &points, const MeshCells &cells,
const MeshFaces &faces, const TimeSteps &timeSteps) :
timeSteps(timeSteps) {
timeSteps(timeSteps) {
for (const auto &point : points) {
double z[]{0.0, 0.0, 0.0, 0.0};
for (int j = 0; j < min(3, int(point.size())); ++j)
......@@ -175,8 +176,8 @@ void CoarseGeometry::Scale(double scalingFactor) {
Logging &operator<<(Logging &s, const CoarseGeometry &M) {
return s << "POINTS: " << M.coordinates.size() << endl << M.coordinates
<< "CELLS: " << M.cellIds.size() << endl << M.cellIds
<< "FACES: " << M.faceIds.size() << endl << M.faceIds
<< "CELLS: " << M.cellIds.size() << endl << M.cellIds
<< "FACES: " << M.faceIds.size() << endl << M.faceIds
#ifdef USE_DATAMESH
<< "VDATA: " << M.vdataList.size() << endl << M.vdataList
<< "CDATA: " << M.cdataList.size() << endl << M.cdataList
......
src/lib2_mesh/mesh/Mesh.hpp
View file @ d639942d
......@@ -24,6 +24,13 @@ struct LevelPair {
LevelPair CoarserInTime() const { return LevelPair{space, time - 1}; }
LevelPair Next() const { return LevelPair{space + 1, time + 1}; }
std::string str() const {
std::stringstream ss;
ss << "[" << space << ", " << time << "]";
return ss.str();
}
};
inline bool operator==(const LevelPair &l1, const LevelPair &l2) {
......
src/lib4_fem/elements/DGElement.cpp
View file @ d639942d
#include "DGElement.hpp"
[[deprecated("Simply Use GlobalToLocal of Cell")]]
Point AcousticDGElement::GlobalToLocal(const Point &z) const {
switch(C.Type()) {
case INTERVAL: {
Point p(z);
p[0] -= C.Corner(0)[0];
double length = abs(C.Corner(0)[0] - C.Corner(1)[0]);
p /= length;
return p;
}
case TRIANGLE: {
Tensor T;
T[0][0] = C.Corner(1)[0] - C.Corner(0)[0];
T[0][1] = C.Corner(2)[0] - C.Corner(0)[0];
T[1][0] = C.Corner(1)[1] - C.Corner(0)[1];
T[1][1] = C.Corner(2)[1] - C.Corner(0)[1];
T[0][2] = 0;
T[1][2] = 0;
T[2][0] = 0;
T[2][1] = 0;
T[2][2] = 1;
Tensor IT = Invert(T);
return IT * (z - C.Corner(0));
}
case TETRAHEDRON: Exit("not implemented");
case QUADRILATERAL: {
if (inverseTensor == false) {
IT = new Tensor;
(*IT)[0][0] = C.Corner(1)[0] - C.Corner(0)[0];
(*IT)[0][1] = C.Corner(3)[0] - C.Corner(0)[0];
(*IT)[1][0] = C.Corner(1)[1] - C.Corner(0)[1];
(*IT)[1][1] = C.Corner(3)[1] - C.Corner(0)[1];
(*IT)[0][2] = 0;
(*IT)[1][2] = 0;
(*IT)[2][0] = 0;
(*IT)[2][1] = 0;
(*IT)[2][2] = 1;
IT->Invert();
inverseTensor = true;
}
return (*IT) * (z - C.Corner(0));
}
case HEXAHEDRON: Exit("not implemented");
default: Exit("not implemented");
}
Exit("Not implemented!");
return Origin;
return C.GlobalToLocal(z);
}
AcousticDGElement::AcousticDGElement(const VectorMatrixBase &base,
......
src/lib4_fem/plot/VtuPlot.cpp
View file @ d639942d
......@@ -86,6 +86,8 @@ void VtuPlot::pointDataStream(std::ostream &out, int k) {
out << aArray;
}
out << ">" << endl;
for (const auto &pData: pointData[plotMesh().t(k)]) {
if (pData.first == "ID") continue;
......@@ -101,6 +103,22 @@ void VtuPlot::pointDataStream(std::ostream &out, int k) {
}
out << endl << indentString(4) << "</DataArray>" << endl;
}
if (isDeformed) {
out << indentString(4)
<< R"(<DataArray type="Float32" Name="U)"
<< R"(" NumberOfComponents="3)"
<< R"(" format="ascii">)";
for (int i = 0; i < points[plotMesh().t(k)].size(); i++) {
for (int j = 0; j < 3; j++) {
Point displacement = displacements[plotMesh().t(k)].at(i);
out << lineBreakString(i + j, lineBreakZ, 5) << to_string(displacement[j]);
}
}
out << endl << indentString(4) << "</DataArray>" << endl;
}
out << indentString(3) << "</PointData>" << endl;
}
......@@ -124,8 +142,8 @@ void VtuPlot::cellDataStream(std::ostream &out, int k) {
if (!activeArrays.second.empty()) {
string aArray = " ";
aArray +=
(cellDataSize[plotMesh().slice(k)].at(activeArrays.second) > 1 ? R"(Vectors=")"
: R"(Scalars=")");
(cellDataSize[plotMesh().slice(k)].at(activeArrays.second) > 1 ? R"(Vectors=")"
: R"(Scalars=")");
aArray += activeArrays.second;
aArray += R"(")";
out << aArray;
......@@ -252,12 +270,12 @@ void VtuPlot::initCells(const Mesh &init) {
try {
if (i == 0) {
cellConnectivity[center.t()] +=
lineBreakString(cells[center.t()].size() - 1, lineBreakC, 5);
lineBreakString(cells[center.t()].size() - 1, lineBreakC, 5);
cellConnectivity[center.t()] +=
to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
} else {
cellConnectivity[center.t()] +=
" " + to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
" " + to_string((int) pointData[corner.t()]["ID"].at(corner).GetData(0));
}
}
catch (...) {
......@@ -266,11 +284,11 @@ void VtuPlot::initCells(const Mesh &init) {
}
}
cellTypes[center.t()] +=
lineBreakString(cells[center.t()].size() - 1, lineBreakZ, 5);
lineBreakString(cells[center.t()].size() - 1, lineBreakZ, 5);
cellTypes[center.t()] += to_string(VtkCellType(cType));
totalOffset[center.t()] += m;
cellOffset[center.t()] +=
lineBreakString(cells[center.t()].size() - 1, lineBreakZ, 5);
lineBreakString(cells[center.t()].size() - 1, lineBreakZ, 5);
if (init.IsSTMesh())
cellOffset[center.t()] += to_string(totalOffset[center.t()] / 2);
else
......@@ -338,7 +356,7 @@ void VtuPlot::addData(const string &name, const Vector &data, std::vector<int> i
}
}
if (typeid(discToPlot) == typeid(LagrangeDiscretization) ||
typeid(discToPlot) == typeid(MultiPartDiscretization)) {
typeid(discToPlot) == typeid(MultiPartDiscretization)) {
int d = ((const LagrangeDiscretization &) discToPlot).Degree();
if (d > 0) addPointData(name, data, indices);
else addCellData(name, data, indices);
......@@ -347,25 +365,25 @@ void VtuPlot::addData(const string &name, const Vector &data, std::vector<int> i
} else if (typeid(discToPlot) == typeid(DoFDiscretization<VertexDoF>)) {
addPointData(name, data, indices);
} else if (typeid(discToPlot) == typeid(DGDiscretization) ||
typeid(discToPlot) == typeid(DoFDiscretization<DGDoF>)) {
typeid(discToPlot) == typeid(DoFDiscretization<DGDoF>)) {
addCellData(name, data, indices);
} else if (typeid(discToPlot) == typeid(ArgyrisDiscretization)) {
addArgyrisData(name, data, indices);
} else if (typeid(discToPlot) == typeid(DivergenceFreeDiscretization)) {
addDivergenceFreeData(name, data);
#if USE_SPACETIME
} else if (typeid(discToPlot) == typeid(STDiscretization_DGDG)) {
addSpaceTimeData(name, data, indices);
} else if (typeid(discToPlot) == typeid(STDiscretization_DGDG_GLGL) ||
typeid(discToPlot) == typeid(STDiscretization_PGDG)) {
Warning("Plotting not implemented for"
" STDiscretization_DGDG_GLGL or"
" STDiscretization_PGDG")
} else if (typeid(discToPlot) == typeid(STDiscretization_DGDG)) {
addSpaceTimeData(name, data, indices);
} else if (typeid(discToPlot) == typeid(STDiscretization_DGDG_GLGL) ||
typeid(discToPlot) == typeid(STDiscretization_PGDG)) {
Warning("Plotting not implemented for"
" STDiscretization_DGDG_GLGL or"
" STDiscretization_PGDG")
#endif
} else {
Warning("Plot type for Discretization "
+ std::string(typeid(discToPlot).name())
+ " not defined.")
+ std::string(typeid(discToPlot).name())
+ " not defined.")
}
}
......@@ -517,6 +535,7 @@ void VtuPlot::DeformGrid(const Vector &displacement) {
} else {
deformLagrangeGrid(displacement);
}
isDeformed = true;
}
void VtuPlot::Erase() {
......@@ -538,7 +557,7 @@ void VtuPlot::ClearData() {
dataToDrop.clear();
for (const auto &cData: cellData[0.0]) {
if (std::string{"IDSubdomainProcLoad"}.find(cData.first)
!= std::string::npos)
!= std::string::npos)
continue;
dataToDrop.push_back(cData.first);
}
......@@ -616,7 +635,7 @@ void VtuPlot::deformDGGrid(const Vector &displacement) {
DGVectorFieldElement E(displacement, c);
for (int i = 0; i < c.Corners(); ++i) {
exBuffer.Send(0) << c.Corner(i);
VectorField u = E.VectorValue(c.Corner(i), displacement);
VectorField u = E.VectorValue(c.LocalCorner(i), displacement);
for (int j = 0; j < dim; ++j)
exBuffer.Send(0) << u[j];
}
......@@ -645,9 +664,9 @@ void VtuPlot::deformDGGrid(const Vector &displacement) {
}
// for (auto &p: points) {
for (int i = 0; i < points.size(); i++) {
for (int i = 0; i < points[0.0].size(); i++) {
displacements[0.0][i] =
deformationData[points[0.0][i]] / deformationCount[points[0.0][i]];
deformationData[points[0.0][i]] / deformationCount[points[0.0][i]];
}
// }
......@@ -680,7 +699,7 @@ void VtuPlot::deformLagrangeGrid(const Vector &displacement) {
deformationData->insert({z, d});
}