Commits (5)
...@@ -131,6 +131,21 @@ std::pair<double, double> Mesh::MeshWidth() const { ...@@ -131,6 +131,21 @@ std::pair<double, double> Mesh::MeshWidth() const {
PPM->Max(h_max, commSplit)); PPM->Max(h_max, commSplit));
} }
void Mesh::ExtractMeshWidth() {
auto temp = MeshWidth();
minMeshWidth = temp.first;
maxMeshWidth = temp.second;
}
double Mesh::MaxMeshWidth() const {
return maxMeshWidth;
}
double Mesh::MinMeshWidth() const {
return minMeshWidth;
}
Mesh::~Mesh() { while (CellCount()) RemoveCell(cells()); } Mesh::~Mesh() { while (CellCount()) RemoveCell(cells()); }
Logging &operator<<(Logging &s, const Mesh &M) { Logging &operator<<(Logging &s, const Mesh &M) {
...@@ -769,3 +784,6 @@ int Mesh::ProcSetsCountGeometryWithoutInfty() const { ...@@ -769,3 +784,6 @@ int Mesh::ProcSetsCountGeometryWithoutInfty() const {
} }
...@@ -15,311 +15,322 @@ ...@@ -15,311 +15,322 @@
class MapGeometry { class MapGeometry {
public: public:
virtual Point operator()(const Point &) const = 0; virtual Point operator()(const Point &) const = 0;
virtual const char *MeshName() const = 0; virtual const char *MeshName() const = 0;
}; };
class MapIdentity : public MapGeometry { class MapIdentity : public MapGeometry {
public: public:
Point operator()(const Point &p) const override { return p; } Point operator()(const Point &p) const override { return p; }
const char *MeshName() const override { return "Meshname not known"; } const char *MeshName() const override { return "Meshname not known"; }
}; };
class Mesh { class Mesh {
int verbose = 1; int verbose = 1;
int level;
int plevel;
int level; int commSplit;
int plevel; int dimension;
int commSplit; double scaling;
int dimension; bool identifyBnd = false;
double scaling; std::string meshName;
bool identifyBnd = false; std::string distName;
std::string meshName; std::string olapName;
std::string distName; double maxMeshWidth = 0.0;
std::string olapName; double minMeshWidth = 0.0;
protected: protected:
Cells meshCells{}; Cells meshCells{};
Edges meshEdges{}; Edges meshEdges{};
Faces meshFaces{}; Faces meshFaces{};
Vertices meshVertices{}; Vertices meshVertices{};
Cells meshOverlapCells{}; Cells meshOverlapCells{};
BoundaryFaces meshBndFaces{}; BoundaryFaces meshBndFaces{};
void InsertVertex(const Point &p); void InsertVertex(const Point &p);
bool RemoveVertex(const Point &p); bool RemoveVertex(const Point &p);
void InsertEdge(const Point &left, const Point &right); void InsertEdge(const Point &left, const Point &right);
bool RemoveEdge(const Point &edgeCenter); bool RemoveEdge(const Point &edgeCenter);
bool RefineEdge(const procset &p, Mesh &N); bool RefineEdge(const procset &p, Mesh &N);
void RemoveFace(const Point &F, const face &f, const Point &C); void RemoveFace(const Point &F, const face &f, const Point &C);
void RemoveFace(const cell &c, const Point &F); void RemoveFace(const cell &c, const Point &F);
void RefineFace(const procset &p, Mesh &N); void RefineFace(const procset &p, Mesh &N);
std::map<CELLTYPE, int> occurringCellTypes; // todo set dimension with this std::map<CELLTYPE, int> occurringCellTypes; // todo set dimension with this
public: public:
ProcSets procSets{}; ProcSets procSets{};
IdentifySets identifySets{}; IdentifySets identifySets{};
bool Parallel; // todo make private and set with information about level and plevel bool Parallel; // todo make private and set with information about level and plevel
Mesh(double scaling = 1.0, int commSplit = 0) : identifyBnd(false), scaling(scaling),
commSplit(commSplit) {
procSets.SetCommSplit(commSplit);
ExtractMeshWidth();
}
Mesh(const std::string &meshName, const std::string &distName,
const std::string &olapName, int plevel, int level,
int commSplit, double scaling)
: meshName(meshName), distName(distName), olapName(olapName),
plevel(plevel), level(level), commSplit(commSplit), scaling(scaling) {
config.get("MeshVerbose", verbose);
procSets.SetCommSplit(commSplit);
ExtractMeshWidth();
}
Mesh(const Mesh &mesh, const MapGeometry &);
Mesh(double scaling = 1.0, int commSplit = 0) : identifyBnd(false), scaling(scaling), std::string Name() const { return meshName; }
commSplit(commSplit) {
procSets.SetCommSplit(commSplit);
}
Mesh(const std::string &meshName, const std::string &distName, void Fill(const CoarseGeometry &geometry, const MapGeometry &Map);
const std::string &olapName, int plevel, int level,
int commSplit, double scaling)
: meshName(meshName), distName(distName), olapName(olapName),
plevel(plevel), level(level), commSplit(commSplit), scaling(scaling) {
config.get("MeshVerbose", verbose);
procSets.SetCommSplit(commSplit);
}
Mesh(const Mesh &mesh, const MapGeometry &); void Fill(const string &geoName, const MapGeometry &Map);
std::string Name() const { return meshName; } void Fill(const CoarseGeometry &geometry);
void Fill(const CoarseGeometry &geometry, const MapGeometry &Map); virtual void Fill(const string &geoName);
void Fill(const string &geoName, const MapGeometry &Map); void Fill(const Mesh &originMesh);
void Fill(const CoarseGeometry &geometry); virtual std::unique_ptr<Mesh> Refine();
virtual void Fill(const string &geoName); cell InsertCell(CELLTYPE tp, int flag, const vector<Point> &x);
void Fill(const Mesh &originMesh); void InsertOverlapCell(cell c);
virtual std::unique_ptr<Mesh> Refine(); void InsertOverlapCell(Cell *C);
cell InsertCell(CELLTYPE tp, int flag, const vector<Point> &x); void RemoveOverlapCell(const Point &center);
void InsertOverlapCell(cell c); void InsertFace(const Point &faceCenter, const Point &cellCenter);
void InsertOverlapCell(Cell *C); void InsertFace(const Point &faceCenter, const Face &faceValue);
void RemoveOverlapCell(const Point &center); void RemoveFace(const Point &faceCenter);
void InsertFace(const Point &faceCenter, const Point &cellCenter); virtual void RemoveCell(cell c);
void InsertFace(const Point &faceCenter, const Face &faceValue); vertex vertices() const;
void RemoveFace(const Point &faceCenter); vertex vertices_end() const;
virtual void RemoveCell(cell c); vertex find_vertex(const Point &point) const;
vertex vertices() const; edge edges() const;
vertex vertices_end() const; edge edges_end() const;
vertex find_vertex(const Point &point) const; edge find_edge(const Point &center) const;
edge edges() const; face faces() const;
edge edges_end() const; face faces_end() const;
edge find_edge(const Point &center) const; face find_face(const Point &center) const;
face faces() const; bnd_face bnd_faces() const;
face faces_end() const; bnd_face bnd_faces_end() const;
face find_face(const Point &center) const; bnd_face find_bnd_face(const Point &z) const;
bnd_face bnd_faces() const; cell cells() const;
bnd_face bnd_faces_end() const; cell cells_end() const;
bnd_face find_bnd_face(const Point &z) const; cell find_cell(const Point &center) const;
cell cells() const; cell find_cell_considering_overlap(const face &f, const Point &z) const;
cell cells_end() const; cell find_neighbour_cell(const cell &c, int f) const;
cell find_cell(const Point &center) const; cell overlap() const;
cell find_cell_considering_overlap(const face &f, const Point &z) const; cell overlap_end() const;
cell find_neighbour_cell(const cell &c, int f) const; cell find_overlap_cell(const Point &center) const;
cell overlap() const; procset procsets() const;
cell overlap_end() const; procset procsets_end() const;
cell find_overlap_cell(const Point &center) const; procset find_procset(const Point &z) const;
procset procsets() const; // Todo remove, only call via midPoint
template<class C>
procset find_procset(const C &c) const;
procset procsets_end() const; identifyset identifysets() const;
procset find_procset(const Point &z) const; identifyset identifysets_end() const;
// Todo remove, only call via midPoint identifyset find_identifyset(const Point &z) const;
template<class C>
procset find_procset(const C &c) const;
identifyset identifysets() const; int VertexCount() const;
identifyset identifysets_end() const; int VertexCountGeometry() const;
identifyset find_identifyset(const Point &z) const; int EdgeCount() const;
int VertexCount() const; int EdgeCountGeometry() const;
int VertexCountGeometry() const; int FaceCount() const;
int EdgeCount() const; int FaceCountGeometry() const;
int EdgeCountGeometry() const; int CellCount() const;
int FaceCount() const; int CellCountGeometry() const;
int FaceCountGeometry() const; int find_neighbour_face_id(const Point &f_c, const cell &cf) const;
int CellCount() const; bool onBndDG(const cell &c, int f) const;
int CellCountGeometry() const; int ProcSetsCount() const;
int find_neighbour_face_id(const Point &f_c, const cell &cf) const; int ProcSetsCountWithoutInfty() const;
bool onBndDG(const cell &c, int f) const; int ProcSetsCountGeometry() const;
int ProcSetsCount() const; int ProcSetsCountGeometryWithoutInfty() const;
int ProcSetsCountWithoutInfty() const; int OverlapCount() const;
int ProcSetsCountGeometry() const; int OverlapCountGeometry() const;
int ProcSetsCountGeometryWithoutInfty() const; int BoundaryFaceCount() const { return meshBndFaces.size(); }
int OverlapCount() const; void InsertBoundaryFace(const Point &z, int part);
int OverlapCountGeometry() const; void InsertBoundaryFace(const Point &z, const BoundaryFace &B);
int BoundaryFaceCount() const { return meshBndFaces.size(); } void RemoveBoundaryFace(const Point &x);
void InsertBoundaryFace(const Point &z, int part); auto &bnd_ref() { return meshBndFaces.ref(); }
void InsertBoundaryFace(const Point &z, const BoundaryFace &B); int BoundaryFacePart(const Point &bndFaceCenter) const;
void RemoveBoundaryFace(const Point &x); void FinishParallel();
auto &bnd_ref() { return meshBndFaces.ref(); } void Finish();
int BoundaryFacePart(const Point &bndFaceCenter) const; int dim() const { return dimension; }
void FinishParallel(); int CommSplit() const { return commSplit; }
void Finish(); void SetDim(int dim) { dimension = dim; }
int dim() const { return dimension; } bool identify() const { return identifyBnd; }
int CommSplit() const { return commSplit; } bool parallel() const { return Parallel; }
void SetDim(int dim) { dimension = dim; } bool onBnd(const Point &x) const;
bool identify() const { return identifyBnd; } double Scaling() const { return scaling; }
bool parallel() const { return Parallel; } string MeshWidthStr() const;
bool onBnd(const Point &x) const; void PrintInfo() const;
double Scaling() const { return scaling; } Point Neighbor(const Cell &c, int cFace) const {
face f = find_face(c.Face(cFace));
return f.Right();
}
string MeshWidthStr() const; std::pair<double, double> MeshWidth() const;
void PrintInfo() const; void ExtractMeshWidth();
Point Neighbor(const Cell &c, int cFace) const { double MaxMeshWidth() const;
face f = find_face(c.Face(cFace));
return f.Right();
}
std::pair<double, double> MeshWidth() const; double MinMeshWidth() const;
const std::map<CELLTYPE, int> &GetOccurringCellTypes() const { const std::map<CELLTYPE, int> &GetOccurringCellTypes() const {
return occurringCellTypes; return occurringCellTypes;
} }
// === Data Mesh Stuff === // // === Data Mesh Stuff === //
/// Dummy method for DataMeshes. See DataMesh.hpp for documentation. /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
virtual bool StoresData() const { return false; } virtual bool StoresData() const { return false; }
/// Dummy method for DataMeshes. See DataMesh.hpp for documentation. /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
virtual int GetChamber() const { return -1; } virtual int GetChamber() const { return -1; }
/// Dummy method for DataMeshes. See DataMesh.hpp for documentation. /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
virtual DataContainer &DataAt(const Point &z) const { return dummyContainer; } virtual DataContainer &DataAt(const Point &z) const { return dummyContainer; }
/// Dummy method for DataMeshes. See DataMesh.h for documentation. /// Dummy method for DataMeshes. See DataMesh.h for documentation.
virtual vector<struct DataContainer> DataIn(const cell &c) const { virtual vector<struct DataContainer> DataIn(const cell &c) const {
vector<DataContainer> data(c.Corners()); vector<DataContainer> data(c.Corners());
for (int i = 0; i < c.Corners(); i++) for (int i = 0; i < c.Corners(); i++)
data[i] = dummyContainer; data[i] = dummyContainer;
return data; return data;
} }
/// Dummy method for DataMeshes. See DataMesh.hpp for documentation. /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
virtual DataContainer &DataAverage(const cell &c) const { return dummyContainer; } virtual DataContainer &DataAverage(const cell &c) const { return dummyContainer; }
/// Dummy method for DataMeshes. See DataMesh.hpp for documentation. /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
virtual cell virtual cell
InsertCell(CELLTYPE tp, InsertCell(CELLTYPE tp,
int flag, int flag,
const vector<Point> &x, const vector<Point> &x,
vector<DataContainer> &data) { vector<DataContainer> &data) {
return InsertCell(tp, flag, x); return InsertCell(tp, flag, x);
}; };
virtual void InsertData(Point p, DataContainer &data) { virtual void InsertData(Point p, DataContainer &data) {
} }
virtual ~Mesh(); virtual ~Mesh();
friend Logging &operator<<(Logging &s, const Mesh &M); friend Logging &operator<<(Logging &s, const Mesh &M);
}; };
class BFParts { class BFParts {
int bnd[6]; int bnd[6];
int n; int n;
bool onbnd; bool onbnd;
public: public:
BFParts(const Mesh &M, const cell &c); BFParts(const Mesh &M, const cell &c);
int operator[](int i) const { return bnd[i]; } int operator[](int i) const { return bnd[i]; }
const int *operator()() const { const int *operator()() const {
if (onbnd) return bnd; if (onbnd) return bnd;
return 0; return 0;
} }
bool onBnd() const { return onbnd; } bool onBnd() const { return onbnd; }
int size() const { return n; } int size() const { return n; }
friend Logging &operator<<(Logging &s, const BFParts &BF); friend Logging &operator<<(Logging &s, const BFParts &BF);
}; };
void CoarseMeshBoundary(Mesh &M, const CoarseGeometry &CG, void CoarseMeshBoundary(Mesh &M, const CoarseGeometry &CG,
...@@ -329,17 +340,17 @@ void CoarseMeshBoundary(Mesh &M, const CoarseGeometry &CG, ...@@ -329,17 +340,17 @@ void CoarseMeshBoundary(Mesh &M, const CoarseGeometry &CG,
std::map<string, int> &surf_to_bc); std::map<string, int> &surf_to_bc);
static Point project(const Point &P) { static Point project(const Point &P) {
Point D = P - Point(10, 0, P[2]); Point D = P - Point(10, 0, P[2]);
double d = norm(D); double d = norm(D);
Point T = D / d + Point(10, 0, P[2]); Point T = D / d + Point(10, 0, P[2]);
return T; return T;
} }
static Point project2D(const Point &P) { static Point project2D(const Point &P) {
Point D = P - Point(10, 0); Point D = P - Point(10, 0);
double d = norm(D); double d = norm(D);
Point T = D / d + Point(10, 0); Point T = D / d + Point(10, 0);
return T; return T;
} }
#endif // of #ifndef _MESH_H_ #endif // of #ifndef _MESH_H_
...@@ -124,3 +124,11 @@ VtuPlot &mpp::deformed_plot(const string &name, const Vector &deformation) { ...@@ -124,3 +124,11 @@ VtuPlot &mpp::deformed_plot(const string &name, const Vector &deformation) {
std::pair<std::string, PlotStatus> mpp::save_plot(const string &filename) { std::pair<std::string, PlotStatus> mpp::save_plot(const string &filename) {
return {filename, Plotting::status}; return {filename, Plotting::status};
} }
void mpp::plot_mesh(const Mesh &mesh) {
auto name = mesh.Name();
if(name.empty()) name = "Mesh";
Plotting::InstanceM().AddPlot(name, mesh);
Plotting::InstanceM().SavePlot(name);
Plotting::InstanceM().RemovePlot(name);
}
...@@ -82,6 +82,8 @@ namespace mpp{ ...@@ -82,6 +82,8 @@ namespace mpp{
std::pair<std::string, PlotStatus> save_plot(const string& filename); std::pair<std::string, PlotStatus> save_plot(const string& filename);
constexpr PlotStatus endp = PlotStatus::CLEAR; constexpr PlotStatus endp = PlotStatus::CLEAR;
void plot_mesh(const Mesh& mesh);
static string intAsString(int i) { static string intAsString(int i) {
char buffer[256]; char buffer[256];
sprintf(buffer, "%04d", i); sprintf(buffer, "%04d", i);
...@@ -90,9 +92,4 @@ namespace mpp{ ...@@ -90,9 +92,4 @@ namespace mpp{
} }
} }
#define plot_vector(v) mpp::plot(#v) << v << mpp::save_plot(#v)
#define plot_deformed(v, u) mpp::deformed_plot(#v, u) << v << mpp::save_plot(#v)
#define plot_time_vector(v, i) mpp::plot(#v) << v << mpp::save_plot(#v+mpp::intAsString(i))
#define plot_time_deformed(v, u, i) mpp::deformed_plot(#v, u) << v << mpp::save_plot(#v+mpp::intAsString(i))
#endif //PLOTTING_HPP #endif //PLOTTING_HPP
...@@ -97,8 +97,8 @@ TimeIntegrator *TimeIntegratorCreator::Create() { ...@@ -97,8 +97,8 @@ TimeIntegrator *TimeIntegratorCreator::Create() {
return new GenericImplicitMidpointRule(linearSolvers.at(0), linearSolvers.at(1), return new GenericImplicitMidpointRule(linearSolvers.at(0), linearSolvers.at(1),
linearSolvers.at(2), linearSolvers.at(3)); linearSolvers.at(2), linearSolvers.at(3));
case CRANK_NICOLSON: case CRANK_NICOLSON:
return new CrankNicolson(linearSolvers.at(0), linearSolvers.at(1), linearSolvers.at(2), return new CrankNicolson(linearSolvers.at(0), linearSolvers.at(1),
linearSolvers.at(3)); linearSolvers.at(2), linearSolvers.at(3));
case DI_RK_ORDER2: case DI_RK_ORDER2:
return new DiagonalImplicitRungeKutta2(linearSolvers.at(0), linearSolvers.at(1), return new DiagonalImplicitRungeKutta2(linearSolvers.at(0), linearSolvers.at(1),
linearSolvers.at(2), linearSolvers.at(3)); linearSolvers.at(2), linearSolvers.at(3));
......