src/lib2_mesh/mesh/Mesh.cpp

@@ -131,6 +131,21 @@ std::pair<double, double> Mesh::MeshWidth() const {
                                      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()); }
 
 Logging &operator<<(Logging &s, const Mesh &M) {
@@ -769,3 +784,6 @@ int Mesh::ProcSetsCountGeometryWithoutInfty() const {
 }
 
 
+
+
+

src/lib2_mesh/mesh/Mesh.hpp

@@ -15,311 +15,322 @@
 
 class MapGeometry {
 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 {
 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 {
-    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:
-    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:
-    ProcSets procSets{};
-    IdentifySets identifySets{};
-    bool Parallel; // todo make private and set with information about level and plevel
+  ProcSets procSets{};
+  IdentifySets identifySets{};
+  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),
-                                                    commSplit(commSplit) {
-        procSets.SetCommSplit(commSplit);
-    }
+  std::string Name() const { return meshName; }
 
-    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);
-    }
+  void Fill(const CoarseGeometry &geometry, const MapGeometry &Map);
 
-    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
-    template<class C>
-    procset find_procset(const C &c) const;
+  identifyset find_identifyset(const Point &z) 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 {
-        face f = find_face(c.Face(cFace));
-        return f.Right();
-    }
+  double MaxMeshWidth() const;
 
-    std::pair<double, double> MeshWidth() const;
+  double MinMeshWidth() const;
 
-    const std::map<CELLTYPE, int> &GetOccurringCellTypes() const {
-        return occurringCellTypes;
-    }
+  const std::map<CELLTYPE, int> &GetOccurringCellTypes() const {
+    return occurringCellTypes;
+  }
 
-    // === Data Mesh Stuff === //
-    /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
-    virtual bool StoresData() const { return false; }
+  // === Data Mesh Stuff === //
+  /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
+  virtual bool StoresData() const { return false; }
 
-    /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
-    virtual int GetChamber() const { return -1; }
+  /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
+  virtual int GetChamber() const { return -1; }
 
-    /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
-    virtual DataContainer &DataAt(const Point &z) const { return dummyContainer; }
+  /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
+  virtual DataContainer &DataAt(const Point &z) const { return dummyContainer; }
 
-    /// Dummy method for DataMeshes. See DataMesh.h for documentation.
-    virtual vector<struct DataContainer> DataIn(const cell &c) const {
-        vector<DataContainer> data(c.Corners());
-        for (int i = 0; i < c.Corners(); i++)
-            data[i] = dummyContainer;
-        return data;
-    }
+  /// Dummy method for DataMeshes. See DataMesh.h for documentation.
+  virtual vector<struct DataContainer> DataIn(const cell &c) const {
+    vector<DataContainer> data(c.Corners());
+    for (int i = 0; i < c.Corners(); i++)
+      data[i] = dummyContainer;
+    return data;
+  }
 
-    /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
-    virtual DataContainer &DataAverage(const cell &c) const { return dummyContainer; }
+  /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
+  virtual DataContainer &DataAverage(const cell &c) const { return dummyContainer; }
 
-    /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
-    virtual cell
-    InsertCell(CELLTYPE tp,
-               int flag,
-               const vector<Point> &x,
-               vector<DataContainer> &data) {
-        return InsertCell(tp, flag, x);
-    };
+  /// Dummy method for DataMeshes. See DataMesh.hpp for documentation.
+  virtual cell
+  InsertCell(CELLTYPE tp,
+             int flag,
+             const vector<Point> &x,
+             vector<DataContainer> &data) {
+    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 {
-    int bnd[6];
-    int n;
-    bool onbnd;
+  int bnd[6];
+  int n;
+  bool onbnd;
 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 {
-        if (onbnd) return bnd;
-        return 0;
-    }
+  const int *operator()() const {
+    if (onbnd) return bnd;
+    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,
@@ -329,17 +340,17 @@ void CoarseMeshBoundary(Mesh &M, const CoarseGeometry &CG,
                         std::map<string, int> &surf_to_bc);
 
 static Point project(const Point &P) {
-    Point D = P - Point(10, 0, P[2]);
-    double d = norm(D);
-    Point T = D / d + Point(10, 0, P[2]);
-    return T;
+  Point D = P - Point(10, 0, P[2]);
+  double d = norm(D);
+  Point T = D / d + Point(10, 0, P[2]);
+  return T;
 }
 
 static Point project2D(const Point &P) {
-    Point D = P - Point(10, 0);
-    double d = norm(D);
-    Point T = D / d + Point(10, 0);
-    return T;
+  Point D = P - Point(10, 0);
+  double d = norm(D);
+  Point T = D / d + Point(10, 0);
+  return T;
 }
 
 #endif // of #ifndef _MESH_H_

src/lib4_vector/plot/Plotting.cpp

@@ -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) {
   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);
+}

src/lib4_vector/plot/Plotting.hpp

@@ -82,6 +82,8 @@ namespace mpp{
   std::pair<std::string, PlotStatus>  save_plot(const string& filename);
   constexpr PlotStatus endp = PlotStatus::CLEAR;
 
+  void plot_mesh(const Mesh& mesh);
+
   static string intAsString(int i) {
     char buffer[256];
     sprintf(buffer, "%04d", i);
@@ -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

src/lib6_solvers/timeintegrator/TimeIntegrator.cpp

@@ -97,8 +97,8 @@ TimeIntegrator *TimeIntegratorCreator::Create() {
         return new GenericImplicitMidpointRule(linearSolvers.at(0), linearSolvers.at(1),
                                                linearSolvers.at(2), linearSolvers.at(3));
       case CRANK_NICOLSON:
-        return new CrankNicolson(linearSolvers.at(0), linearSolvers.at(1), linearSolvers.at(2),
-                                 linearSolvers.at(3));
+        return new CrankNicolson(linearSolvers.at(0), linearSolvers.at(1),
+                                 linearSolvers.at(2), linearSolvers.at(3));
       case DI_RK_ORDER2:
         return new DiagonalImplicitRungeKutta2(linearSolvers.at(0), linearSolvers.at(1),
                                                linearSolvers.at(2), linearSolvers.at(3));