Mesh.cpp 24.9 KB
Newer Older
1
#include "Mesh.hpp"
2
#include "Distribution.hpp"
3
4
5
6
#include "Config.hpp"
#include "Parallel.hpp"
#include "ExchangeBuffer.hpp"
#include "M_IOFiles.hpp"
7
8
9
10
11
12
13

#include <iostream>


using namespace std;

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

void Mesh::RemoveFace(const cell &c, const Point &F) {
21
22
23
24
    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());
25
26
27
}

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

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

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

void Mesh::Finish() {
42
    Parallel = PPM->Or(procsets() != procsets_end(), commSplit);
43
    cell c = cells();
44
45
    if (c == cells_end()) dimension = 0;
    else dimension = c.dim();
46
    for (cell c2 = cells(); c2 != cells_end(); c2++) {
47
48
        if (c2.dim() != dimension) Exit(
            "No good Mesh_ " + to_string(c2.dim()) + " - " + to_string(dimension))
49
    }
50
    dimension = PPM->Max(dimension, commSplit);
51
52
53
54
55
56
57
58
59
60
61
    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;
62
            identifySets.Identify(b(), b.Part());
63
            for (int j = 0; j < c.FaceCorners(i); ++j)
64
                identifySets.Identify(c.FaceCorner(i, j), b.Part());
65
            for (int j = 0; j < c.FaceEdges(i); ++j)
66
                identifySets.Identify(c.FaceEdge(i, j), b.Part());
67
68
            i = c.Faces();
        }
69
    }
70
    identifyBnd = PPM->Or(identifysets() != identifysets_end(), commSplit);
71
72
73
74
75
76
    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) {
77
            if (p[i] == PPM->Proc(commSplit)) continue;
78
79
80
81
            E.Send(p[i]) << short(is.size()) << is();
            for (int j = 0; j < is.size(); ++j)
                E.Send(p[i]) << is[j];
        }
82
    }
83
    for (short q = 0; q < PPM->Size(commSplit); ++q)
84
85
86
        if (E.Send(q).size())
            E.Send(q) << short(0);
    E.Communicate();
87
    for (short q = 0; q < PPM->Size(commSplit); ++q) {
88
89
90
91
92
93
94
95
96
        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;
97
                identifySets.Insert(x, y);
98
99
100
            }
            E.Receive(q) >> m;
        }
101
    }
102
    E.ClearBuffers();
103
104
105
}

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

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

std::pair<double, double> Mesh::MeshWidth() const {
116
117
118
    double h_min = infty;
    double h_max = 0;
    for (cell c = cells(); c != cells_end(); ++c)
119
        if (this->dim() != 1) {
120
121
122
123
124
125
126
127
128
129
            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;
        }
130
131
    return std::pair<double, double>(PPM->Min(h_min, commSplit),
                                     PPM->Max(h_max, commSplit));
132
133
}

134
135
136
137
138
139
140
141
142
143
144
145
146
147
148

void Mesh::ExtractMeshWidth() {
  auto temp = MeshWidth();
  minMeshWidth = temp.first;
  maxMeshWidth = temp.second;
}

double Mesh::MaxMeshWidth() const {
  return maxMeshWidth;
}

double Mesh::MinMeshWidth() const {
  return minMeshWidth;
}

149
150
151
Mesh::~Mesh() { while (CellCount()) RemoveCell(cells()); }

Logging &operator<<(Logging &s, const Mesh &M) {
152
153
154
155
    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()
156
157
             << "OverlapCells: " << M.meshOverlapCells.size() << endl
             << M.meshOverlapCells.ref()
158
159
             << "BoundaryFaces: " << M.BoundaryFaceCount()
             << endl << M.meshBndFaces.ref()
160
             << "ProcSets: " << M.procSets.size() << endl << M.procSets.ref()
161
162
             << "IdentifySets: " << M.identifySets.size() << endl
             << M.identifySets.ref();
163
164
165
}

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

172
Logging &operator<<(Logging &s, const BFParts &BF) {
173
174
    for (int i = 0; i < BF.size(); ++i) s << BF[i] << " ";
    return s << endl;
175
176
177
178
179
}

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

bool Mesh::RefineEdge(const procset &p, Mesh &N) {
180
181
182
    edge e = find_edge(p());
    if (e != edges_end()) {
        Point x = Move(e());
183
184
        N.procSets.Copy(p, 0.5 * (e.Left() + x));
        N.procSets.Copy(p, 0.5 * (e.Right() + x));
185
186
187
        return true;
    }
    return false;
188
189
190
}

void Mesh::RefineFace(const procset &p, Mesh &N) {
191
192
193
194
195
196
197
198
199
200
201
202
203
    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)));
204
205
206
207
208
209
210
211
        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));
212
213
    } else if (c.FaceCorners(i) == 3) {
        Point E2 = Move(0.5 * (c.FaceCorner(i, 2) + c.FaceCorner(i, 0)));
214
215
216
217
218
219
220
        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));
221
    }
222
223
224
225
226
227
228
}

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) {
229
    if (!PPM->Master(M.CommSplit())) return;
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
    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);
            }
        }
257
    }
258
259
260
261
262
    for (auto &it : surf) {
        int bc_id = surf_to_bc[it.first];
        for (auto &m : it.second)
            M.InsertBoundaryFace(m, bc_id);
    }
263
264
265
}

void Mesh::Fill(const Mesh &originMesh) {
266
267
268
269
270
271
272
273
274
275
    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())) {
276
            procSets.Add(p(), PPM->Proc(commSplit));
277
            for (int i = 0; i < p.size(); ++i)
278
                if (p[i] != PPM->Proc(commSplit))
279
280
281
                    E.Send(p[i]) << p();
        }
    E.Communicate();
282
    for (short q = 0; q < PPM->Size(commSplit); ++q)
283
284
285
286
287
        while (E.Receive(q).size() < E.ReceiveSize(q)) {
            Point z;
            E.Receive(q) >> z;
            procset p = find_procset(z);
            if (p != procsets_end())
288
                procSets.Add(z, q);
289
        }
290
    procSets.RemoveSingle();
291
    Parallel = PPM->Or(procsets() != procsets_end(), commSplit);
292
}
293

294
// todo check!
295
Mesh::Mesh(const Mesh &mesh, const MapGeometry &map) :
296
297
    meshName(mesh.meshName), distName(mesh.distName), olapName(mesh.olapName),
    plevel(mesh.plevel), level(mesh.level), commSplit(mesh.commSplit),
298
299
    scaling(mesh.scaling),
    dimension(mesh.dimension), Parallel(mesh.Parallel), identifyBnd(mesh.identifyBnd) {
300
    procSets.SetCommSplit(commSplit);
301
    for (cell c = mesh.cells(); c != mesh.cells_end(); ++c) {
302
        vector<Point> x(c.size());
303
        for (int k = 0; k < x.size(); ++k) x[k] = map(c[k]);
304
305
        cell C = InsertCell(c.Type(), c.Subdomain(), x);
        for (int k = 0; k < c.Corners(); ++k) {
306
            procset p = mesh.find_procset(c[k]);
307
            if (p != mesh.procsets_end()) procSets.Insert(C[k], p->second);
308
        }
309
        for (int i = 0; i < c.Edges(); ++i) {
310
            for (int k = 0; k < c.Edges(); ++k) {
311
                procset p = mesh.find_procset(c.Edge(k));
312
                if (p != mesh.procsets_end()) procSets.Insert(C.Edge(k), p->second);
313
314
            }
            for (int k = 0; k < c.Faces(); ++k) {
315
                procset p = mesh.find_procset(c.Face(k));
316
                if (p != mesh.procsets_end()) procSets.Insert(C.Face(k), p->second);
317
318
                bnd_face b = mesh.find_bnd_face(c.Face(k));
                if (b != mesh.bnd_faces_end())
319
320
321
                    meshBndFaces.Insert(C.Face(k), b->second);
            }
        }
322
    }
323
324
}

325
326
327
328
329
330
331
332
void Mesh::Fill(const CoarseGeometry &geometry) {
    Fill(geometry, MapIdentity());
}

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

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

336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
    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());
352
        Point x = Origin;
353
354
355
        for (int k = 0; k < n; ++k) x += Map(z[i.id(k)]);
        x *= (1.0 / n);
        InsertBoundaryFace(x, i.flag());
356
357
358
359
    }
}

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

362
363
364
365
366
367
    SimpleCoarseGeometry CG(geoName);
    if (Scaling() != 1) {
        mout << "  scaling geometry by " << Scaling() << endl;
        CG.Scale(Scaling());
    }
    Fill(CG, Map);
368
369
370
}

void Mesh::InsertVertex(const Point &p) {
371
372
373
374
375
    auto check = find_vertex(p);
    if (check == vertices_end()) {
        meshVertices.Insert(p, new Vertex(0));
    } else
        check->second->inc();
376
377
378
}

bool Mesh::RemoveVertex(const Point &p) {
379
380
381
382
383
    auto v = find_vertex(p);
    if (v->second->dec() > 0)
        return false;
    meshVertices.Remove(p);
    return true;
384
385
386
}

vertex Mesh::vertices() const {
387
    return vertex(meshVertices.Begin());
388
389
390
}

vertex Mesh::vertices_end() const {
391
    return vertex(meshVertices.End());
392
393
394
}

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

cell Mesh::InsertCell(CELLTYPE tp, int flag, const vector<Point> &x) {
399
400
401
402
403
404
405
406
407
    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());

408
409
410
    auto t = occurringCellTypes.find(tp);
    if (t == occurringCellTypes.end()) {
        occurringCellTypes[tp] = 1;
411
412
413
    } else {
        t->second++;
    }
414

415
    return c;
416
417
418
}

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

423
424
    for (int i = 0; i < c.Corners(); ++i)
        if (RemoveVertex(c.Corner(i)))
425
            procSets.RemoveIfOn(c.Corner(i), PPM->Proc(commSplit));
426
427
    for (int i = 0; i < c.Edges(); ++i)
        if (RemoveEdge(c.Edge(i)))
428
            procSets.RemoveIfOn(c.Edge(i), PPM->Proc(commSplit));
429
430
431
    for (int i = 0; i < c.Faces(); ++i)
        RemoveFace(c, c.Face(i));
    meshCells.Remove(c());
432
433
434
}

cell Mesh::cells() const {
435
    return cell(meshCells.Begin());
436
437
438
}

cell Mesh::cells_end() const {
439
    return cell(meshCells.End());
440
441
442
}

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

cell Mesh::find_cell_considering_overlap(const face &f, const Point &z) const {
447
448
449
450
451
452
453
    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;
454
455
456
}

cell Mesh::find_neighbour_cell(const cell &c, int f) const {
457
458
459
460
461
    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());
462
463
464
}

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

void Mesh::InsertEdge(const Point &left, const Point &right) {
471
472
473
474
475
476
    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();
477
478
479
}

bool Mesh::RemoveEdge(const Point &edgeCenter) {
480
481
482
483
484
    auto e = meshEdges.Find(edgeCenter);
    if (e->second->dec() > 0)
        return false;
    meshEdges.Remove(edgeCenter);
    return true;
485
486
487
}

edge Mesh::edges() const {
488
    return edge(meshEdges.Begin());
489
490
491
}

edge Mesh::edges_end() const {
492
    return edge(meshEdges.End());
493
494
495
}

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

void Mesh::InsertFace(const Point &faceCenter, const Point &cellCenter) {
500
501
502
503
504
505
506
    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));
        }
507
508
509
}

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

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

face Mesh::faces() const {
518
    return face(meshFaces.Begin());
519
520
521
}

face Mesh::faces_end() const {
522
    return face(meshFaces.End());
523
524
525
}

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

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

cell Mesh::overlap() const {
534
    return cell(meshOverlapCells.Begin());
535
536
537
}

cell Mesh::overlap_end() const {
538
    return cell(meshOverlapCells.End());
539
540
541
}

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

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

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

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

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

561
562
563
bnd_face Mesh::bnd_faces_end() const {
    return bnd_face(meshBndFaces.End());
}
564
565

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

569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
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()));
}

586
587
588
589
590
591
592
593
594
595
596
597
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));
}

598
int Mesh::BoundaryFacePart(const Point &bndFaceCenter) const {
599
600
601
602
    auto b = meshBndFaces.Find(bndFaceCenter);
    if (b == meshBndFaces.End())
        return -1;
    return b->second->Part();
603
604
605
}

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

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

612
    refinedMesh.SetDim(dimension);
613
    for (procset p = procsets(); p != procsets_end(); ++p) {
614
        refinedMesh.procSets.Copy(p);
615
616
        if (!RefineEdge(p, refinedMesh))
            RefineFace(p, refinedMesh);
617
    }
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
    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);
635
                                if (dimension == 3) x[8 + k] = project(P);
636
637
638
639
640
641
                                else x[4 + k] = project2D(P);
                            }
                }
                delete t;
            }
            C[i] = refinedMesh.InsertCell(R[i].type(), c.Subdomain(), x);
642
        }
643
644
645
646
647
648
649
650
651
652
653
        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();
654
//    check_mesh(refinedMesh);
655

656
    return meshPtr;
657
658
}

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

666
void Mesh::PrintInfo() const {
667
668
669
670
671
672
673
674
675
    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);
    }
676

677
    mout.PrintInfo(
678
        "Mesh", verbose,
679
680
681
        PrintInfoEntry("Mesh Name", meshName, 1),
        PrintInfoEntry("Distribution", distName, 1),
        PrintInfoEntry("Overlap", olapName, 1),
682
683
684
685
686
687
688
689
690
691
        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),
692
693
        PrintInfoEntry("Min cells on proc", PPM->Min(CellCount(), commSplit), 2),
        PrintInfoEntry("Max cells on proc", PPM->Max(CellCount(), commSplit), 2),
694
        PrintInfoEntry("Communicator split", commSplit, 2),
695
696
697
698
699
700
701
        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)
702
    );
703
704
}

705
706
707
int Mesh::CellCount() const {
    return meshCells.size();
}
708

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

717
718
719
int Mesh::FaceCount() const {
    return meshFaces.size();
}
720

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

729
730
731
int Mesh::EdgeCount() const {
    return meshEdges.size();
}
732

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

741
742
743
int Mesh::VertexCount() const {
    return meshVertices.size();
}
744

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

753
int Mesh::OverlapCount() const {
754
    return meshOverlapCells.size();
755
}
756

757
int Mesh::OverlapCountGeometry() const {
758
    return PPM->SumOnCommSplit(OverlapCount(), commSplit);
759
}
760
761

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

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

770
771
772
int Mesh::ProcSetsCountGeometry() const {
    int count = 0;
    for (procset p = this->procsets(); p != this->procsets_end(); ++p)
773
        if (this->procSets.master(p()))
774
775
776
777
778
779
780
            count++;
    return PPM->SumOnCommSplit(count, commSplit);
}

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


787
788
789