All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
domainDecompositionMesh.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2020 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 InClass
25  domainDecomposition
26 
27 Description
28  Private member of domainDecomposition.
29  Decomposes the mesh into bits
30 
31 \*---------------------------------------------------------------------------*/
32 
33 #include "domainDecomposition.H"
34 #include "IOstreams.H"
35 #include "boolList.H"
36 #include "cyclicPolyPatch.H"
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 void Foam::domainDecomposition::append(labelList& lst, const label elem)
41 {
42  label sz = lst.size();
43  lst.setSize(sz+1);
44  lst[sz] = elem;
45 }
46 
47 
48 void Foam::domainDecomposition::addInterProcFace
49 (
50  const label facei,
51  const label ownerProc,
52  const label nbrProc,
53 
54  List<Map<label>>& nbrToInterPatch,
55  List<DynamicList<DynamicList<label>>>& interPatchFaces
56 ) const
57 {
58  Map<label>::iterator patchiter = nbrToInterPatch[ownerProc].find(nbrProc);
59 
60  // Introduce turning index only for internal faces (are duplicated).
61  label ownerIndex = facei+1;
62  label nbrIndex = -(facei+1);
63 
64  if (patchiter != nbrToInterPatch[ownerProc].end())
65  {
66  // Existing interproc patch. Add to both sides.
67  label toNbrProcPatchi = patchiter();
68  interPatchFaces[ownerProc][toNbrProcPatchi].append(ownerIndex);
69 
70  if (isInternalFace(facei))
71  {
72  label toOwnerProcPatchi = nbrToInterPatch[nbrProc][ownerProc];
73  interPatchFaces[nbrProc][toOwnerProcPatchi].append(nbrIndex);
74  }
75  }
76  else
77  {
78  // Create new interproc patches.
79  label toNbrProcPatchi = nbrToInterPatch[ownerProc].size();
80  nbrToInterPatch[ownerProc].insert(nbrProc, toNbrProcPatchi);
81  DynamicList<label> oneFace;
82  oneFace.append(ownerIndex);
83  interPatchFaces[ownerProc].append(oneFace);
84 
85  if (isInternalFace(facei))
86  {
87  label toOwnerProcPatchi = nbrToInterPatch[nbrProc].size();
88  nbrToInterPatch[nbrProc].insert(ownerProc, toOwnerProcPatchi);
89  oneFace.clear();
90  oneFace.append(nbrIndex);
91  interPatchFaces[nbrProc].append(oneFace);
92  }
93  }
94 }
95 
96 
97 void Foam::domainDecomposition::decomposeMesh(const fileName& dict)
98 {
99  // Decide which cell goes to which processor
100  distributeCells(dict);
101 
102  // Distribute the cells according to the given processor label
103 
104  // calculate the addressing information for the original mesh
105  Info<< "\nCalculating original mesh data" << endl;
106 
107  // set references to the original mesh
108  const polyBoundaryMesh& patches = boundaryMesh();
109  const faceList& fcs = faces();
110  const labelList& owner = faceOwner();
111  const labelList& neighbour = faceNeighbour();
112 
113  // loop through the list of processor labels for the cell and add the
114  // cell shape to the list of cells for the appropriate processor
115 
116  Info<< "\nDistributing cells to processors" << endl;
117 
118  // Cells per processor
119  procCellAddressing_ = invertOneToMany(nProcs_, cellToProc_);
120 
121  Info<< "\nDistributing faces to processors" << endl;
122 
123  // Loop through all internal faces and decide which processor they belong to
124  // First visit all internal faces. If cells at both sides belong to the
125  // same processor, the face is an internal face. If they are different,
126  // it belongs to both processors.
127 
128  procFaceAddressing_.setSize(nProcs_);
129 
130  // Internal faces
131  forAll(neighbour, facei)
132  {
133  if (cellToProc_[owner[facei]] == cellToProc_[neighbour[facei]])
134  {
135  // Face internal to processor. Notice no turning index.
136  procFaceAddressing_[cellToProc_[owner[facei]]].append(facei+1);
137  }
138  }
139 
140  // for all processors, set the size of start index and patch size
141  // lists to the number of patches in the mesh
142  forAll(procPatchSize_, proci)
143  {
144  procPatchSize_[proci].setSize(patches.size());
145  procPatchStartIndex_[proci].setSize(patches.size());
146  }
147 
148  forAll(patches, patchi)
149  {
150  // Reset size and start index for all processors
151  forAll(procPatchSize_, proci)
152  {
153  procPatchSize_[proci][patchi] = 0;
154  procPatchStartIndex_[proci][patchi] =
155  procFaceAddressing_[proci].size();
156  }
157 
158  const label patchStart = patches[patchi].start();
159 
160  if (!isA<cyclicPolyPatch>(patches[patchi]))
161  {
162  // Normal patch. Add faces to processor where the cell
163  // next to the face lives
164 
165  const labelUList& patchFaceCells =
166  patches[patchi].faceCells();
167 
168  forAll(patchFaceCells, facei)
169  {
170  const label curProc = cellToProc_[patchFaceCells[facei]];
171 
172  // add the face without turning index
173  procFaceAddressing_[curProc].append(patchStart+facei+1);
174 
175  // increment the number of faces for this patch
176  procPatchSize_[curProc][patchi]++;
177  }
178  }
179  else
180  {
181  const cyclicPolyPatch& pp = refCast<const cyclicPolyPatch>
182  (
183  patches[patchi]
184  );
185  // cyclic: check opposite side on this processor
186  const labelUList& patchFaceCells = pp.faceCells();
187 
188  const labelUList& nbrPatchFaceCells =
189  pp.nbrPatch().faceCells();
190 
191  forAll(patchFaceCells, facei)
192  {
193  const label curProc = cellToProc_[patchFaceCells[facei]];
194  const label nbrProc = cellToProc_[nbrPatchFaceCells[facei]];
195  if (curProc == nbrProc)
196  {
197  // add the face without turning index
198  procFaceAddressing_[curProc].append(patchStart+facei+1);
199  // increment the number of faces for this patch
200  procPatchSize_[curProc][patchi]++;
201  }
202  }
203  }
204  }
205 
206 
207  // Done internal bits of the new mesh and the ordinary patches.
208 
209 
210  // Per processor, from neighbour processor to the inter-processor patch
211  // that communicates with that neighbour
212  List<Map<label>> procNbrToInterPatch(nProcs_);
213 
214  // Per processor the faces per inter-processor patch
215  List<DynamicList<DynamicList<label>>> interPatchFaces(nProcs_);
216 
217  // Processor boundaries from internal faces
218  forAll(neighbour, facei)
219  {
220  label ownerProc = cellToProc_[owner[facei]];
221  label nbrProc = cellToProc_[neighbour[facei]];
222 
223  if (ownerProc != nbrProc)
224  {
225  // inter - processor patch face found.
226  addInterProcFace
227  (
228  facei,
229  ownerProc,
230  nbrProc,
231 
232  procNbrToInterPatch,
233  interPatchFaces
234  );
235  }
236  }
237 
238  // Add the proper processor faces to the sub information. For faces
239  // originating from internal faces this is always -1.
240  List<labelListList> subPatchIDs(nProcs_);
241  List<labelListList> subPatchStarts(nProcs_);
242  forAll(interPatchFaces, proci)
243  {
244  label nInterfaces = interPatchFaces[proci].size();
245 
246  subPatchIDs[proci].setSize(nInterfaces, labelList(1, label(-1)));
247  subPatchStarts[proci].setSize(nInterfaces, labelList(1, label(0)));
248  }
249 
250 
251  // Special handling needed for the case that multiple processor cyclic
252  // patches are created on each local processor domain, e.g. if a 3x3 case
253  // is decomposed using the decomposition:
254  //
255  // | 1 | 0 | 2 |
256  // cyclic left | 2 | 0 | 1 | cyclic right
257  // | 2 | 0 | 1 |
258  //
259  // - processors 1 and 2 will both have pieces of both cyclic left- and
260  // right sub-patches present
261  // - the interface patch faces are stored in a single list, where each
262  // sub-patch is referenced into the list using a patch start index and
263  // size
264  // - if the patches are in order (in the boundary file) of left, right
265  // - processor 1 will send: left, right
266  // - processor 1 will need to receive in reverse order: right, left
267  // - similarly for processor 2
268  // - the sub-patches are therefore generated in 4 passes of the patch lists
269  // 1. add faces from owner patch where local proc i < nbr proc i
270  // 2. add faces from nbr patch where local proc i < nbr proc i
271  // 3. add faces from owner patch where local proc i > nbr proc i
272  // 4. add faces from nbr patch where local proc i > nbr proc i
273 
274  processInterCyclics
275  (
276  patches,
277  interPatchFaces,
278  procNbrToInterPatch,
279  subPatchIDs,
280  subPatchStarts,
281  true,
282  lessOp<label>()
283  );
284 
285  processInterCyclics
286  (
287  patches,
288  interPatchFaces,
289  procNbrToInterPatch,
290  subPatchIDs,
291  subPatchStarts,
292  false,
293  lessOp<label>()
294  );
295 
296  processInterCyclics
297  (
298  patches,
299  interPatchFaces,
300  procNbrToInterPatch,
301  subPatchIDs,
302  subPatchStarts,
303  false,
304  greaterOp<label>()
305  );
306 
307  processInterCyclics
308  (
309  patches,
310  interPatchFaces,
311  procNbrToInterPatch,
312  subPatchIDs,
313  subPatchStarts,
314  true,
315  greaterOp<label>()
316  );
317 
318 
319  // Sort inter-proc patch by neighbour
320  labelList order;
321  forAll(procNbrToInterPatch, proci)
322  {
323  label nInterfaces = procNbrToInterPatch[proci].size();
324 
325  procNeighbourProcessors_[proci].setSize(nInterfaces);
326  procProcessorPatchSize_[proci].setSize(nInterfaces);
327  procProcessorPatchStartIndex_[proci].setSize(nInterfaces);
328  procProcessorPatchSubPatchIDs_[proci].setSize(nInterfaces);
329  procProcessorPatchSubPatchStarts_[proci].setSize(nInterfaces);
330 
331  // Info<< "Processor " << proci << endl;
332 
333  // Get sorted neighbour processors
334  const Map<label>& curNbrToInterPatch = procNbrToInterPatch[proci];
335  labelList nbrs = curNbrToInterPatch.toc();
336 
337  sortedOrder(nbrs, order);
338 
339  DynamicList<DynamicList<label>>& curInterPatchFaces =
340  interPatchFaces[proci];
341 
342  forAll(nbrs, i)
343  {
344  const label nbrProc = nbrs[i];
345  const label interPatch = curNbrToInterPatch[nbrProc];
346 
347  procNeighbourProcessors_[proci][i] = nbrProc;
348  procProcessorPatchSize_[proci][i] =
349  curInterPatchFaces[interPatch].size();
350  procProcessorPatchStartIndex_[proci][i] =
351  procFaceAddressing_[proci].size();
352 
353  // Add size as last element to substarts and transfer
354  append
355  (
356  subPatchStarts[proci][interPatch],
357  curInterPatchFaces[interPatch].size()
358  );
359  procProcessorPatchSubPatchIDs_[proci][i].transfer
360  (
361  subPatchIDs[proci][interPatch]
362  );
363  procProcessorPatchSubPatchStarts_[proci][i].transfer
364  (
365  subPatchStarts[proci][interPatch]
366  );
367 
368  // Info<< " nbr:" << nbrProc << endl;
369  // Info<< " interpatch:" << interPatch << endl;
370  // Info<< " size:" << procProcessorPatchSize_[proci][i] << endl;
371  // Info<< " start:" << procProcessorPatchStartIndex_[proci][i]
372  // << endl;
373  // Info<< " subPatches:"
374  // << procProcessorPatchSubPatchIDs_[proci][i]
375  // << endl;
376  // Info<< " subStarts:"
377  // << procProcessorPatchSubPatchStarts_[proci][i] << endl;
378 
379  // And add all the face labels for interPatch
380  DynamicList<label>& interPatchFaces =
381  curInterPatchFaces[interPatch];
382 
383  forAll(interPatchFaces, j)
384  {
385  procFaceAddressing_[proci].append(interPatchFaces[j]);
386  }
387  interPatchFaces.clearStorage();
388  }
389  curInterPatchFaces.clearStorage();
390  procFaceAddressing_[proci].shrink();
391  }
392 
393 
396 // forAll(procPatchStartIndex_, proci)
397 // {
398 // Info<< "Processor:" << proci << endl;
399 //
400 // Info<< " total faces:" << procFaceAddressing_[proci].size()
401 // << endl;
402 //
403 // const labelList& curProcPatchStartIndex = procPatchStartIndex_[proci];
404 //
405 // forAll(curProcPatchStartIndex, patchi)
406 // {
407 // Info<< " patch:" << patchi
408 // << "\tstart:" << curProcPatchStartIndex[patchi]
409 // << "\tsize:" << procPatchSize_[proci][patchi]
410 // << endl;
411 // }
412 // }
413 // Info<< endl;
414 //
415 // forAll(procNeighbourProcessors_, proci)
416 // {
417 // Info<< "Processor " << proci << endl;
418 //
419 // forAll(procNeighbourProcessors_[proci], i)
420 // {
421 // Info<< " nbr:" << procNeighbourProcessors_[proci][i] << endl;
422 // Info<< " size:" << procProcessorPatchSize_[proci][i] << endl;
423 // Info<< " start:" << procProcessorPatchStartIndex_[proci][i]
424 // << endl;
425 // }
426 // }
427 // Info<< endl;
428 //
429 // forAll(procFaceAddressing_, proci)
430 // {
431 // Info<< "Processor:" << proci << endl;
432 //
433 // Info<< " faces:" << procFaceAddressing_[proci] << endl;
434 // }
435 
436 
437 
438  Info<< "\nDistributing points to processors" << endl;
439  // For every processor, loop through the list of faces for the processor.
440  // For every face, loop through the list of points and mark the point as
441  // used for the processor. Collect the list of used points for the
442  // processor.
443 
444  forAll(procPointAddressing_, proci)
445  {
446  boolList pointLabels(nPoints(), false);
447 
448  // Get reference to list of used faces
449  const labelList& procFaceLabels = procFaceAddressing_[proci];
450 
451  forAll(procFaceLabels, facei)
452  {
453  // Because of the turning index, some labels may be negative
454  const labelList& facePoints = fcs[mag(procFaceLabels[facei]) - 1];
455 
456  forAll(facePoints, pointi)
457  {
458  // Mark the point as used
459  pointLabels[facePoints[pointi]] = true;
460  }
461  }
462 
463  // Collect the used points
464  labelList& procPointLabels = procPointAddressing_[proci];
465 
466  procPointLabels.setSize(pointLabels.size());
467 
468  label nUsedPoints = 0;
469 
470  forAll(pointLabels, pointi)
471  {
472  if (pointLabels[pointi])
473  {
474  procPointLabels[nUsedPoints] = pointi;
475 
476  nUsedPoints++;
477  }
478  }
479 
480  // Reset the size of used points
481  procPointLabels.setSize(nUsedPoints);
482  }
483 }
484 
485 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
void sortedOrder(const UList< T > &, labelList &order)
Generate the (stable) sort order for the list.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:112
labelList pointLabels(nPoints, -1)
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1175
HashTable< label, label, Hash< label > >::iterator iterator
Definition: Map.H:56
labelListList invertOneToMany(const label len, const labelUList &)
Invert one-to-many map. Unmapped elements will be size 0.
Definition: ListOps.C:67
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
List< face > faceList
Definition: faceListFwd.H:43
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
label size() const
Return number of elements in table.
Definition: HashTableI.H:65
UList< label > labelUList
Definition: UList.H:65
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
void decomposeMesh(const fileName &dict)
Decompose mesh.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
List< label > labelList
A List of labels.
Definition: labelList.H:56
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1156
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342