faMeshDecomposition.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2016-2017 Wikki Ltd
9  Copyright (C) 2018-2019 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "faMeshDecomposition.H"
30 #include "Time.H"
31 #include "dictionary.H"
32 #include "labelIOList.H"
33 #include "processorFaPatch.H"
34 #include "faMesh.H"
35 #include "OSspecific.H"
36 #include "Map.H"
37 #include "globalMeshData.H"
38 #include "decompositionModel.H"
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void Foam::faMeshDecomposition::distributeFaces()
43 {
44  Info<< "\nCalculating distribution of faces" << endl;
45 
46  cpuTime decompositionTime;
47 
48  for (label procI = 0; procI < nProcs(); procI++)
49  {
50  Time processorDb
51  (
53  time().rootPath(),
54  time().caseName()/("processor" + Foam::name(procI))
55  );
56 
57  fvMesh procMesh
58  (
59  IOobject
60  (
62  processorDb.timeName(),
63  processorDb
64  )
65  );
66 
67  // If faMesh's fvPatch is a part of the global face zones, faces of that
68  // patch will be present on all processors. Because of that, looping
69  // through faceProcAddressing will decompose global faMesh faces to the
70  // very last processor regardless of where fvPatch is really decomposed.
71  // Since global faces which do not belong to specific processor are
72  // located at the end of the faceProcAddressing, cutting it at
73  // i = owner.size() will correctly decompose faMesh faces.
74  // Vanja Skuric, 2016-04-21
75  if (hasGlobalFaceZones_)
76  {
78  (
80  (
81  IOobject
82  (
83  "faceProcAddressing",
84  "constant",
85  procMesh.meshSubDir,
86  procMesh,
89  )
90  )
91  );
92 
93  const label ownerSize =
94  (
96  (
97  IOobject
98  (
99  "owner",
100  "constant",
101  procMesh.meshSubDir,
102  procMesh,
105  )
106  )
107  ).size();
108 
109  labelHashSet faceProcAddressingHash(ownerSize);
110 
111  for (int i = 0; i < ownerSize; ++i)
112  {
113  faceProcAddressingHash.insert(faceProcAddressing[i]);
114  }
115 
116  forAll(faceLabels(), faceI)
117  {
118  if (faceProcAddressingHash.found(faceLabels()[faceI] + 1))
119  {
120  faceToProc_[faceI] = procI;
121  }
122  }
123  }
124  else
125  {
126  labelHashSet faceProcAddressingHash
127  (
129  (
130  IOobject
131  (
132  "faceProcAddressing",
133  "constant",
134  procMesh.meshSubDir,
135  procMesh,
138  )
139  )
140  );
141 
142  forAll(faceLabels(), faceI)
143  {
144  if (faceProcAddressingHash.found(faceLabels()[faceI] + 1))
145  {
146  faceToProc_[faceI] = procI;
147  }
148  }
149  }
150  }
151 
152  Info<< "\nFinished decomposition in "
153  << decompositionTime.elapsedCpuTime()
154  << " s" << endl;
155 }
156 
157 
158 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
159 
161 (
162  const fvMesh& mesh,
163  const fileName& decompDictFile
164 )
165 :
166  faMesh(mesh),
167  decompDictFile_(decompDictFile),
168  nProcs_
169  (
170  decompositionMethod::nDomains
171  (
172  decompositionModel::New
173  (
174  mesh,
175  decompDictFile
176  )
177  )
178  ),
179  distributed_(false),
180  hasGlobalFaceZones_(false),
181  faceToProc_(nFaces()),
182  procFaceLabels_(nProcs_),
183  procMeshEdgesMap_(nProcs_),
184  procNInternalEdges_(nProcs_, Zero),
185  procPatchEdgeLabels_(nProcs_),
186  procPatchPointAddressing_(nProcs_),
187  procPatchEdgeAddressing_(nProcs_),
188  procEdgeAddressing_(nProcs_),
189  procFaceAddressing_(nProcs_),
190  procBoundaryAddressing_(nProcs_),
191  procPatchSize_(nProcs_),
192  procPatchStartIndex_(nProcs_),
193  procNeighbourProcessors_(nProcs_),
194  procProcessorPatchSize_(nProcs_),
195  procProcessorPatchStartIndex_(nProcs_),
196  globallySharedPoints_(0),
197  cyclicParallel_(false)
198 {
199  const decompositionModel& model = decompositionModel::New
200  (
201  mesh,
202  decompDictFile
203  );
204 
205  model.readIfPresent("distributed", distributed_);
206  hasGlobalFaceZones_ = model.found("globalFaceZones");
207 }
208 
209 
210 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
211 
213 {
214  // Decide which cell goes to which processor
215  distributeFaces();
216 
217  Info<< "\nDistributing faces to processors" << endl;
218 
219  // Memory management
220  {
221  List<SLList<label>> procFaceList(nProcs());
222 
223  forAll(faceToProc_, faceI)
224  {
225  if (faceToProc_[faceI] >= nProcs())
226  {
227  FatalErrorIn("Finite area mesh decomposition")
228  << "Impossible processor label " << faceToProc_[faceI]
229  << "for face " << faceI
230  << abort(FatalError);
231  }
232  else
233  {
234  procFaceList[faceToProc_[faceI]].append(faceI);
235  }
236  }
237 
238  // Convert linked lists into normal lists
239  forAll(procFaceList, procI)
240  {
241  procFaceAddressing_[procI] = procFaceList[procI];
242  }
243  }
244 
245 
246  // Find processor mesh faceLabels and ...
247 
248  for (label procI = 0; procI < nProcs(); procI++)
249  {
250  Time processorDb
251  (
253  time().rootPath(),
254  time().caseName()/("processor" + Foam::name(procI))
255  );
256 
257  fvMesh procFvMesh
258  (
259  IOobject
260  (
262  processorDb.timeName(),
263  processorDb
264  )
265  );
266 
267  labelIOList fvPointProcAddressing
268  (
269  IOobject
270  (
271  "pointProcAddressing",
272  "constant",
273  procFvMesh.meshSubDir,
274  procFvMesh,
277  )
278  );
279 
280  Map<label> fvFaceProcAddressingHash;
281 
282  {
283  labelIOList fvFaceProcAddressing
284  (
285  IOobject
286  (
287  "faceProcAddressing",
288  "constant",
289  procFvMesh.meshSubDir,
290  procFvMesh,
293  )
294  );
295 
296  forAll(fvFaceProcAddressing, faceI)
297  {
298  fvFaceProcAddressingHash.insert
299  (
300  fvFaceProcAddressing[faceI], faceI
301  );
302  }
303  };
304 
305  const labelList& curProcFaceAddressing = procFaceAddressing_[procI];
306 
307  labelList& curFaceLabels = procFaceLabels_[procI];
308 
309  curFaceLabels = labelList(curProcFaceAddressing.size(), -1);
310 
311  forAll(curProcFaceAddressing, faceI)
312  {
313  curFaceLabels[faceI] =
314  fvFaceProcAddressingHash.find
315  (
316  faceLabels()[curProcFaceAddressing[faceI]] + 1
317  )();
318  }
319 
320  // create processor finite area mesh
321  faMesh procMesh
322  (
323  procFvMesh,
324  procFaceLabels_[procI]
325  );
326 
327  const indirectPrimitivePatch& patch = this->patch();
328  const Map<label>& map = patch.meshPointMap();
329 
330  EdgeMap<label> edgesHash;
331 
332  label edgeI = -1;
333 
334  const label nIntEdges = patch.nInternalEdges();
335 
336  for (label curEdge = 0; curEdge < nIntEdges; curEdge++)
337  {
338  edgesHash.insert(patch.edges()[curEdge], ++edgeI);
339  }
340 
341  forAll(boundary(), patchI)
342  {
343  // Include emptyFaPatch
344 
345  const label size = boundary()[patchI].labelList::size();
346 
347  for(int eI=0; eI<size; eI++)
348  {
349  edgesHash.insert
350  (
351  patch.edges()[boundary()[patchI][eI]],
352  ++edgeI
353  );
354  }
355  }
356 
357 
358  const indirectPrimitivePatch& procPatch = procMesh.patch();
359  const vectorField& procPoints = procPatch.localPoints();
360  const labelList& procMeshPoints = procPatch.meshPoints();
361  const edgeList& procEdges = procPatch.edges();
362 
363  labelList& curPatchPointAddressing = procPatchPointAddressing_[procI];
364  curPatchPointAddressing.setSize(procPoints.size(), -1);
365 
366  forAll(procPoints, pointI)
367  {
368  curPatchPointAddressing[pointI] =
369  map[fvPointProcAddressing[procMeshPoints[pointI]]];
370  }
371 
372  labelList& curPatchEdgeAddressing = procPatchEdgeAddressing_[procI];
373  curPatchEdgeAddressing.setSize(procEdges.size(), -1);
374 
375  forAll(procEdges, edgeI)
376  {
377  edge curGlobalEdge = procEdges[edgeI];
378  curGlobalEdge[0] = curPatchPointAddressing[curGlobalEdge[0]];
379  curGlobalEdge[1] = curPatchPointAddressing[curGlobalEdge[1]];
380 
381  curPatchEdgeAddressing[edgeI] = edgesHash.find(curGlobalEdge)();
382  }
383 
384  Map<label>& curMap = procMeshEdgesMap_[procI];
385 
386  curMap = Map<label>(2*procEdges.size());
387 
388  forAll(curPatchEdgeAddressing, edgeI)
389  {
390  curMap.insert(curPatchEdgeAddressing[edgeI], edgeI);
391  }
392 
393  procNInternalEdges_[procI] = procPatch.nInternalEdges();
394  }
395 
396 
397  Info << "\nDistributing edges to processors" << endl;
398 
399  // Loop through all internal edges and decide which processor they
400  // belong to. First visit all internal edges.
401 
402  // set references to the original mesh
403  const faBoundaryMesh& patches = boundary();
404  const edgeList& edges = this->edges();
405  const labelList& owner = edgeOwner();
406  const labelList& neighbour = edgeNeighbour();
407 
408  // Memory management
409  {
410  List<SLList<label>> procEdgeList(nProcs());
411 
412  forAll(procEdgeList, procI)
413  {
414  for(label i=0; i<procNInternalEdges_[procI]; i++)
415  {
416  procEdgeList[procI].append(procPatchEdgeAddressing_[procI][i]);
417  }
418  }
419 
420 
421  // Detect inter-processor boundaries
422 
423  // Neighbour processor for each subdomain
424  List<SLList<label>> interProcBoundaries(nProcs());
425 
426  // Edge labels belonging to each inter-processor boundary
427  List<SLList<SLList<label>>> interProcBEdges(nProcs());
428 
429  List<SLList<label>> procPatchIndex(nProcs());
430 
431  forAll(neighbour, edgeI)
432  {
433  if (faceToProc_[owner[edgeI]] != faceToProc_[neighbour[edgeI]])
434  {
435  // inter - processor patch edge found. Go through the list of
436  // inside boundaries for the owner processor and try to find
437  // this inter-processor patch.
438 
439  bool interProcBouFound = false;
440 
441  const label ownProc = faceToProc_[owner[edgeI]];
442  const label neiProc = faceToProc_[neighbour[edgeI]];
443 
444  auto curInterProcBdrsOwnIter =
445  interProcBoundaries[ownProc].cbegin();
446 
447  auto curInterProcBEdgesOwnIter =
448  interProcBEdges[ownProc].begin();
449 
450  // WARNING: Synchronous SLList iterators
451 
452  for
453  (
454  ;
455  curInterProcBdrsOwnIter.good()
456  && curInterProcBEdgesOwnIter.good();
457  ++curInterProcBdrsOwnIter,
458  ++curInterProcBEdgesOwnIter
459  )
460  {
461  if (curInterProcBdrsOwnIter() == neiProc)
462  {
463  // the inter - processor boundary exists. Add the face
464  interProcBouFound = true;
465 
466  bool neighbourFound = false;
467 
468  curInterProcBEdgesOwnIter().append(edgeI);
469 
470  auto curInterProcBdrsNeiIter =
471  interProcBoundaries[neiProc].cbegin();
472 
473  auto curInterProcBEdgesNeiIter =
474  interProcBEdges[neiProc].begin();
475 
476  // WARNING: Synchronous SLList iterators
477 
478  for
479  (
480  ;
481  curInterProcBdrsNeiIter.good()
482  && curInterProcBEdgesNeiIter.good();
483  ++curInterProcBdrsNeiIter,
484  ++curInterProcBEdgesNeiIter
485  )
486  {
487  if (curInterProcBdrsNeiIter() == ownProc)
488  {
489  // boundary found. Add the face
490  neighbourFound = true;
491 
492  curInterProcBEdgesNeiIter().append(edgeI);
493  }
494 
495  if (neighbourFound) break;
496  }
497 
498  if (interProcBouFound && !neighbourFound)
499  {
501  ("faDomainDecomposition::decomposeMesh()")
502  << "Inconsistency in inter - "
503  << "processor boundary lists for processors "
504  << ownProc << " and " << neiProc
505  << abort(FatalError);
506  }
507  }
508 
509  if (interProcBouFound) break;
510  }
511 
512  if (!interProcBouFound)
513  {
514  // inter - processor boundaries do not exist and need to
515  // be created
516 
517  // set the new addressing information
518 
519  // owner
520  interProcBoundaries[ownProc].append(neiProc);
521  interProcBEdges[ownProc].append(SLList<label>(edgeI));
522 
523  // neighbour
524  interProcBoundaries[neiProc].append(ownProc);
525  interProcBEdges[neiProc]
526  .append
527  (
528  SLList<label>(edgeI)
529  );
530  }
531  }
532  }
533 
534 
535  // Loop through patches. For cyclic boundaries detect inter-processor
536  // edges; for all other, add edges to the edge list and remember start
537  // and size of all patches.
538 
539  // for all processors, set the size of start index and patch size
540  // lists to the number of patches in the mesh
541  forAll(procPatchSize_, procI)
542  {
543  procPatchSize_[procI].setSize(patches.size());
544  procPatchStartIndex_[procI].setSize(patches.size());
545  }
546 
547  forAll(patches, patchI)
548  {
549  // Reset size and start index for all processors
550  forAll(procPatchSize_, procI)
551  {
552  procPatchSize_[procI][patchI] = 0;
553  procPatchStartIndex_[procI][patchI] =
554  procEdgeList[procI].size();
555  }
556 
557  const label patchStart = patches[patchI].start();
558 
559 // if (!isA<cyclicFaPatch>(patches[patchI]))
560  if (true)
561  {
562  // Normal patch. Add edges to processor where the face
563  // next to the edge lives
564 
565  const labelListList& eF = patch().edgeFaces();
566 
567  const label size = patches[patchI].labelList::size();
568 
569  labelList patchEdgeFaces(size, -1);
570 
571  for(int eI=0; eI<size; eI++)
572  {
573  patchEdgeFaces[eI] = eF[patches[patchI][eI]][0];
574  }
575 
576  forAll(patchEdgeFaces, edgeI)
577  {
578  const label curProc = faceToProc_[patchEdgeFaces[edgeI]];
579 
580  // add the face
581  procEdgeList[curProc].append(patchStart + edgeI);
582 
583  // increment the number of edges for this patch
584  procPatchSize_[curProc][patchI]++;
585  }
586  }
587  else
588  {
589  // Cyclic patch special treatment
590 
591  const faPatch& cPatch = patches[patchI];
592 
593  const label cycOffset = cPatch.size()/2;
594 
595  // Set reference to faceCells for both patches
596  const labelList::subList firstEdgeFaces
597  (
598  cPatch.edgeFaces(),
599  cycOffset
600  );
601 
602  const labelList::subList secondEdgeFaces
603  (
604  cPatch.edgeFaces(),
605  cycOffset,
606  cycOffset
607  );
608 
609  forAll(firstEdgeFaces, edgeI)
610  {
611  if
612  (
613  faceToProc_[firstEdgeFaces[edgeI]]
614  != faceToProc_[secondEdgeFaces[edgeI]]
615  )
616  {
617  // This edge becomes an inter-processor boundary edge
618  // inter - processor patch edge found. Go through
619  // the list of inside boundaries for the owner
620  // processor and try to find this inter-processor
621  // patch.
622 
623  cyclicParallel_ = true;
624 
625  bool interProcBouFound = false;
626 
627  const label ownProc =
628  faceToProc_[firstEdgeFaces[edgeI]];
629  const label neiProc =
630  faceToProc_[secondEdgeFaces[edgeI]];
631 
632  auto curInterProcBdrsOwnIter =
633  interProcBoundaries[ownProc].cbegin();
634 
635  auto curInterProcBEdgesOwnIter =
636  interProcBEdges[ownProc].begin();
637 
638  // WARNING: Synchronous SLList iterators
639 
640  for
641  (
642  ;
643  curInterProcBdrsOwnIter.good()
644  && curInterProcBEdgesOwnIter.good();
645  ++curInterProcBdrsOwnIter,
646  ++curInterProcBEdgesOwnIter
647  )
648  {
649  if (curInterProcBdrsOwnIter() == neiProc)
650  {
651  // the inter - processor boundary exists.
652  // Add the face
653  interProcBouFound = true;
654 
655  bool neighbourFound = false;
656 
657  curInterProcBEdgesOwnIter()
658  .append(patchStart + edgeI);
659 
660  auto curInterProcBdrsNeiIter
661  = interProcBoundaries[neiProc].cbegin();
662 
663  auto curInterProcBEdgesNeiIter =
664  interProcBEdges[neiProc].begin();
665 
666  // WARNING: Synchronous SLList iterators
667 
668  for
669  (
670  ;
671  curInterProcBdrsNeiIter.good()
672  && curInterProcBEdgesNeiIter.good();
673  ++curInterProcBdrsNeiIter,
674  ++curInterProcBEdgesNeiIter
675  )
676  {
677  if (curInterProcBdrsNeiIter() == ownProc)
678  {
679  // boundary found. Add the face
680  neighbourFound = true;
681 
682  curInterProcBEdgesNeiIter()
683  .append
684  (
685  patchStart
686  + cycOffset
687  + edgeI
688  );
689  }
690 
691  if (neighbourFound) break;
692  }
693 
694  if (interProcBouFound && !neighbourFound)
695  {
697  (
698  "faDomainDecomposition::decomposeMesh()"
699  ) << "Inconsistency in inter-processor "
700  << "boundary lists for processors "
701  << ownProc << " and " << neiProc
702  << " in cyclic boundary matching"
703  << abort(FatalError);
704  }
705  }
706 
707  if (interProcBouFound) break;
708  }
709 
710  if (!interProcBouFound)
711  {
712  // inter - processor boundaries do not exist
713  // and need to be created
714 
715  // set the new addressing information
716 
717  // owner
718  interProcBoundaries[ownProc].append(neiProc);
719  interProcBEdges[ownProc]
720  .append(SLList<label>(patchStart + edgeI));
721 
722  // neighbour
723  interProcBoundaries[neiProc].append(ownProc);
724  interProcBEdges[neiProc]
725  .append
726  (
727  SLList<label>
728  (
729  patchStart
730  + cycOffset
731  + edgeI
732  )
733  );
734  }
735  }
736  else
737  {
738  // This cyclic edge remains on the processor
739  label ownProc = faceToProc_[firstEdgeFaces[edgeI]];
740 
741  // add the edge
742  procEdgeList[ownProc].append(patchStart + edgeI);
743 
744  // increment the number of edges for this patch
745  procPatchSize_[ownProc][patchI]++;
746 
747  // Note: I cannot add the other side of the cyclic
748  // boundary here because this would violate the order.
749  // They will be added in a separate loop below
750  }
751  }
752 
753  // Ordering in cyclic boundaries is important.
754  // Add the other half of cyclic edges for cyclic boundaries
755  // that remain on the processor
756  forAll(secondEdgeFaces, edgeI)
757  {
758  if
759  (
760  faceToProc_[firstEdgeFaces[edgeI]]
761  == faceToProc_[secondEdgeFaces[edgeI]]
762  )
763  {
764  // This cyclic edge remains on the processor
765  label ownProc = faceToProc_[firstEdgeFaces[edgeI]];
766 
767  // add the second edge
768  procEdgeList[ownProc].append
769  (patchStart + cycOffset + edgeI);
770 
771  // increment the number of edges for this patch
772  procPatchSize_[ownProc][patchI]++;
773  }
774  }
775  }
776  }
777 
778  // Convert linked lists into normal lists
779  // Add inter-processor boundaries and remember start indices
780  forAll(procEdgeList, procI)
781  {
782  // Get internal and regular boundary processor faces
783  SLList<label>& curProcEdges = procEdgeList[procI];
784 
785  // Get reference to processor edge addressing
786  labelList& curProcEdgeAddressing = procEdgeAddressing_[procI];
787 
788  labelList& curProcNeighbourProcessors =
789  procNeighbourProcessors_[procI];
790 
791  labelList& curProcProcessorPatchSize =
792  procProcessorPatchSize_[procI];
793 
794  labelList& curProcProcessorPatchStartIndex =
795  procProcessorPatchStartIndex_[procI];
796 
797  // calculate the size
798  label nEdgesOnProcessor = curProcEdges.size();
799 
800  for (const auto& bedges : interProcBEdges[procI])
801  {
802  nEdgesOnProcessor += bedges.size();
803  }
804 
805  curProcEdgeAddressing.setSize(nEdgesOnProcessor);
806 
807  // Fill in the list. Calculate turning index.
808  // Turning index will be -1 only for some edges on processor
809  // boundaries, i.e. the ones where the current processor ID
810  // is in the face which is a edge neighbour.
811  // Turning index is stored as the sign of the edge addressing list
812 
813  label nEdges = 0;
814 
815  // Add internal and boundary edges
816  // Remember to increment the index by one such that the
817  // turning index works properly.
818  for (const label procEdgei : curProcEdges)
819  {
820  curProcEdgeAddressing[nEdges] = procEdgei;
821 // curProcEdgeAddressing[nEdges] = procEdgei + 1;
822  ++nEdges;
823  }
824 
825  // Add inter-processor boundary edges. At the beginning of each
826  // patch, grab the patch start index and size
827 
828  curProcNeighbourProcessors.setSize
829  (
830  interProcBoundaries[procI].size()
831  );
832 
833  curProcProcessorPatchSize.setSize
834  (
835  interProcBoundaries[procI].size()
836  );
837 
838  curProcProcessorPatchStartIndex.setSize
839  (
840  interProcBoundaries[procI].size()
841  );
842 
843  label nProcPatches = 0;
844 
845  auto curInterProcBdrsIter =
846  interProcBoundaries[procI].cbegin();
847 
848  auto curInterProcBEdgesIter =
849  interProcBEdges[procI].cbegin();
850 
851  for
852  (
853  ;
854  curInterProcBdrsIter.good()
855  && curInterProcBEdgesIter.good();
856  ++curInterProcBdrsIter,
857  ++curInterProcBEdgesIter
858  )
859  {
860  curProcNeighbourProcessors[nProcPatches] =
861  curInterProcBdrsIter();
862 
863  // Get start index for processor patch
864  curProcProcessorPatchStartIndex[nProcPatches] = nEdges;
865 
866  label& curSize =
867  curProcProcessorPatchSize[nProcPatches];
868 
869  curSize = 0;
870 
871  // add faces for this processor boundary
872 
873  for (const label edgei : *curInterProcBEdgesIter)
874  {
875  // add the edges
876 
877  // Remember to increment the index by one such that the
878  // turning index works properly.
879  if (faceToProc_[owner[edgei]] == procI)
880  {
881  curProcEdgeAddressing[nEdges] = edgei;
882 // curProcEdgeAddressing[nEdges] = edgei + 1;
883  }
884  else
885  {
886  // turning edge
887  curProcEdgeAddressing[nEdges] = edgei;
888 // curProcEdgeAddressing[nEdges] = -(edgei + 1);
889  }
890 
891  // increment the size
892  ++curSize;
893 
894  ++nEdges;
895  }
896 
897  ++nProcPatches;
898  }
899  }
900  }
901 
902  Info << "\nCalculating processor boundary addressing" << endl;
903  // For every patch of processor boundary, find the index of the original
904  // patch. Mis-alignment is caused by the fact that patches with zero size
905  // are omitted. For processor patches, set index to -1.
906  // At the same time, filter the procPatchSize_ and procPatchStartIndex_
907  // lists to exclude zero-size patches
908  forAll(procPatchSize_, procI)
909  {
910  // Make a local copy of old lists
911  const labelList oldPatchSizes = procPatchSize_[procI];
912 
913  const labelList oldPatchStarts = procPatchStartIndex_[procI];
914 
915  labelList& curPatchSizes = procPatchSize_[procI];
916 
917  labelList& curPatchStarts = procPatchStartIndex_[procI];
918 
919  const labelList& curProcessorPatchSizes =
920  procProcessorPatchSize_[procI];
921 
922  labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
923 
924  curBoundaryAddressing.setSize
925  (
926  oldPatchSizes.size()
927  + curProcessorPatchSizes.size()
928  );
929 
930  label nPatches = 0;
931 
932  forAll(oldPatchSizes, patchI)
933  {
934  //- Do not suppress zero sized patches since make parallel
935  // actions inside patches near impossible.
936  //if (oldPatchSizes[patchI] > 0)
937  {
938  curBoundaryAddressing[nPatches] = patchI;
939 
940  curPatchSizes[nPatches] = oldPatchSizes[patchI];
941 
942  curPatchStarts[nPatches] = oldPatchStarts[patchI];
943 
944  nPatches++;
945  }
946  }
947 
948  // reset to the size of live patches
949  curPatchSizes.setSize(nPatches);
950  curPatchStarts.setSize(nPatches);
951 
952  forAll(curProcessorPatchSizes, procPatchI)
953  {
954  curBoundaryAddressing[nPatches] = -1;
955 
956  nPatches++;
957  }
958 
959  curBoundaryAddressing.setSize(nPatches);
960  }
961 
962 
963  // Gather data about globally shared points
964 
965  labelList globallySharedPoints_(0);
966 
967  // Memory management
968  {
969  labelList pointsUsage(nPoints(), Zero);
970 
971  // Globally shared points are the ones used by more than 2 processors
972  // Size the list approximately and gather the points
973  labelHashSet gSharedPoints
974  (
975  min(100, nPoints()/1000)
976  );
977 
978  // Loop through all the processors and mark up points used by
979  // processor boundaries. When a point is used twice, it is a
980  // globally shared point
981 
982  for (label procI = 0; procI < nProcs(); procI++)
983  {
984  // Get list of edge labels
985  const labelList& curEdgeLabels = procEdgeAddressing_[procI];
986 
987  // Get start of processor faces
988  const labelList& curProcessorPatchStarts =
989  procProcessorPatchStartIndex_[procI];
990 
991  const labelList& curProcessorPatchSizes =
992  procProcessorPatchSize_[procI];
993 
994  // Reset the lookup list
995  pointsUsage = 0;
996 
997  forAll(curProcessorPatchStarts, patchI)
998  {
999  const label curStart = curProcessorPatchStarts[patchI];
1000  const label curEnd = curStart + curProcessorPatchSizes[patchI];
1001 
1002  for
1003  (
1004  label edgeI = curStart;
1005  edgeI < curEnd;
1006  edgeI++
1007  )
1008  {
1009  // Mark the original edge as used
1010  // Remember to decrement the index by one (turning index)
1011  const label curE = curEdgeLabels[edgeI];
1012 
1013  const edge& e = edges[curE];
1014 
1015  forAll(e, pointI)
1016  {
1017  if (pointsUsage[e[pointI]] == 0)
1018  {
1019  // Point not previously used
1020  pointsUsage[e[pointI]] = patchI + 1;
1021  }
1022  else if (pointsUsage[e[pointI]] != patchI + 1)
1023  {
1024  // Point used by some other patch = global point!
1025  gSharedPoints.insert(e[pointI]);
1026  }
1027  }
1028  }
1029  }
1030  }
1031 
1032  // Grab the result from the hash list
1033  globallySharedPoints_ = gSharedPoints.sortedToc();
1034  }
1035 
1036 
1037  // Edge label for faPatches
1038 
1039  for (label procI = 0; procI < nProcs(); procI++)
1040  {
1041  fileName processorCasePath
1042  (
1043  time().caseName()/("processor" + Foam::name(procI))
1044  );
1045 
1046  // create a database
1047  Time processorDb
1048  (
1050  time().rootPath(),
1051  processorCasePath
1052  );
1053 
1054 
1055  // read finite volume mesh
1056  fvMesh procFvMesh
1057  (
1058  IOobject
1059  (
1061  processorDb.timeName(),
1062  processorDb
1063  )
1064  );
1065 
1066  // create finite area mesh
1067  faMesh procMesh
1068  (
1069  procFvMesh,
1070  procFaceLabels_[procI]
1071  );
1072 
1073 
1074  const labelList& curEdgeAddressing = procEdgeAddressing_[procI];
1075 
1076  const labelList& curPatchStartIndex = procPatchStartIndex_[procI];
1077  const labelList& curPatchSize = procPatchSize_[procI];
1078 
1079  const labelList& curProcessorPatchStartIndex =
1080  procProcessorPatchStartIndex_[procI];
1081 
1082  const labelList& curProcessorPatchSize =
1083  procProcessorPatchSize_[procI];
1084 
1085  labelListList& curPatchEdgeLabels = procPatchEdgeLabels_[procI];
1086  curPatchEdgeLabels =
1088  (
1089  curPatchSize.size()
1090  + curProcessorPatchSize.size()
1091  );
1092 
1093  forAll(curPatchSize, patchI)
1094  {
1095  labelList& curEdgeLabels = curPatchEdgeLabels[patchI];
1096  curEdgeLabels.setSize(curPatchSize[patchI], -1);
1097 
1098  label edgeI = 0;
1099 
1100  for
1101  (
1102  int i=curPatchStartIndex[patchI];
1103  i<(curPatchStartIndex[patchI]+curPatchSize[patchI]);
1104  i++
1105  )
1106  {
1107  curEdgeLabels[edgeI] =
1108  procMeshEdgesMap_[procI][curEdgeAddressing[i]];
1109  edgeI++;
1110  }
1111  }
1112 
1113  forAll(curProcessorPatchSize, patchI)
1114  {
1115  labelList& curEdgeLabels =
1116  curPatchEdgeLabels[curPatchSize.size() + patchI];
1117  curEdgeLabels.setSize(curProcessorPatchSize[patchI], -1);
1118 
1119  label edgeI = 0;
1120 
1121  for
1122  (
1123  int i=curProcessorPatchStartIndex[patchI];
1124  i<(curProcessorPatchStartIndex[patchI]
1125  +curProcessorPatchSize[patchI]);
1126  i++
1127  )
1128  {
1129  curEdgeLabels[edgeI] =
1130  procMeshEdgesMap_[procI][curEdgeAddressing[i]];
1131  edgeI++;
1132  }
1133  }
1134  }
1135 }
1136 
1137 
1139 {
1140  Info<< "\nConstructing processor FA meshes" << endl;
1141 
1142 
1143  // Make a lookup map for globally shared points
1144  Map<label> sharedPointLookup(2*globallySharedPoints_.size());
1145 
1146  forAll(globallySharedPoints_, pointi)
1147  {
1148  sharedPointLookup.insert(globallySharedPoints_[pointi], pointi);
1149  }
1150 
1151  label totProcEdges = 0;
1152  label maxProcPatches = 0;
1153  label maxProcEdges = 0;
1154 
1155  // Write out the meshes
1156  for (label procI = 0; procI < nProcs(); procI++)
1157  {
1158  // Create processor mesh without a boundary
1159 
1160  fileName processorCasePath
1161  (
1162  time().caseName()/("processor" + Foam::name(procI))
1163  );
1164 
1165  // create a database
1166  Time processorDb
1167  (
1169  time().rootPath(),
1170  processorCasePath
1171  );
1172 
1173  // read finite volume mesh
1174  fvMesh procFvMesh
1175  (
1176  IOobject
1177  (
1179  processorDb.timeName(),
1180  processorDb
1181  )
1182  );
1183 
1184  labelIOList fvBoundaryProcAddressing
1185  (
1186  IOobject
1187  (
1188  "boundaryProcAddressing",
1189  "constant",
1190  procFvMesh.meshSubDir,
1191  procFvMesh,
1194  )
1195  );
1196 
1197 
1198  // create finite area mesh
1199  faMesh procMesh
1200  (
1201  procFvMesh,
1202  procFaceLabels_[procI]
1203  );
1204 
1205  // Create processor boundary patches
1206  const labelList& curBoundaryAddressing =
1207  procBoundaryAddressing_[procI];
1208 
1209  const labelList& curPatchSizes = procPatchSize_[procI];
1210 
1211  const labelList& curNeighbourProcessors =
1212  procNeighbourProcessors_[procI];
1213 
1214  const labelList& curProcessorPatchSizes =
1215  procProcessorPatchSize_[procI];
1216 
1217  const labelListList& curPatchEdgeLabels =
1218  procPatchEdgeLabels_[procI];
1219 
1220  const faPatchList& meshPatches = boundary();
1221 
1222  List<faPatch*> procPatches
1223  (
1224  curPatchSizes.size()
1225  + curProcessorPatchSizes.size(),
1226  reinterpret_cast<faPatch*>(NULL)
1227  );
1228 
1229  label nPatches = 0;
1230 
1231  forAll(curPatchSizes, patchi)
1232  {
1233  const labelList& curEdgeLabels = curPatchEdgeLabels[nPatches];
1234 
1235  const label ngbPolyPatchIndex =
1236  fvBoundaryProcAddressing.find
1237  (
1238  meshPatches[curBoundaryAddressing[patchi]]
1239  .ngbPolyPatchIndex()
1240  );
1241 
1242  procPatches[nPatches] =
1243  meshPatches[curBoundaryAddressing[patchi]].clone
1244  (
1245  procMesh.boundary(),
1246  curEdgeLabels,
1247  nPatches,
1248  ngbPolyPatchIndex
1249  ).ptr();
1250 
1251  nPatches++;
1252  }
1253 
1254  forAll(curProcessorPatchSizes, procPatchI)
1255  {
1256  const labelList& curEdgeLabels = curPatchEdgeLabels[nPatches];
1257 
1258  procPatches[nPatches] =
1259  new processorFaPatch
1260  (
1261  word("procBoundary") + Foam::name(procI)
1262  + word("to")
1263  + Foam::name(curNeighbourProcessors[procPatchI]),
1264  curEdgeLabels,
1265  nPatches,
1266  procMesh.boundary(),
1267  -1,
1268  procI,
1269  curNeighbourProcessors[procPatchI]
1270  );
1271 
1272  nPatches++;
1273  }
1274 
1275  // Add boundary patches
1276  procMesh.addFaPatches(procPatches);
1277 
1278  // Set the precision of the points data to 10
1280 
1281  procMesh.write();
1282 
1283  Info<< endl
1284  << "Processor " << procI << nl
1285  << " Number of faces = " << procMesh.nFaces()
1286  << endl;
1287 
1288  label nBoundaryEdges = 0;
1289  label nProcPatches = 0;
1290  label nProcEdges = 0;
1291 
1292  forAll(procMesh.boundary(), patchi)
1293  {
1294  if
1295  (
1296  procMesh.boundary()[patchi].type()
1297  == processorFaPatch::typeName
1298  )
1299  {
1300  const processorFaPatch& ppp =
1301  refCast<const processorFaPatch>
1302  (
1303  procMesh.boundary()[patchi]
1304  );
1305 
1306  Info<< " Number of edges shared with processor "
1307  << ppp.neighbProcNo() << " = " << ppp.size() << endl;
1308 
1309  nProcPatches++;
1310  nProcEdges += ppp.size();
1311  }
1312  else
1313  {
1314  nBoundaryEdges += procMesh.boundary()[patchi].size();
1315  }
1316  }
1317 
1318  Info<< " Number of processor patches = " << nProcPatches << nl
1319  << " Number of processor edges = " << nProcEdges << nl
1320  << " Number of boundary edges = " << nBoundaryEdges << endl;
1321 
1322  totProcEdges += nProcEdges;
1323  maxProcPatches = max(maxProcPatches, nProcPatches);
1324  maxProcEdges = max(maxProcEdges, nProcEdges);
1325 
1326  // create and write the addressing information
1327  labelIOList pointProcAddressing
1328  (
1329  IOobject
1330  (
1331  "pointProcAddressing",
1332  "constant",
1333  procMesh.meshSubDir,
1334  procFvMesh,
1337  ),
1338  procPatchPointAddressing_[procI]
1339  );
1340  pointProcAddressing.write();
1341 
1342  labelIOList edgeProcAddressing
1343  (
1344  IOobject
1345  (
1346  "edgeProcAddressing",
1347  "constant",
1348  procMesh.meshSubDir,
1349  procFvMesh,
1352  ),
1353  procEdgeAddressing_[procI]
1354  );
1355  edgeProcAddressing.write();
1356 
1358  (
1359  IOobject
1360  (
1361  "faceProcAddressing",
1362  "constant",
1363  procMesh.meshSubDir,
1364  procFvMesh,
1367  ),
1368  procFaceAddressing_[procI]
1369  );
1370  faceProcAddressing.write();
1371 
1372  labelIOList boundaryProcAddressing
1373  (
1374  IOobject
1375  (
1376  "boundaryProcAddressing",
1377  "constant",
1378  procMesh.meshSubDir,
1379  procFvMesh,
1382  ),
1383  procBoundaryAddressing_[procI]
1384  );
1385  boundaryProcAddressing.write();
1386  }
1387 
1388  Info<< nl
1389  << "Number of processor edges = " << totProcEdges/2 << nl
1390  << "Max number of processor patches = " << maxProcPatches << nl
1391  << "Max number of faces between processors = " << maxProcEdges
1392  << endl;
1393 
1394  return true;
1395 }
1396 
1397 
1398 // ************************************************************************* //
Foam::IOobject::NO_WRITE
Definition: IOobject.H:130
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:71
Foam::faMesh::name
const word & name() const
Name function is needed to disambiguate those inherited.
Definition: faMesh.H:415
faceProcAddressing
PtrList< labelIOList > & faceProcAddressing
Definition: checkFaceAddressingComp.H:9
OSspecific.H
Functions used by OpenFOAM that are specific to POSIX compliant operating systems and need to be repl...
Foam::edgeList
List< edge > edgeList
A List of edges.
Definition: edgeList.H:63
nPatches
label nPatches
Definition: readKivaGrid.H:396
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
globalMeshData.H
Foam::IOobject::rootPath
const fileName & rootPath() const
Definition: IOobject.C:469
processorFaPatch.H
Foam::faMeshDecomposition::decomposeMesh
void decomposeMesh()
Decompose mesh.
Foam::faMesh::faceLabels
const labelList & faceLabels() const
Return faMesh face labels.
Definition: faMesh.H:424
Foam::IOobject::IOobject
IOobject(const IOobject &)=default
Copy construct.
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:350
faMesh.H
Foam::GeoMesh< polyMesh >::mesh_
const polyMesh & mesh_
Reference to Mesh.
Definition: GeoMesh.H:56
Foam::labelIOList
IOList< label > labelIOList
Label container classes.
Definition: labelIOList.H:44
Foam::Time::controlDictName
static word controlDictName
The default control dictionary name (normally "controlDict")
Definition: Time.H:226
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
nProcPatches
const label nProcPatches
Definition: convertProcessorPatches.H:171
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
nPoints
label nPoints
Definition: gmvOutputHeader.H:2
Foam::faMeshDecomposition::faMeshDecomposition
faMeshDecomposition(const fvMesh &mesh, const fileName &decompDictFile="")
Construct from components.
Foam::faMeshDecomposition::nProcs
label nProcs() const
Number of processor in decomposition.
Definition: faMeshDecomposition.H:155
Map.H
Foam::cpuTime
cpuTimeCxx cpuTime
Selection of preferred clock mechanism for the elapsed cpu time.
Definition: cpuTimeFwd.H:41
faMeshDecomposition.H
Foam::Info
messageStream Info
Information stream (uses stdout - output is on the master only)
Foam::name
word name(const complex &c)
Return string representation of complex.
Definition: complex.C:76
Foam::faMesh::time
const Time & time() const
Return reference to time.
Definition: faMesh.C:897
Foam::IOobject::caseName
const fileName & caseName() const
Definition: IOobject.C:475
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::List< label >::subList
SubList< label > subList
Declare type of subList.
Definition: List.H:114
Foam::FatalError
error FatalError
mesh
dynamicFvMesh & mesh
Definition: createDynamicFvMesh.H:6
Foam::abort
errorManip< error > abort(error &err)
Definition: errorManip.H:144
Foam::indirectPrimitivePatch
PrimitivePatch< IndirectList< face >, const pointField & > indirectPrimitivePatch
A PrimitivePatch with an IndirectList for the faces, const reference for the point field.
Definition: indirectPrimitivePatch.H:49
Foam::New
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
Global function forwards to reuseTmpDimensionedField::New.
Definition: DimensionedFieldReuseFunctions.H:105
Foam::decompositionModel::New
static const decompositionModel & New(const polyMesh &mesh, const fileName &decompDictFile="")
Read (optionally from absolute path) and register on mesh.
Definition: decompositionModel.C:116
Time.H
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
Foam::nl
constexpr char nl
Definition: Ostream.H:385
Foam::IOstream::defaultPrecision
static unsigned int defaultPrecision()
Return the default precision.
Definition: IOstream.H:333
Foam::foamVersion::patch
const std::string patch
OpenFOAM patch number as a std::string.
labelIOList.H
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
dictionary.H
patches
const polyBoundaryMesh & patches
Definition: convertProcessorPatches.H:65
FatalErrorIn
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:376
decompositionModel.H
Foam::faPatchList
PtrList< faPatch > faPatchList
Definition: faPatchList.H:47
Foam::List::clone
autoPtr< List< T > > clone() const
Clone.
Definition: ListI.H:99
Foam::labelHashSet
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys and label hasher.
Definition: HashSet.H:409
Foam::List::setSize
void setSize(const label newSize)
Alias for resize(const label)
Definition: ListI.H:146
Foam::IOobject::NO_READ
Definition: IOobject.H:123
Foam::faMeshDecomposition::writeDecomposition
bool writeDecomposition()
Write decomposition.
boundary
faceListList boundary
Definition: createBlockMesh.H:4
Foam::IOobject::MUST_READ
Definition: IOobject.H:120