Mesh.cpp

#include "Mesh.hpp"
#include "Distribution.hpp"
#include "Config.hpp"
#include "Parallel.hpp"
#include "ExchangeBuffer.hpp"
#include "M_IOFiles.hpp"

#include <iostream>


using namespace std;

void Mesh::RemoveFace(const Point &F, const face &f, const Point &C) {
    if (find_cell(C) != cells_end()) return;
    if (C == Infty) procSets.RemoveIfOn(F, PPM->Proc(commSplit));
    meshFaces.Remove(F);
//    meshBndFaces.Remove(F);
}

void Mesh::RemoveFace(const cell &c, const Point &F) {
    face f = find_face(F);
    if (f == faces_end()) return;
    if (find_cell(f.Left()) == c) RemoveFace(F, f, f.Right());
    else if (find_cell(f.Right()) == c) RemoveFace(F, f, f.Left());
}

void Mesh::InsertOverlapCell(cell c) {
    meshOverlapCells
        .Insert(c.Center(), CreateCell(c.Type(), c.Subdomain(), c.AsVector()));
}

void Mesh::InsertOverlapCell(Cell *C) {
    meshOverlapCells
        .Insert(C->Center(), CreateCell(C->Type(), C->Subdomain(), C->AsVector()));
}

void Mesh::FinishParallel() {
    Parallel = PPM->Or(procsets() != procsets_end(), commSplit);
}

void Mesh::Finish() {
    Parallel = PPM->Or(procsets() != procsets_end(), commSplit);
    cell c = cells();
    if (c == cells_end()) dimension = 0;
    else dimension = c.dim();
    for (cell c2 = cells(); c2 != cells_end(); c2++) {
        if (c2.dim() != dimension) Exit(
            "No good Mesh_ " + to_string(c2.dim()) + " - " + to_string(dimension))
    }
    dimension = PPM->Max(dimension, commSplit);
    for (bnd_face b = bnd_faces(); b != bnd_faces_end(); ++b) {
        face f = find_face(b());
        if (f == faces_end()) {
            if (StoresData()) continue;
            else Exit("face not found");
        }
        cell c = find_cell(f.Left());
        if (c == cells_end()) c = find_cell(f.Right());
        if (c == cells_end()) continue;
        for (int i = 0; i < c.Faces(); ++i) {
            if (c.Face(i) != b()) continue;
            identifySets.Identify(b(), b.Part());
            for (int j = 0; j < c.FaceCorners(i); ++j)
                identifySets.Identify(c.FaceCorner(i, j), b.Part());
            for (int j = 0; j < c.FaceEdges(i); ++j)
                identifySets.Identify(c.FaceEdge(i, j), b.Part());
            i = c.Faces();
        }
    }
    identifyBnd = PPM->Or(identifysets() != identifysets_end(), commSplit);
    if (!Parallel) return;
    ExchangeBuffer E;
    for (identifyset is = identifysets(); is != identifysets_end(); ++is) {
        procset p = find_procset(is());
        if (p == procsets_end()) continue;
        for (int i = 0; i < p.size(); ++i) {
            if (p[i] == PPM->Proc(commSplit)) continue;
            E.Send(p[i]) << short(is.size()) << is();
            for (int j = 0; j < is.size(); ++j)
                E.Send(p[i]) << is[j];
        }
    }
    for (short q = 0; q < PPM->Size(commSplit); ++q)
        if (E.Send(q).size())
            E.Send(q) << short(0);
    E.Communicate();
    for (short q = 0; q < PPM->Size(commSplit); ++q) {
        if (E.Receive(q).Size() == 0) continue;
        short m;
        E.Receive(q) >> m;
        while (m) {
            Point x;
            E.Receive(q) >> x;
            for (int i = 0; i < m; ++i) {
                Point y;
                E.Receive(q) >> y;
                identifySets.Insert(x, y);
            }
            E.Receive(q) >> m;
        }
    }
    E.ClearBuffers();
}

bool Mesh::onBnd(const Point &x) const {
    return (find_bnd_face(x) != bnd_faces_end());
}

bool Mesh::onBndDG(const cell &c, int f) const {
    face ff = find_face(c.Face(f));
    return (find_bnd_face(c.Face(f)) != bnd_faces_end())
        && ((ff.Right() == Infty) || (ff.Left() == Infty));
}

std::pair<double, double> Mesh::MeshWidth() const {
    double h_min = infty;
    double h_max = 0;
    for (cell c = cells(); c != cells_end(); ++c)
        if (this->dim() != 1) {
            for (int i = 0; i < c.Edges(); ++i) {
                double h = dist(c.EdgeCorner(i, 0), c.EdgeCorner(i, 1));
                if (h < h_min) h_min = h;
                if (h > h_max) h_max = h;
            }
        } else {
            double h = dist(c.Corner(0), c.Corner(1));
            if (h < h_min) h_min = h;
            if (h > h_max) h_max = h;
        }
    return std::pair<double, double>(PPM->Min(h_min, 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()); }

Logging &operator<<(Logging &s, const Mesh &M) {
    return s << "Vertices: " << M.VertexCount() << endl << M.meshVertices.ref()
             << "Edges: " << M.EdgeCount() << endl << M.meshEdges.ref()
             << "Faces: " << M.FaceCount() << endl << M.meshFaces.ref()
             << "Cells: " << M.CellCount() << endl << M.meshCells.ref()
             << "OverlapCells: " << M.meshOverlapCells.size() << endl
             << M.meshOverlapCells.ref()
             << "BoundaryFaces: " << M.BoundaryFaceCount()
             << endl << M.meshBndFaces.ref()
             << "ProcSets: " << M.procSets.size() << endl << M.procSets.ref()
             << "IdentifySets: " << M.identifySets.size() << endl
             << M.identifySets.ref();
}

BFParts::BFParts(const Mesh &M, const cell &c) : n(c.Faces()), onbnd(false) {
    for (int i = 0; i < n; ++i) {
        bnd[i] = M.BoundaryFacePart(c.Face(i));
        if (bnd[i] != -1) onbnd = true;
    }
}

Logging &operator<<(Logging &s, const BFParts &BF) {
    for (int i = 0; i < BF.size(); ++i) s << BF[i] << " ";
    return s << endl;
}

const Point &Move(const Point &x) { return x; }

bool Mesh::RefineEdge(const procset &p, Mesh &N) {
    edge e = find_edge(p());
    if (e != edges_end()) {
        Point x = Move(e());
        N.procSets.Copy(p, 0.5 * (e.Left() + x));
        N.procSets.Copy(p, 0.5 * (e.Right() + x));
        return true;
    }
    return false;
}

void Mesh::RefineFace(const procset &p, Mesh &N) {
    face f = find_face(p());
    if (f == faces_end()) return;
    cell c = find_cell(f.Left());
    if (c == cells_end()) c = find_cell(f.Right());
    if (c == cells_end()) return;
    int i = c.facecorner(p());
    if (i == -1) return;
    Point P = Move(f->first);
    Point E0 = Move(0.5 * (c.FaceCorner(i, 0) + c.FaceCorner(i, 1)));
    Point E1 = Move(0.5 * (c.FaceCorner(i, 1) + c.FaceCorner(i, 2)));
    if (c.FaceCorners(i) == 4) {
        Point E2 = Move(0.5 * (c.FaceCorner(i, 2) + c.FaceCorner(i, 3)));
        Point E3 = Move(0.5 * (c.FaceCorner(i, 3) + c.FaceCorner(i, 0)));
        N.procSets.Copy(p, 0.5 * (E0 + P));
        N.procSets.Copy(p, 0.5 * (E1 + P));
        N.procSets.Copy(p, 0.5 * (E2 + P));
        N.procSets.Copy(p, 0.5 * (E3 + P));
        N.procSets.Copy(p, 0.25 * (c.FaceCorner(i, 0) + E0 + P + E3));
        N.procSets.Copy(p, 0.25 * (c.FaceCorner(i, 1) + E1 + P + E0));
        N.procSets.Copy(p, 0.25 * (c.FaceCorner(i, 2) + E2 + P + E1));
        N.procSets.Copy(p, 0.25 * (c.FaceCorner(i, 3) + E3 + P + E2));
    } else if (c.FaceCorners(i) == 3) {
        Point E2 = Move(0.5 * (c.FaceCorner(i, 2) + c.FaceCorner(i, 0)));
        N.procSets.Copy(p, 0.5 * (E0 + E1));
        N.procSets.Copy(p, 0.5 * (E1 + E2));
        N.procSets.Copy(p, 0.5 * (E0 + E2));
        N.procSets.Copy(p, (1 / 3.0) * (c.FaceCorner(i, 0) + E0 + E2));
        N.procSets.Copy(p, (1 / 3.0) * (c.FaceCorner(i, 1) + E0 + E1));
        N.procSets.Copy(p, (1 / 3.0) * (c.FaceCorner(i, 2) + E1 + E2));
        N.procSets.Copy(p, (1 / 3.0) * (E0 + E1 + E2));
    }
}

void CoarseMeshBoundary(Mesh &M, const CoarseGeometry &CG,
                        map<string, list<int>> &bc_nodes,
                        map<string, int> &key_to_bc,
                        map<string, list<Point>> &surf,
                        map<string, int> &surf_to_bc) {
    if (!PPM->Master(M.CommSplit())) return;
    for (face f = M.faces(); f != M.faces_end(); ++f)
        if (f.Right() == Infty)
            M.InsertBoundaryFace(f(), 0);
    vector<Point> z(CG.coordinates.begin(), CG.coordinates.end());
    for (auto &k2b_it : key_to_bc) {
        int bc_id = k2b_it.second;
        vector<int> node_id(bc_nodes[k2b_it.first].begin(),
                            bc_nodes[k2b_it.first].end());
        for (bnd_face b = M.bnd_faces(); b != M.bnd_faces_end(); ++b) {
            face f = M.find_face(b());
            cell c = M.find_cell(f.Left());
            if (c == M.cells_end()) c = M.find_cell(f.Right());
            if (c == M.cells_end()) continue;
            int cnt = 0;
            for (int i = 0; i < c.Faces(); ++i) {
                if (c.Face(i) != b()) continue;
                for (int j = 0; j < c.FaceCorners(i); ++j) {
                    Point FACECORNER = c.FaceCorner(i, j);
                    for (int k : node_id) {
                        if (FACECORNER == z[k]) ++cnt;
                    }
                }
            }
            if (cnt >= 3) {
                M.InsertBoundaryFace(f(), bc_id);
            }
        }
    }
    for (auto &it : surf) {
        int bc_id = surf_to_bc[it.first];
        for (auto &m : it.second)
            M.InsertBoundaryFace(m, bc_id);
    }
}

void Mesh::Fill(const Mesh &originMesh) {
    for (cell c = originMesh.cells(); c != originMesh.cells_end(); ++c)
        InsertCell(c.Type(), c.Subdomain(), c.AsVector());
    for (bnd_face f = originMesh.bnd_faces(); f != originMesh.bnd_faces_end(); ++f)
        if (find_face(f()) != faces_end())
            InsertBoundaryFace(f(), f.Part());
    ExchangeBuffer E;
    for (procset p = originMesh.procsets(); p != originMesh.procsets_end(); ++p)
        if ((find_vertex(p()) != vertices_end())
            || (find_edge(p()) != edges_end())
            || (find_face(p()) != faces_end())) {
            procSets.Add(p(), PPM->Proc(commSplit));
            for (int i = 0; i < p.size(); ++i)
                if (p[i] != PPM->Proc(commSplit))
                    E.Send(p[i]) << p();
        }
    E.Communicate();
    for (short q = 0; q < PPM->Size(commSplit); ++q)
        while (E.Receive(q).size() < E.ReceiveSize(q)) {
            Point z;
            E.Receive(q) >> z;
            procset p = find_procset(z);
            if (p != procsets_end())
                procSets.Add(z, q);
        }
    procSets.RemoveSingle();
    Parallel = PPM->Or(procsets() != procsets_end(), commSplit);
}

// todo check!
Mesh::Mesh(const Mesh &mesh, const MapGeometry &map) :
    meshName(mesh.meshName), distName(mesh.distName), olapName(mesh.olapName),
    plevel(mesh.plevel), level(mesh.level), commSplit(mesh.commSplit),
    scaling(mesh.scaling),
    dimension(mesh.dimension), Parallel(mesh.Parallel), identifyBnd(mesh.identifyBnd) {
    procSets.SetCommSplit(commSplit);
    for (cell c = mesh.cells(); c != mesh.cells_end(); ++c) {
        vector<Point> x(c.size());
        for (int k = 0; k < x.size(); ++k) x[k] = map(c[k]);
        cell C = InsertCell(c.Type(), c.Subdomain(), x);
        for (int k = 0; k < c.Corners(); ++k) {
            procset p = mesh.find_procset(c[k]);
            if (p != mesh.procsets_end()) procSets.Insert(C[k], p->second);
        }
        for (int i = 0; i < c.Edges(); ++i) {
            for (int k = 0; k < c.Edges(); ++k) {
                procset p = mesh.find_procset(c.Edge(k));
                if (p != mesh.procsets_end()) procSets.Insert(C.Edge(k), p->second);
            }
            for (int k = 0; k < c.Faces(); ++k) {
                procset p = mesh.find_procset(c.Face(k));
                if (p != mesh.procsets_end()) procSets.Insert(C.Face(k), p->second);
                bnd_face b = mesh.find_bnd_face(c.Face(k));
                if (b != mesh.bnd_faces_end())
                    meshBndFaces.Insert(C.Face(k), b->second);
            }
        }
    }
}

void Mesh::Fill(const CoarseGeometry &geometry) {
    Fill(geometry, MapIdentity());
}

void Mesh::Fill(const string &geoName) {
    Fill(geoName, MapIdentity());
}

void Mesh::Fill(const CoarseGeometry &geometry, const MapGeometry &Map) {
    if (!PPM->Master(commSplit)) return;

    vector<Point> z(geometry.coordinates);
    for (const auto &Cell_Id : geometry.cellIds) {
        vector<Point> x(Cell_Id.size());
        for (int k = 0; k < Cell_Id.size(); ++k) {
            x[k] = Map(z[Cell_Id.id(k)]);
        }
        InsertCell(GetCellType(Cell_Id.type(), x), Cell_Id.flag(), x);
    }

    if (geometry.faceIds.empty())
        for (face f = faces(); f != faces_end(); ++f)
            if (f.Right() == Infty)
                InsertBoundaryFace(f(), 1);
    for (const auto &i : geometry.faceIds) {
        int n = i.size();
        Assert(n == i.type());
        Point x = Origin;
        for (int k = 0; k < n; ++k) x += Map(z[i.id(k)]);
        x *= (1.0 / n);
        InsertBoundaryFace(x, i.flag());
    }
}

void Mesh::Fill(const string &geoName, const MapGeometry &Map) {
    if (!PPM->Master(commSplit)) return;

    SimpleCoarseGeometry CG(geoName);
    if (Scaling() != 1) {
        mout << "  scaling geometry by " << Scaling() << endl;
        CG.Scale(Scaling());
    }
    Fill(CG, Map);
}

void Mesh::InsertVertex(const Point &p) {
    auto check = find_vertex(p);
    if (check == vertices_end()) {
        meshVertices.Insert(p, new Vertex(0));
    } else
        check->second->inc();
}

bool Mesh::RemoveVertex(const Point &p) {
    auto v = find_vertex(p);
    if (v->second->dec() > 0)
        return false;
    meshVertices.Remove(p);
    return true;
}

vertex Mesh::vertices() const {
    return vertex(meshVertices.Begin());
}

vertex Mesh::vertices_end() const {
    return vertex(meshVertices.End());
}

vertex Mesh::find_vertex(const Point &p) const {
    return vertex(meshVertices.Find(p));
}

cell Mesh::InsertCell(CELLTYPE tp, int flag, const vector<Point> &x) {
    Cell *insertedCell = CreateCell(tp, flag, x);
    cell c = cell(meshCells.Insert(insertedCell->Center(), insertedCell));
    for (int i = 0; i < c.Corners(); ++i)
        InsertVertex(c.Corner(i));
    for (int i = 0; i < c.Edges(); ++i)
        InsertEdge(c.EdgeCorner(i, 0), c.EdgeCorner(i, 1));
    for (int i = 0; i < c.Faces(); ++i)
        InsertFace(c.Face(i), c());

    auto t = occurringCellTypes.find(tp);
    if (t == occurringCellTypes.end()) {
        occurringCellTypes[tp] = 1;
    } else {
        t->second++;
    }

    return c;
}

void Mesh::RemoveCell(cell c) {
    auto t = occurringCellTypes.find(c.Type());
    if (--t->second == 0)
        occurringCellTypes.erase(c.Type());

    for (int i = 0; i < c.Corners(); ++i)
        if (RemoveVertex(c.Corner(i)))
            procSets.RemoveIfOn(c.Corner(i), PPM->Proc(commSplit));
    for (int i = 0; i < c.Edges(); ++i)
        if (RemoveEdge(c.Edge(i)))
            procSets.RemoveIfOn(c.Edge(i), PPM->Proc(commSplit));
    for (int i = 0; i < c.Faces(); ++i)
        RemoveFace(c, c.Face(i));
    meshCells.Remove(c());
}

cell Mesh::cells() const {
    return cell(meshCells.Begin());
}

cell Mesh::cells_end() const {
    return cell(meshCells.End());
}

cell Mesh::find_cell(const Point &center) const {
    return cell(meshCells.Find(center));
}

cell Mesh::find_cell_considering_overlap(const face &f, const Point &z) const {
    cell c = find_cell(z);
    if (c == cells_end())
        c = find_overlap_cell(z);
    if (c == overlap_end()) {
        Exit("Error: no cell found")
    }
    return c;
}

cell Mesh::find_neighbour_cell(const cell &c, int f) const {
    face ff = find_face(c.Face(f));
    if (c() != ff.Right())
        return find_cell_considering_overlap(ff, ff.Right());
    else
        return find_cell_considering_overlap(ff, ff.Left());
}

int Mesh::find_neighbour_face_id(const Point &f_c, const cell &cf) const {
    for (int f1 = 0; f1 < cf.Faces(); ++f1)
        if (f_c == find_face(cf.Face(f1))()) return f1;
    Exit("Error: no face found")
}

void Mesh::InsertEdge(const Point &left, const Point &right) {
    Point middle = 0.5 * (left + right);
    auto e = meshEdges.Find(middle);
    if (e == meshEdges.End())
        meshEdges.Insert(middle, new Edge(left, right));
    else
        e->second->inc();
}

bool Mesh::RemoveEdge(const Point &edgeCenter) {
    auto e = meshEdges.Find(edgeCenter);
    if (e->second->dec() > 0)
        return false;
    meshEdges.Remove(edgeCenter);
    return true;
}

edge Mesh::edges() const {
    return edge(meshEdges.Begin());
}

edge Mesh::edges_end() const {
    return edge(meshEdges.End());
}

edge Mesh::find_edge(const Point &center) const {
    return edge(meshEdges.Find(center));
}

void Mesh::InsertFace(const Point &faceCenter, const Point &cellCenter) {
    auto f = meshFaces.Find(faceCenter);
    if (f == meshFaces.End())
        meshFaces.Insert(faceCenter, new Face(cellCenter));
    else if (f->second->Right() == Infty)
        if (f->second->Left() != cellCenter) {
            meshFaces.Replace(faceCenter, new Face(f->second->Left(), cellCenter));
        }
}

void Mesh::InsertFace(const Point &faceCenter, const Face &faceValue) {
    meshFaces.Insert(faceCenter, new Face(faceValue));
}

void Mesh::RemoveFace(const Point &faceCenter) {
    meshFaces.Remove(faceCenter);
}

face Mesh::faces() const {
    return face(meshFaces.Begin());
}

face Mesh::faces_end() const {
    return face(meshFaces.End());
}

face Mesh::find_face(const Point &center) const {
    return face(meshFaces.Find(center));
}

void Mesh::RemoveOverlapCell(const Point &center) {
    meshOverlapCells.Remove(center);
}

cell Mesh::overlap() const {
    return cell(meshOverlapCells.Begin());
}

cell Mesh::overlap_end() const {
    return cell(meshOverlapCells.End());
}

cell Mesh::find_overlap_cell(const Point &center) const {
    return cell(meshOverlapCells.Find(center));
}

void Mesh::InsertBoundaryFace(const Point &z, int part) {
    meshBndFaces.Insert(z, new BoundaryFace(part));
}

void Mesh::InsertBoundaryFace(const Point &z, const BoundaryFace &B) {
    meshBndFaces.Insert(z, new BoundaryFace(B));
}

void Mesh::RemoveBoundaryFace(const Point &x) {
    meshBndFaces.Remove(x);
}

bnd_face Mesh::bnd_faces() const {
    return bnd_face(meshBndFaces.Begin());
}

bnd_face Mesh::bnd_faces_end() const {
    return bnd_face(meshBndFaces.End());
}

bnd_face Mesh::find_bnd_face(const Point &z) const {
    return bnd_face(meshBndFaces.Find(z));
}

procset Mesh::procsets() const {
    return procset(procSets.Begin());
}

procset Mesh::procsets_end() const {
    return procset(procSets.End());
}

procset Mesh::find_procset(const Point &z) const {
    return procset(procSets.Find(z));
}

template<class C>
procset Mesh::find_procset(const C &c) const {
    return procset(procSets.Find(c()));
}

identifyset Mesh::identifysets() const {
    return identifyset(identifySets.Begin());
}

identifyset Mesh::identifysets_end() const {
    return identifyset(identifySets.End());
}

identifyset Mesh::find_identifyset(const Point &z) const {
    return identifyset(identifySets.Find(z));
}

int Mesh::BoundaryFacePart(const Point &bndFaceCenter) const {
    auto b = meshBndFaces.Find(bndFaceCenter);
    if (b == meshBndFaces.End())
        return -1;
    return b->second->Part();
}

std::unique_ptr<Mesh> Mesh::Refine() {
    dout(1) << "Refine Regular" << endl;

    auto meshPtr = std::make_unique<Mesh>(meshName, distName, olapName,
                                          plevel, level + 1, commSplit, scaling);
    auto &refinedMesh = *meshPtr;

    refinedMesh.SetDim(dimension);
    for (procset p = procsets(); p != procsets_end(); ++p) {
        refinedMesh.procSets.Copy(p);
        if (!RefineEdge(p, refinedMesh))
            RefineFace(p, refinedMesh);
    }
    vector<Point> z;
    vector<Point> x;
    vector<cell> C(8);
    for (cell c = cells(); c != cells_end(); ++c) {
        BFParts bnd(*this, c);
        const vector<Rule> &R = c.Refine(z);
        for (int i = 0; i < R.size(); ++i) {
            R[i](z, x);
            if (bnd.onBnd()) {
                Cell *t = CreateCell(R[i].type(), c.Subdomain(), x);
                for (short face = 0; face < t->Faces(); ++face) {
                    int id = R[i].face[face];
                    if (id != -1)
                        if (bnd[id] == 333)
                            for (int edge = 0; edge < t->FaceEdges(id); ++edge) {
                                int k = t->faceedge(face, edge);
                                Point P = t->Edge(k);
                                if (dimension == 3) x[8 + k] = project(P);
                                else x[4 + k] = project2D(P);
                            }
                }
                delete t;
            }
            C[i] = refinedMesh.InsertCell(R[i].type(), c.Subdomain(), x);
        }
        BFParts b(*this, c);
        if (b.onBnd())
            for (int k = 0; k < R.size(); ++k)
                for (int j = 0; j < C[k].Faces(); ++j) {
                    int i = R[k].face[j];
                    if (i != -1)
                        if (b[i] != -1)
                            refinedMesh.InsertBoundaryFace(C[k].Face(j), b[i]);
                }
    }
    refinedMesh.Finish();
//    check_mesh(refinedMesh);

    return meshPtr;
}

string Mesh::MeshWidthStr() const {
    auto h = this->MeshWidth();
    std::string strMeshWidth = "[" + std::to_string(h.first) +
        ", " + std::to_string(h.second) + "]";
    return strMeshWidth;
}

void Mesh::PrintInfo() const {
    ExchangeBuffer exBuffer(commSplit);
    exBuffer.Send(0) << MeshInfoOnProc(*this);
    exBuffer.Communicate();
    MeshInfoOnProcs infos;
    for (int p = 0; p < PPM->Size(commSplit); p++) {
        MeshInfoOnProc info;
        exBuffer.Receive(p) >> info;
        infos.push_back(info);
    }

    mout.PrintInfo(
        "Mesh", verbose,
        PrintInfoEntry("Mesh Name", meshName, 1),
        PrintInfoEntry("Distribution", distName, 1),
        PrintInfoEntry("Overlap", olapName, 1),
        PrintInfoEntry("Level", level, 1),
        PrintInfoEntry("pLevel", plevel, 1),
        PrintInfoEntry("Dimension", dim(), 1),
        PrintInfoEntry("Cells", CellCountGeometry(), 1),
        PrintInfoEntry("Overlap", OverlapCountGeometry(), 1),
        PrintInfoEntry("Vertices", VertexCountGeometry(), 1),
        PrintInfoEntry("Edges", EdgeCountGeometry(), 1),
        PrintInfoEntry("Faces", FaceCountGeometry(), 1),
        PrintInfoEntry("ProcSets", ProcSetsCountGeometryWithoutInfty(), 1),
        PrintInfoEntry("Mesh width [min, max]", MeshWidthStr(), 1),
        PrintInfoEntry("Min cells on proc", PPM->Min(CellCount(), commSplit), 2),
        PrintInfoEntry("Max cells on proc", PPM->Max(CellCount(), commSplit), 2),
        PrintInfoEntry("Communicator split", commSplit, 2),
        PrintInfoEntry("Processes", infos.ProcString(), 2),
        PrintInfoEntry("Cell distribution", infos.CellCountString(), 2),
        PrintInfoEntry("Overlap distribution", infos.OverlapCountString(), 2),
        PrintInfoEntry("Vertex distribution", infos.VertexCountString(), 2),
        PrintInfoEntry("Edges distribution", infos.EdgesCountString(), 2),
        PrintInfoEntry("Faces distribution", infos.FacesCountString(), 2),
        PrintInfoEntry("ProcSets distribution", infos.ProcSetsCountString(), 2)
    );
}

int Mesh::CellCount() const {
    return meshCells.size();
}

int Mesh::CellCountGeometry() const {
    int count = 0;
    for (cell c = this->cells(); c != this->cells_end(); ++c)
        if (this->procSets.master(c()))
            count++;
    return PPM->SumOnCommSplit(count, commSplit);
}

int Mesh::FaceCount() const {
    return meshFaces.size();
}

int Mesh::FaceCountGeometry() const {
    int count = 0;
    for (face f = this->faces(); f != this->faces_end(); ++f)
        if (this->procSets.master(f()))
            count++;
    return PPM->SumOnCommSplit(count, commSplit);
}

int Mesh::EdgeCount() const {
    return meshEdges.size();
}

int Mesh::EdgeCountGeometry() const {
    int count = 0;
    for (edge e = this->edges(); e != this->edges_end(); ++e)
        if (this->procSets.master(e()))
            count++;
    return PPM->SumOnCommSplit(count, commSplit);
}

int Mesh::VertexCount() const {
    return meshVertices.size();
}

int Mesh::VertexCountGeometry() const {
    int count = 0;
    for (vertex v = this->vertices(); v != this->vertices_end(); ++v)
        if (this->procSets.master(v()))
            count++;
    return PPM->SumOnCommSplit(count, commSplit);
}

int Mesh::OverlapCount() const {
    return meshOverlapCells.size();
}

int Mesh::OverlapCountGeometry() const {
    return PPM->SumOnCommSplit(OverlapCount(), commSplit);
}

int Mesh::ProcSetsCount() const {
    return procSets.size();
}

int Mesh::ProcSetsCountWithoutInfty() const {
    if ((procSets.size() - 1) == -1) return 0;
    else return procSets.size() - 1;
}

int Mesh::ProcSetsCountGeometry() const {
    int count = 0;
    for (procset p = this->procsets(); p != this->procsets_end(); ++p)
        if (this->procSets.master(p()))
            count++;
    return PPM->SumOnCommSplit(count, commSplit);
}

int Mesh::ProcSetsCountGeometryWithoutInfty() const {
    int count = 0;
    for (procset p = this->procsets(); p != this->procsets_end(); ++p)
        if (this->procSets.master(p()) && p() != Infty)
            count++;
    return PPM->SumOnCommSplit(count, commSplit);
}