1 #ifdef PARALLEL
2 
3 
4 #include <meshing.hpp>
5 #include "paralleltop.hpp"
6 
7 
8 namespace netgen
9 {
10 
ParallelMeshTopology(const Mesh & amesh)11   ParallelMeshTopology :: ParallelMeshTopology (const Mesh & amesh)
12     : mesh(amesh)
13   {
14     is_updated = false;
15   }
16 
17 
~ParallelMeshTopology()18   ParallelMeshTopology :: ~ParallelMeshTopology ()
19   {
20     ;
21   }
22 
23 
Reset()24   void ParallelMeshTopology :: Reset ()
25   {
26     *testout << "ParallelMeshTopology::Reset" << endl;
27 
28     if ( mesh.GetCommunicator().Size() == 1 ) return;
29 
30     size_t ned = mesh.GetTopology().GetNEdges();
31     size_t nfa = mesh.GetTopology().GetNFaces();
32 
33     if (glob_edge.Size() != ned)
34       {
35 	glob_edge.SetSize(ned);
36 	glob_face.SetSize(nfa);
37 	glob_edge = -1;
38 	glob_face = -1;
39 
40 	loc2distedge.ChangeSize (ned);
41 	loc2distface.ChangeSize (nfa);
42       }
43 
44     if (glob_vert.Size() != mesh.GetNV())
45       {
46 	SetNV(mesh.GetNV());
47 	SetNE(mesh.GetNE());
48       }
49   }
50 
51 
Print() const52   void ParallelMeshTopology ::  Print() const
53   {
54     ;
55   }
56 
57 
EnumeratePointsGlobally()58   void ParallelMeshTopology :: EnumeratePointsGlobally ()
59   {
60     auto comm = mesh.GetCommunicator();
61     auto rank = comm.Rank();
62 
63     size_t oldnv = glob_vert.Size();
64     size_t nv = loc2distvert.Size();
65     *testout << "enumerate globally, loc2distvert.size = " << loc2distvert.Size()
66              << ", glob_vert.size = " << glob_vert.Size() << endl;
67 
68 
69     if (rank == 0)
70       nv = 0;
71 
72     // IntRange newvr(oldnv, nv); // new vertex range
73     auto new_pir = Range(PointIndex(oldnv+PointIndex::BASE),
74                          PointIndex(nv+PointIndex::BASE));
75 
76     glob_vert.SetSize (nv);
77     for (auto pi : new_pir)
78       L2G(pi) = -1;
79 
80     int num_master_points = 0;
81 
82     for (auto pi : new_pir)
83       {
84         auto dps = GetDistantProcs(pi);
85         // check sorted:
86         for (int j = 0; j+1 < dps.Size(); j++)
87           if (dps[j+1] < dps[j]) cout << "wrong sort" << endl;
88 
89         if (dps.Size() == 0 || dps[0] > comm.Rank())
90           L2G(pi) = num_master_points++;
91       }
92 
93     *testout << "nummaster = " << num_master_points << endl;
94 
95     Array<int> first_master_point(comm.Size());
96     comm.AllGather (num_master_points, first_master_point);
97     auto max_oldv = comm.AllReduce (Max (glob_vert.Range(0, oldnv)), MPI_MAX);
98     if (comm.AllReduce (oldnv, MPI_SUM) == 0)
99       max_oldv = PointIndex::BASE-1;
100 
101     size_t num_glob_points = max_oldv+1;
102     for (int i = 0; i < comm.Size(); i++)
103       {
104         int cur = first_master_point[i];
105         first_master_point[i] = num_glob_points;
106         num_glob_points += cur;
107       }
108 
109     for (auto pi : new_pir)
110       if (L2G(pi) != -1)
111         L2G(pi) += first_master_point[comm.Rank()];
112 
113     // ScatterDofData (global_nums);
114 
115     Array<int> nsend(comm.Size()), nrecv(comm.Size());
116     nsend = 0;
117     nrecv = 0;
118 
119     /** Count send/recv size **/
120     for (auto pi : new_pir)
121       if (auto dps = GetDistantProcs(pi); dps.Size())
122         {
123           if (rank < dps[0])
124             for (auto p : dps)
125               nsend[p]++;
126           else
127             nrecv[dps[0]]++;
128         }
129 
130     Table<PointIndex> send_data(nsend);
131     Table<PointIndex> recv_data(nrecv);
132 
133     /** Fill send_data **/
134     nsend = 0;
135     for (auto pi : new_pir)
136       if (auto dps = GetDistantProcs(pi); dps.Size())
137         if (rank < dps[0])
138           for (auto p : dps)
139             send_data[p][nsend[p]++] = L2G(pi);
140 
141     Array<MPI_Request> requests;
142     for (int i = 0; i < comm.Size(); i++)
143       {
144         if (nsend[i])
145           requests.Append (comm.ISend (send_data[i], i, 200));
146         if (nrecv[i])
147           requests.Append (comm.IRecv (recv_data[i], i, 200));
148       }
149 
150     MyMPI_WaitAll (requests);
151 
152     Array<int> cnt(comm.Size());
153     cnt = 0;
154 
155     for (auto pi : new_pir)
156       if (auto dps = GetDistantProcs(pi); dps.Size())
157         if (int master = dps[0]; master < comm.Rank())
158           L2G(pi) = recv_data[master][cnt[master]++];
159 
160     // reorder following global ordering:
161     Array<int> index0(glob_vert.Size());
162     for (int pi : Range(index0))
163       index0[pi] = pi;
164     QuickSortI (glob_vert, index0);
165 
166     if (rank != 0)
167       {
168         Array<PointIndex, PointIndex> inv_index(index0.Size());
169         for (int i = 0; i < index0.Size(); i++)
170           inv_index[index0[i]+PointIndex::BASE] = i+PointIndex::BASE;
171 
172         for (auto & el : mesh.VolumeElements())
173           for (PointIndex & pi : el.PNums())
174             pi = inv_index[pi];
175         for (auto & el : mesh.SurfaceElements())
176           for (PointIndex & pi : el.PNums())
177             pi = inv_index[pi];
178         for (auto & el : mesh.LineSegments())
179           for (PointIndex & pi : el.PNums())
180             pi = inv_index[pi];
181 
182         // auto hpoints (mesh.Points());
183         Array<MeshPoint, PointIndex> hpoints { mesh.Points() };
184         for (PointIndex pi : Range(mesh.Points()))
185           mesh.Points()[inv_index[pi]] = hpoints[pi];
186 
187         if (mesh.mlbetweennodes.Size() == mesh.Points().Size())
188           {
189             NgArray<PointIndices<2>,PointIndex::BASE> hml { mesh.mlbetweennodes };
190             for (PointIndex pi : Range(mesh.Points()))
191               mesh.mlbetweennodes[inv_index[pi]] = hml[pi];
192           }
193 
194         // *testout << "index0 = " << endl << index0 << endl;
195         // *testout << "loc2distvertold = " << endl;
196         // for (auto i : Range(index0))
197         // *testout << "l " << i << " globi "<< glob_vert[i]  << " dist = " << loc2distvert[i] << endl;
198 
199         DynamicTable<int> oldtable = std::move(loc2distvert);
200         loc2distvert = DynamicTable<int> (oldtable.Size());
201         for (size_t i = 0; i < oldtable.Size(); i++)
202           for (auto val : oldtable[index0[i]])
203             loc2distvert.Add (i, val);
204 
205         Array<int> hglob_vert(glob_vert);
206         for (int i = 0; i < index0.Size(); i++)
207           glob_vert[i] = hglob_vert[index0[i]];
208 
209         // *testout << "loc2distvertnew = " << endl;
210         // for (auto i : Range(index0))
211         // *testout << "l " << i << " globi "<< glob_vert[i]  << " dist = " << loc2distvert[i] << endl;
212       }
213 
214     /*
215     for (size_t i = 0; i+1 < glob_vert.Size(); i++)
216       if (glob_vert[i] > glob_vert[i+1])
217         cout << "wrong ordering of globvert" << endl;
218     */
219     if (glob_vert.Size() > 1)
220       for (auto i : Range(glob_vert).Modify(0,-1))
221         if (glob_vert[i] > glob_vert[i+1])
222           cout << "wrong ordering of globvert" << endl;
223   }
224 
225   /*
226   void ParallelMeshTopology :: SetDistantFaceNum (int dest, int locnum)
227   {
228     for ( int i = 0; i < loc2distface[locnum-1].Size(); i+=1 )
229       if ( loc2distface[locnum-1][i] == dest )
230 	return;
231     loc2distface.Add(locnum-1, dest);
232   }
233 
234   void ParallelMeshTopology :: SetDistantPNum (int dest, int locnum)
235   {
236     for ( int i = 0;  i < loc2distvert[locnum-1].Size(); i+=1 )
237       if ( loc2distvert[locnum-1][i] == dest )
238 	return;
239     loc2distvert.Add (locnum-1, dest);
240   }
241 
242 
243   void ParallelMeshTopology :: SetDistantEdgeNum (int dest, int locnum)
244   {
245     for ( int i = 0; i < loc2distedge[locnum-1].Size(); i+=1 )
246       if ( loc2distedge[locnum-1][i] == dest )
247 	return;
248     loc2distedge.Add (locnum-1, dest);
249   }
250   */
251 
SetNV_Loc2Glob(int anv)252   void ParallelMeshTopology :: SetNV_Loc2Glob (int anv)
253   {
254     glob_vert.SetSize(anv);
255     glob_vert = -1;
256   }
257 
SetNV(int anv)258   void ParallelMeshTopology :: SetNV (int anv)
259   {
260     // glob_vert.SetSize(anv);
261     // glob_vert = -1;
262     // loc2distvert.ChangeSize (anv);
263 
264     DynamicTable<int> oldtable(loc2distvert.Size());
265     for (size_t i = 0; i < loc2distvert.Size(); i++)
266       for (auto val : loc2distvert[i])
267         oldtable.Add (i, val);
268     loc2distvert = DynamicTable<int> (anv);
269     for (size_t i = 0; i < min(size_t(anv), oldtable.Size()); i++)
270       for (auto val : oldtable[i])
271         loc2distvert.Add (i, val);
272   }
273 
SetNE(int ane)274   void ParallelMeshTopology :: SetNE ( int ane )
275   {
276     glob_el.SetSize (ane);
277     glob_el = -1;
278   }
279 
SetNSE(int anse)280   void ParallelMeshTopology :: SetNSE ( int anse )
281   {
282     glob_surfel.SetSize(anse);
283     glob_surfel = -1;
284   }
285 
SetNSegm(int anseg)286   void ParallelMeshTopology :: SetNSegm ( int anseg )
287   {
288     glob_segm.SetSize (anseg);
289     glob_segm = -1;
290   }
291 
292 
293 
294 
295 
296 
UpdateCoarseGridGlobal()297   void ParallelMeshTopology :: UpdateCoarseGridGlobal ()
298   {
299     // cout << "updatecoarsegridglobal called" << endl;
300     if (id == 0)
301       PrintMessage ( 3, "UPDATE GLOBAL COARSEGRID STARTS" );
302 
303     int timer = NgProfiler::CreateTimer ("UpdateCoarseGridGlobal");
304     NgProfiler::RegionTimer reg(timer);
305 
306     *testout << "ParallelMeshTopology :: UpdateCoarseGridGlobal" << endl;
307 
308     const MeshTopology & topology = mesh.GetTopology();
309     auto comm = mesh.GetCommunicator();
310 
311     if ( id == 0 )
312       {
313 	NgArray<NgArray<int>*> sendarrays(ntasks);
314 	for (int dest = 1; dest < ntasks; dest++)
315 	  sendarrays[dest] = new NgArray<int>;
316 
317 	NgArray<int> edges, faces;
318 	for (int el = 1; el <= mesh.GetNE(); el++)
319 	  {
320 	    topology.GetElementFaces (el, faces);
321 	    topology.GetElementEdges (el, edges);
322 	    const Element & volel = mesh.VolumeElement (el);
323 
324 	    // NgArray<int> & sendarray = *sendarrays[volel.GetPartition()];
325             NgArray<int> & sendarray = *sendarrays[mesh.vol_partition[el-1]];
326 
327 	    for ( int i = 0; i < edges.Size(); i++ )
328 	      sendarray.Append (edges[i]);
329 	    for ( int i = 0; i < faces.Size(); i++ )
330 	      sendarray.Append (faces[i]);
331 	  }
332 
333 	for (int el = 1; el <= mesh.GetNSE(); el++)
334 	  {
335 	    topology.GetSurfaceElementEdges (el, edges);
336 	    const Element2d & surfel = mesh.SurfaceElement (el);
337 	    // NgArray<int> & sendarray = *sendarrays[surfel.GetPartition()];
338             NgArray<int> & sendarray = *sendarrays[mesh.surf_partition[el-1]];
339 
340 	    for ( int i = 0; i < edges.Size(); i++ )
341 	      sendarray.Append (edges[i]);
342 	    sendarray.Append (topology.GetSurfaceElementFace (el));
343 	  }
344 
345 	Array<MPI_Request> sendrequests;
346 	for (int dest = 1; dest < ntasks; dest++)
347 	  // sendrequests.Append (MyMPI_ISend (*sendarrays[dest], dest, MPI_TAG_MESH+10, comm));
348           sendrequests.Append (comm.ISend (FlatArray<int>(*sendarrays[dest]), dest, MPI_TAG_MESH+10));
349 	MyMPI_WaitAll (sendrequests);
350 
351 	for (int dest = 1; dest < ntasks; dest++)
352 	  delete sendarrays[dest];
353       }
354 
355     else
356 
357       {
358 	NgArray<int> recvarray;
359 	MyMPI_Recv (recvarray, 0, MPI_TAG_MESH+10, comm);
360 
361 	int ii = 0;
362 
363 	NgArray<int> faces, edges;
364 
365 	for (int volel = 1; volel <= mesh.GetNE(); volel++)
366 	  {
367 	    topology.GetElementEdges ( volel, edges);
368 	    for ( int i = 0; i  < edges.Size(); i++)
369 	      SetLoc2Glob_Edge ( edges[i], recvarray[ii++]);
370 
371 	    topology.GetElementFaces( volel, faces);
372 	    for ( int i = 0; i  < faces.Size(); i++)
373 	      SetLoc2Glob_Face ( faces[i], recvarray[ii++]);
374 	  }
375 
376 	for (int surfel = 1; surfel <= mesh.GetNSE(); surfel++)
377 	  {
378 	    topology.GetSurfaceElementEdges (surfel, edges);
379 	    for (int i = 0; i  < edges.Size(); i++)
380 	      SetLoc2Glob_Edge (edges[i], recvarray[ii++]);
381 	    int face = topology.GetSurfaceElementFace (surfel);
382 	    SetLoc2Glob_Face ( face, recvarray[ii++]);
383 	  }
384       }
385 
386     is_updated = true;
387   }
388 
389 
IdentifyVerticesAfterRefinement()390   void ParallelMeshTopology :: IdentifyVerticesAfterRefinement()
391   {
392     static Timer t("ParallelTopology::UpdateCoarseGrid"); RegionTimer r(t);
393 
394     NgMPI_Comm comm = mesh.GetCommunicator();
395     int id = comm.Rank();
396     int ntasks = comm.Size();
397 
398     if (ntasks == 1) return;
399 
400     Reset();
401     static int timer = NgProfiler::CreateTimer ("UpdateCoarseGrid");
402     NgProfiler::RegionTimer reg(timer);
403 
404 
405     (*testout) << "UPDATE COARSE GRID PARALLEL TOPOLOGY " << endl;
406     if (id == 0)
407       PrintMessage (1, "update parallel topology");
408 
409 
410     const MeshTopology & topology = mesh.GetTopology();
411 
412     NgArray<int> cnt_send(ntasks);
413 
414     int maxsize = comm.AllReduce (mesh.mlbetweennodes.Size(), MPI_MAX);
415     // update new vertices after mesh-refinement
416     if (maxsize > 0)
417       {
418         int newnv = mesh.mlbetweennodes.Size();
419 
420         loc2distvert.ChangeSize(mesh.mlbetweennodes.Size());
421 
422 	bool changed = true;
423 	while (changed)
424 	  {
425 	    changed = false;
426 
427 	    // build exchange vertices
428 	    cnt_send = 0;
429 	    for (PointIndex pi : mesh.Points().Range())
430 	      for (int dist : GetDistantProcs(pi))
431 		cnt_send[dist]++;
432 	    TABLE<int> dest2vert(cnt_send);
433 	    for (PointIndex pi : mesh.Points().Range())
434 	      for (int dist : GetDistantProcs(pi))
435 		dest2vert.Add (dist, pi);
436 
437 	    for (PointIndex pi = PointIndex::BASE; pi < newnv+PointIndex::BASE; pi++)
438               if (auto [v1,v2] = mesh.mlbetweennodes[pi]; v1.IsValid())
439                 {
440                   auto procs1 = GetDistantProcs(v1);
441                   auto procs2 = GetDistantProcs(v2);
442                   for (int p : procs1)
443                     if (procs2.Contains(p))
444                       cnt_send[p]++;
445                 }
446 
447 	    TABLE<int> dest2pair(cnt_send);
448 
449             for (PointIndex pi : mesh.mlbetweennodes.Range())
450               if (auto [v1,v2] = mesh.mlbetweennodes[pi]; v1.IsValid())
451                 {
452                   auto procs1 = GetDistantProcs(v1);
453                   auto procs2 = GetDistantProcs(v2);
454                   for (int p : procs1)
455                     if (procs2.Contains(p))
456                       dest2pair.Add (p, pi);
457                 }
458 
459 	    cnt_send = 0;
460 	    for (PointIndex pi : mesh.mlbetweennodes.Range())
461               if (auto [v1,v2] = mesh.mlbetweennodes[pi]; v1.IsValid())
462                 {
463                   auto procs1 = GetDistantProcs(v1);
464                   auto procs2 = GetDistantProcs(v2);
465 
466                   for (int p : procs1)
467                     if (procs2.Contains(p))
468                       cnt_send[p]+=2;
469                 }
470 
471 	    TABLE<int> send_verts(cnt_send);
472 
473 	    NgArray<int, PointIndex::BASE> loc2exchange(mesh.GetNV());
474 
475 	    for (int dest = 0; dest < ntasks; dest++)
476 	      if (dest != id)
477 		{
478 		  loc2exchange = -1;
479 		  int cnt = 0;
480 		  for (PointIndex pi : dest2vert[dest])
481 		    loc2exchange[pi] = cnt++;
482 
483 		  for (PointIndex pi : dest2pair[dest])
484                     if (auto [v1,v2] = mesh.mlbetweennodes[pi]; v1.IsValid())
485                       {
486                         auto procs1 = GetDistantProcs(v1);
487                         auto procs2 = GetDistantProcs(v2);
488 
489                         if (procs1.Contains(dest) && procs2.Contains(dest))
490                           {
491                             send_verts.Add (dest, loc2exchange[v1]);
492                             send_verts.Add (dest, loc2exchange[v2]);
493                           }
494                       }
495 		}
496 
497 	    TABLE<int> recv_verts(ntasks);
498 	    MyMPI_ExchangeTable (send_verts, recv_verts, MPI_TAG_MESH+9, comm);
499 
500 	    for (int dest = 0; dest < ntasks; dest++)
501 	      if (dest != id)
502 		{
503 		  loc2exchange = -1;
504 		  int cnt = 0;
505 
506 		  for (PointIndex pi : dest2vert[dest])
507 		    loc2exchange[pi] = cnt++;
508 
509 		  NgFlatArray<int> recvarray = recv_verts[dest];
510 		  for (int ii = 0; ii < recvarray.Size(); ii+=2)
511 		    for (PointIndex pi : dest2pair[dest])
512 		      {
513 			PointIndex v1 = mesh.mlbetweennodes[pi][0];
514 			PointIndex v2 = mesh.mlbetweennodes[pi][1];
515 			if (v1.IsValid())
516 			  {
517 			    INDEX_2 re(recvarray[ii], recvarray[ii+1]);
518 			    INDEX_2 es(loc2exchange[v1], loc2exchange[v2]);
519 			    // if (es == re && !IsExchangeVert(dest, pi))
520                             if (es == re && !GetDistantProcs(pi).Contains(dest))
521 			      {
522 				// SetDistantPNum(dest, pi);
523                                 AddDistantProc (pi, dest);
524 				changed = true;
525 			      }
526 			  }
527 		      }
528 		}
529 
530             changed = comm.AllReduce (changed, MPI_LOR);
531 	  }
532       }
533 
534     NgArray<int> sendarray, recvarray;
535     // cout << "UpdateCoarseGrid - edges" << endl;
536 
537     // static int timerv = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex vertices");
538     static int timere = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex edges");
539     static int timerf = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex faces");
540 
541 
542     NgProfiler::StartTimer (timere);
543 
544     // build exchange vertices
545     cnt_send = 0;
546     for (PointIndex pi : mesh.Points().Range())
547       for (int dist : GetDistantProcs(pi))
548 	cnt_send[dist]++;
549     TABLE<int> dest2vert(cnt_send);
550     for (PointIndex pi : mesh.Points().Range())
551       for (int dist : GetDistantProcs(pi))
552 	dest2vert.Add (dist, pi);
553 
554     // MPI_Group_free(&MPI_LocalGroup);
555     // MPI_Comm_free(&MPI_LocalComm);
556   }
557 
558 
UpdateCoarseGrid()559   void ParallelMeshTopology :: UpdateCoarseGrid ()
560   {
561     static Timer t("ParallelTopology::UpdateCoarseGrid"); RegionTimer r(t);
562     // cout << "UpdateCoarseGrid" << endl;
563     // if (is_updated) return;
564 
565     NgMPI_Comm comm = mesh.GetCommunicator();
566     int id = comm.Rank();
567     int ntasks = comm.Size();
568 
569     if (ntasks == 1) return;
570 
571     Reset();
572     static int timer = NgProfiler::CreateTimer ("UpdateCoarseGrid");
573     NgProfiler::RegionTimer reg(timer);
574 
575 
576     (*testout) << "UPDATE COARSE GRID PARALLEL TOPOLOGY " << endl;
577     if (id == 0)
578       PrintMessage (1, "update parallel topology");
579 
580 
581     // UpdateCoarseGridGlobal();
582 
583 
584     /*
585     // MPI_Barrier (MPI_COMM_WORLD);
586 
587     MPI_Group MPI_GROUP_comm;
588     MPI_Group MPI_LocalGroup;
589     MPI_Comm MPI_LocalComm;
590 
591     int process_ranks[] = { 0 };
592     MPI_Comm_group (comm, &MPI_GROUP_comm);
593     MPI_Group_excl (MPI_GROUP_comm, 1, process_ranks, &MPI_LocalGroup);
594     MPI_Comm_create (comm, MPI_LocalGroup, &MPI_LocalComm);
595 
596     if (id == 0)
597       {
598         // SetNV(0);
599         // EnumeratePointsGlobally();
600         return;
601       }
602     */
603 
604     const MeshTopology & topology = mesh.GetTopology();
605 
606     NgArray<int> cnt_send(ntasks);
607 
608     // NgArray<int> sendarray, recvarray;
609     // cout << "UpdateCoarseGrid - edges" << endl;
610 
611     // static int timerv = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex vertices");
612     static int timere = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex edges");
613     static int timerf = NgProfiler::CreateTimer ("UpdateCoarseGrid - ex faces");
614 
615 
616     NgProfiler::StartTimer (timere);
617 
618 
619     int nfa = topology . GetNFaces();
620     int ned = topology . GetNEdges();
621 
622     // build exchange vertices
623     cnt_send = 0;
624     for (PointIndex pi : mesh.Points().Range())
625       for (int dist : GetDistantProcs(pi))
626 	cnt_send[dist]++;
627     TABLE<int> dest2vert(cnt_send);
628     for (PointIndex pi : mesh.Points().Range())
629       for (int dist : GetDistantProcs(pi))
630 	dest2vert.Add (dist, pi);
631 
632     // exchange edges
633     cnt_send = 0;
634     int v1, v2;
635     for (int edge = 1; edge <= ned; edge++)
636       {
637 	topology.GetEdgeVertices (edge, v1, v2);
638         /*
639 	for (int dest = 1; dest < ntasks; dest++)
640 	  // if (IsExchangeVert (dest, v1) && IsExchangeVert (dest, v2))
641           if (GetDistantProcs(v1).Contains(dest) && GetDistantProcs(v2).Contains(dest))
642 	    cnt_send[dest-1]+=1;
643         */
644         for (auto p : GetDistantProcs(v1))
645           if (GetDistantProcs(v2).Contains(p))
646 	    cnt_send[p]+=1;
647       }
648 
649     TABLE<int> dest2edge(cnt_send);
650     for (int & v : cnt_send) v *= 2;
651     TABLE<int> send_edges(cnt_send);
652 
653     for (int edge = 1; edge <= ned; edge++)
654       {
655 	topology.GetEdgeVertices (edge, v1, v2);
656 	for (int dest = 0; dest < ntasks; dest++)
657 	  // if (IsExchangeVert (dest, v1) && IsExchangeVert (dest, v2))
658           if (GetDistantProcs(v1).Contains(dest) && GetDistantProcs(v2).Contains(dest))
659 	    dest2edge.Add (dest, edge);
660       }
661 
662 
663     NgArray<int, PointIndex::BASE> loc2exchange(mesh.GetNV());
664     for (int dest = 0; dest < ntasks; dest++)
665       {
666         loc2exchange = -1;
667         int cnt = 0;
668         for (PointIndex pi : dest2vert[dest])
669 	  loc2exchange[pi] = cnt++;
670 
671 	for (int edge : dest2edge[dest])
672           {
673             topology.GetEdgeVertices (edge, v1, v2);
674             // if (IsExchangeVert (dest, v1) && IsExchangeVert (dest, v2))
675             if (GetDistantProcs(v1).Contains(dest) && GetDistantProcs(v2).Contains(dest))
676               {
677                 send_edges.Add (dest, loc2exchange[v1]);
678                 send_edges.Add (dest, loc2exchange[v2]);
679               }
680           }
681       }
682 
683     // cout << "UpdateCoarseGrid - edges mpi-exchange" << endl;
684     TABLE<int> recv_edges(ntasks);
685     MyMPI_ExchangeTable (send_edges, recv_edges, MPI_TAG_MESH+9, comm);
686     // cout << "UpdateCoarseGrid - edges mpi-exchange done" << endl;
687 
688     for (int dest = 0; dest < ntasks; dest++)
689       {
690 	auto ex2loc = dest2vert[dest];
691 	if (ex2loc.Size() == 0) continue;
692 
693 	INDEX_2_CLOSED_HASHTABLE<int> vert2edge(2*dest2edge[dest].Size()+10);
694 	for (int edge : dest2edge[dest])
695 	  {
696 	    topology.GetEdgeVertices (edge, v1, v2);
697 	    vert2edge.Set(INDEX_2(v1,v2), edge);
698 	  }
699 
700 	NgFlatArray<int> recvarray = recv_edges[dest];
701         for (int ii = 0; ii < recvarray.Size(); ii+=2)
702 	  {
703 	    INDEX_2 re(ex2loc[recvarray[ii]],
704 		       ex2loc[recvarray[ii+1]]);
705 	    if (vert2edge.Used(re))
706 	      // SetDistantEdgeNum(dest, vert2edge.Get(re));
707               AddDistantEdgeProc(vert2edge.Get(re)-1, dest);
708 	  }
709       }
710 
711 
712 
713     NgProfiler::StopTimer (timere);
714 
715     // cout << "UpdateCoarseGrid - faces" << endl;
716     if (mesh.GetDimension() == 3)
717       {
718 	NgProfiler::StartTimer (timerf);
719 	NgArray<int> verts;
720 
721 	// exchange faces
722 	cnt_send = 0;
723 	for (int face = 1; face <= nfa; face++)
724 	  {
725 	    topology.GetFaceVertices (face, verts);
726 	    for (int dest = 0; dest < ntasks; dest++)
727 	      if (dest != id)
728                 /*
729 		if (IsExchangeVert (dest, verts[0]) &&
730 		    IsExchangeVert (dest, verts[1]) &&
731 		    IsExchangeVert (dest, verts[2]))
732                 */
733                 if (GetDistantProcs (verts[0]).Contains(dest) &&
734                     GetDistantProcs (verts[1]).Contains(dest) &&
735                     GetDistantProcs (verts[2]).Contains(dest))
736 		  cnt_send[dest]++;
737 	  }
738 
739 	TABLE<int> dest2face(cnt_send);
740 	for (int face = 1; face <= nfa; face++)
741 	  {
742 	    topology.GetFaceVertices (face, verts);
743 	    for (int dest = 0; dest < ntasks; dest++)
744 	      if (dest != id)
745                 /*
746 		if (IsExchangeVert (dest, verts[0]) &&
747 		    IsExchangeVert (dest, verts[1]) &&
748 		    IsExchangeVert (dest, verts[2]))
749                 */
750                 if (GetDistantProcs (verts[0]).Contains(dest) &&
751                     GetDistantProcs (verts[1]).Contains(dest) &&
752                     GetDistantProcs (verts[2]).Contains(dest))
753 		  dest2face.Add(dest, face);
754 	  }
755 
756 	for (int & c : cnt_send) c*=3;
757 	TABLE<int> send_faces(cnt_send);
758 	NgArray<int, PointIndex::BASE> loc2exchange(mesh.GetNV());
759 	for (int dest = 0; dest < ntasks; dest++)
760 	  if (dest != id)
761 	    {
762 	      if (dest2vert[dest].Size() == 0) continue;
763 
764 	      loc2exchange = -1;
765 	      int cnt = 0;
766 	      for (PointIndex pi : dest2vert[dest])
767 		loc2exchange[pi] = cnt++;
768 
769 	      for (int face : dest2face[dest])
770 		{
771 		  topology.GetFaceVertices (face, verts);
772                   /*
773 		  if (IsExchangeVert (dest, verts[0]) &&
774 		      IsExchangeVert (dest, verts[1]) &&
775 		      IsExchangeVert (dest, verts[2]))
776                   */
777                   if (GetDistantProcs (verts[0]).Contains(dest) &&
778                       GetDistantProcs (verts[1]).Contains(dest) &&
779                       GetDistantProcs (verts[2]).Contains(dest))
780 		    {
781 		      send_faces.Add (dest, loc2exchange[verts[0]]);
782 		      send_faces.Add (dest, loc2exchange[verts[1]]);
783 		      send_faces.Add (dest, loc2exchange[verts[2]]);
784 		    }
785 		}
786 	    }
787 
788 	// cout << "UpdateCoarseGrid - faces mpi-exchange" << endl;
789 	TABLE<int> recv_faces(ntasks);
790 	MyMPI_ExchangeTable (send_faces, recv_faces, MPI_TAG_MESH+9, comm);
791 	// cout << "UpdateCoarseGrid - faces mpi-exchange done" << endl;
792 
793 	for (int dest = 0; dest < ntasks; dest++)
794 	  {
795 	    auto ex2loc = dest2vert[dest];
796 	    if (ex2loc.Size() == 0) continue;
797 
798 	    INDEX_3_CLOSED_HASHTABLE<int> vert2face(2*dest2face[dest].Size()+10);
799 	    for (int face : dest2face[dest])
800 	      {
801 		topology.GetFaceVertices (face, verts);
802 		vert2face.Set(INDEX_3(verts[0], verts[1], verts[2]), face);
803 	      }
804 
805 	    NgFlatArray<int> recvarray = recv_faces[dest];
806 	    for (int ii = 0; ii < recvarray.Size(); ii+=3)
807 	      {
808 		INDEX_3 re(ex2loc[recvarray[ii]],
809 			   ex2loc[recvarray[ii+1]],
810 			   ex2loc[recvarray[ii+2]]);
811 		if (vert2face.Used(re))
812 		  AddDistantFaceProc(vert2face.Get(re)-1, dest);
813 	      }
814 	  }
815 
816 	NgProfiler::StopTimer (timerf);
817       }
818     // cout << "UpdateCoarseGrid - done" << endl;
819     // EnumeratePointsGlobally();
820     is_updated = true;
821 
822     // MPI_Group_free(&MPI_LocalGroup);
823     // MPI_Comm_free(&MPI_LocalComm);
824   }
825 }
826 
827 
828 #endif
829