Merge branch 'master' of ssh://git.code.sf.net/p/foam-extend/foam-extend-3.2
[foam-extend-3.2.git] / applications / utilities / parallelProcessing / decomposePar / domainDecomposition.C
blobf2153aef38b5fc70339dd305f0b985534fe3b746
1 /*---------------------------------------------------------------------------*\
2   =========                 |
3   \\      /  F ield         | foam-extend: Open Source CFD
4    \\    /   O peration     | Version:     3.2
5     \\  /    A nd           | Web:         http://www.foam-extend.org
6      \\/     M anipulation  | For copyright notice see file Copyright
7 -------------------------------------------------------------------------------
8 License
9     This file is part of foam-extend.
11     foam-extend is free software: you can redistribute it and/or modify it
12     under the terms of the GNU General Public License as published by the
13     Free Software Foundation, either version 3 of the License, or (at your
14     option) any later version.
16     foam-extend is distributed in the hope that it will be useful, but
17     WITHOUT ANY WARRANTY; without even the implied warranty of
18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
19     General Public License for more details.
21     You should have received a copy of the GNU General Public License
22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
24 \*---------------------------------------------------------------------------*/
26 #include "domainDecomposition.H"
27 #include "foamTime.H"
28 #include "dictionary.H"
29 #include "labelIOList.H"
30 #include "processorPolyPatch.H"
31 #include "fvMesh.H"
32 #include "OSspecific.H"
33 #include "Map.H"
34 #include "globalMeshData.H"
35 #include "DynamicList.H"
37 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
39 domainDecomposition::domainDecomposition(const IOobject& io)
41     fvMesh(io),
42     decompositionDict_
43     (
44         IOobject
45         (
46             "decomposeParDict",
47             time().system(),
48             *this,
49             IOobject::MUST_READ,
50             IOobject::NO_WRITE
51         )
52     ),
53     nProcs_(readInt(decompositionDict_.lookup("numberOfSubdomains"))),
54     distributed_(false),
55     cellToProc_(nCells()),
56     procPointAddressing_(nProcs_),
57     procFaceAddressing_(nProcs_),
58     nInternalProcFaces_(nProcs_),
59     nLiveProcFaces_(nProcs_),
60     procCellAddressing_(nProcs_),
61     procBoundaryAddressing_(nProcs_),
62     procPatchSize_(nProcs_),
63     procPatchStartIndex_(nProcs_),
64     procNeighbourProcessors_(nProcs_),
65     procProcessorPatchSize_(nProcs_),
66     procProcessorPatchStartIndex_(nProcs_),
67     globallySharedPoints_(0),
68     cyclicParallel_(false)
70     if (decompositionDict_.found("distributed"))
71     {
72         distributed_ = Switch(decompositionDict_.lookup("distributed"));
73     }
77 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
79 domainDecomposition::~domainDecomposition()
83 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
85 bool domainDecomposition::writeDecomposition()
87     Info<< "\nConstructing processor meshes" << endl;
89     // Make a lookup map for globally shared points
90     Map<label> sharedPointLookup(2*globallySharedPoints_.size());
92     forAll (globallySharedPoints_, pointi)
93     {
94         sharedPointLookup.insert(globallySharedPoints_[pointi], pointi);
95     }
98     // Mark point/faces/cells that are in zones.  Bad coding - removed
99     // HJ, 31/Mar/2009
101     label totProcFaces = 0;
102     label maxProcPatches = 0;
103     label maxProcFaces = 0;
106     // Write out the meshes
107     for (label procI = 0; procI < nProcs_; procI++)
108     {
109         // Create processor points
110         const labelList& curPointLabels = procPointAddressing_[procI];
112         // Access list of all points in the mesh.  HJ, 27/Mar/2009
113         const pointField& meshPoints = allPoints();
115         labelList pointLookup(meshPoints.size(), -1);
117         pointField procPoints(curPointLabels.size());
119         forAll (curPointLabels, pointi)
120         {
121             procPoints[pointi] = meshPoints[curPointLabels[pointi]];
123             pointLookup[curPointLabels[pointi]] = pointi;
124         }
126         // Create processor faces
127         const labelList& curFaceLabels = procFaceAddressing_[procI];
129         // Access list of all faces in the mesh.  HJ, 27/Mar/2009
130         const faceList& meshFaces = allFaces();
132         labelList faceLookup(meshFaces.size(), -1);
134         faceList procFaces(curFaceLabels.size());
136         forAll (curFaceLabels, facei)
137         {
138             // Mark the original face as used
139             // Remember to decrement the index by one (turning index)
140             // HJ, 5/Dec/2001
141             label curF = mag(curFaceLabels[facei]) - 1;
143             faceLookup[curF] = facei;
145             // get the original face
146             labelList origFaceLabels;
148             if (curFaceLabels[facei] >= 0)
149             {
150                 // face not turned
151                 origFaceLabels = meshFaces[curF];
152             }
153             else
154             {
155                 origFaceLabels = meshFaces[curF].reverseFace();
156             }
158             // translate face labels into local point list
159             face& procFaceLabels = procFaces[facei];
161             procFaceLabels.setSize(origFaceLabels.size());
163             forAll (origFaceLabels, pointi)
164             {
165                 procFaceLabels[pointi] = pointLookup[origFaceLabels[pointi]];
166             }
167         }
169         // Create cell lookup
170         labelList cellLookup(nCells(), -1);
171         const labelList& curCellLabels = procCellAddressing_[procI];
173         forAll (curCellLabels, cellI)
174         {
175             cellLookup[curCellLabels[cellI]] = cellI;
176         }
178         // Get complete owner-neighour addressing in the mesh
179         const labelList& own = faceOwner();
180         const labelList& nei = faceNeighbour();
182         // Calculate owner and neighbour list
183         // Owner list is sized to number of live faces
184         // Neighbour list is sized to number of internal faces
186         labelList procOwner(nLiveProcFaces_[procI]);
188         // Note: loop over owner, not all faces: sizes are different
189         forAll (procOwner, faceI)
190         {
191             // Remember to decrement the index by one (turning index)
192             // HJ, 28/Mar/2009
193             label curF = mag(curFaceLabels[faceI]) - 1;
195             if (curFaceLabels[faceI] >= 0)
196             {
197                 procOwner[faceI] = cellLookup[own[curF]];
198             }
199             else
200             {
201                 procOwner[faceI] = cellLookup[nei[curF]];
202             }
203         }
205         labelList procNeighbour(nInternalProcFaces_[procI]);
207         // Note: loop over neighbour, not all faces: sizes are different
208         forAll (procNeighbour, faceI)
209         {
210             // Remember to decrement the index by one (turning index)
211             // HJ, 28/Mar/2009
212             label curF = mag(curFaceLabels[faceI]) - 1;
214             if (curFaceLabels[faceI] >= 0)
215             {
216                 procNeighbour[faceI] = cellLookup[nei[curF]];
217             }
218             else
219             {
220                 procNeighbour[faceI] = cellLookup[own[curF]];
221             }
222         }
224         // Create processor cells.  No longer needed: using owner and neighbour
225         // HJ, 28/Mar/2009
226 //         const cellList& meshCells = cells();
228 //         cellList procCells(curCellLabels.size());
230 //         forAll (curCellLabels, cellI)
231 //         {
232 //             const labelList& origCellLabels = meshCells[curCellLabels[cellI]];
234 //             cell& curCell = procCells[cellI];
236 //             curCell.setSize(origCellLabels.size());
238 //             forAll (origCellLabels, cellFaceI)
239 //             {
240 //                 curCell[cellFaceI] = faceLookup[origCellLabels[cellFaceI]];
241 //             }
242 //         }
244         // Create processor mesh without a boundary
246         fileName processorCasePath
247         (
248             time().caseName()/fileName(word("processor") + Foam::name(procI))
249         );
251         // make the processor directory
252         mkDir(time().rootPath()/processorCasePath);
254         // create a database
255         Time processorDb
256         (
257             Time::controlDictName,
258             time().rootPath(),
259             processorCasePath,
260             "system",
261             "constant",
262             true
263         );
265         // Create the mesh
266         polyMesh procMesh
267         (
268             IOobject
269             (
270                 this->polyMesh::name(),  // region name of undecomposed mesh
271                 pointsInstance(),
272                 processorDb
273             ),
274             xferMove(procPoints),
275             xferMove(procFaces),
276             xferMove(procOwner),
277             xferMove(procNeighbour),
278             false          // Do not sync par
279 //   xferMove(procCells)   // Old-fashioned mesh creation using cells.
280                            // Deprecated: using face owner/neighbour
281                            // HJ, 30/Mar/2009
282         );
284         // Create processor boundary patches
285         const labelList& curBoundaryAddressing = procBoundaryAddressing_[procI];
287         const labelList& curPatchSizes = procPatchSize_[procI];
289         const labelList& curPatchStarts = procPatchStartIndex_[procI];
291         const labelList& curNeighbourProcessors =
292             procNeighbourProcessors_[procI];
294         const labelList& curProcessorPatchSizes =
295             procProcessorPatchSize_[procI];
297         const labelList& curProcessorPatchStarts =
298             procProcessorPatchStartIndex_[procI];
300         const polyPatchList& meshPatches = boundaryMesh();
302         List<polyPatch*> procPatches
303         (
304             curPatchSizes.size()
305           + curProcessorPatchSizes.size(),
306             reinterpret_cast<polyPatch*>(0)
307         );
309         label nPatches = 0;
311         forAll (curPatchSizes, patchi)
312         {
313             procPatches[nPatches] =
314                 meshPatches[curBoundaryAddressing[patchi]].clone
315                 (
316                     procMesh.boundaryMesh(),
317                     nPatches,
318                     curPatchSizes[patchi],
319                     curPatchStarts[patchi]
320                 ).ptr();
322             nPatches++;
323         }
325         forAll (curProcessorPatchSizes, procPatchI)
326         {
327             procPatches[nPatches] =
328                 new processorPolyPatch
329                 (
330                     word("procBoundary") + Foam::name(procI)
331                   + word("to")
332                   + Foam::name(curNeighbourProcessors[procPatchI]),
333                     curProcessorPatchSizes[procPatchI],
334                     curProcessorPatchStarts[procPatchI],
335                     nPatches,
336                     procMesh.boundaryMesh(),
337                     procI,
338                     curNeighbourProcessors[procPatchI]
339             );
341             nPatches++;
342         }
344         // Add boundary patches
345         procMesh.addPatches(procPatches);
347         // Create and add zones
349         // Note:
350         // This coding was all wrong, as each point/face/cell may only
351         //  belong to a single zone.
352         // Additionally, ordering of points/faces/cells in the processor mesh
353         // needs to match the ordering in global mesh zones.  Full rewrite.
354         // HJ, 30/Mar/2009
356         // Create zones if needed
357         if
358         (
359             pointZones().size() > 0
360          || faceZones().size() > 0
361          || cellZones().size() > 0
362         )
363         {
364             // Point zones
365             List<pointZone*> procPz(pointZones().size());
367             {
368                 const pointZoneMesh& pz = pointZones();
370                 // Go through all the zoned points and find out if they
371                 // belong to a processor.  If so, add it to the zone as
372                 // necessary
373                 forAll (pz, zoneI)
374                 {
375                     const labelList& zonePoints = pz[zoneI];
377                     labelList procZonePoints(zonePoints.size());
378                     label nZonePoints = 0;
380                     forAll (zonePoints, pointI)
381                     {
382                         const label localIndex =
383                             pointLookup[zonePoints[pointI]];
385                         if (localIndex >= 0)
386                         {
387                             // Point live on processor: add to zone
388                             procZonePoints[nZonePoints] = localIndex;
389                             nZonePoints++;
390                         }
391                     }
393                     // Add the zone
394                     procZonePoints.setSize(nZonePoints);
396                     procPz[zoneI] = new pointZone
397                     (
398                         pz[zoneI].name(),
399                         procZonePoints,
400                         zoneI,
401                         procMesh.pointZones()
402                     );
403                 }
404             }
407             // Face zones
408             List<faceZone*> procFz(faceZones().size());
410             {
411                 const faceZoneMesh& fz = faceZones();
413                 forAll (fz, zoneI)
414                 {
415                     const labelList& zoneFaces = fz[zoneI];
416                     const boolList& flipMap = fz[zoneI].flipMap();
418                     // Go through all the zoned faces and find out if they
419                     // belong to a processor.  If so, add it to the zone as
420                     // necessary
422                     labelList procZoneFaces(zoneFaces.size());
423                     boolList procZoneFaceFlips(zoneFaces.size());
424                     label nZoneFaces = 0;
426                     forAll (zoneFaces, faceI)
427                     {
428                         const label localIndex = faceLookup[zoneFaces[faceI]];
430                         if (localIndex >= 0)
431                         {
432                             // Face is present on the processor
434                             // Add the face to the zone
435                             procZoneFaces[nZoneFaces] = localIndex;
437                             // Grab the flip
438                             bool flip = flipMap[faceI];
440                             if (curFaceLabels[localIndex] < 0)
441                             {
442                                 flip = !flip;
443                             }
445                             procZoneFaceFlips[nZoneFaces] = flip;
446                             nZoneFaces++;
447                         }
448                     }
450                     // Add the zone
451                     procZoneFaces.setSize(nZoneFaces);
452                     procZoneFaceFlips.setSize(nZoneFaces);
454                     procFz[zoneI] = new faceZone
455                     (
456                         fz[zoneI].name(),
457                         procZoneFaces,
458                         procZoneFaceFlips,
459                         zoneI,
460                         procMesh.faceZones()
461                     );
462                 }
463             }
465             // Cell zones
466             List<cellZone*> procCz(cellZones().size());
468             {
469                 const cellZoneMesh& cz = cellZones();
471                 // Go through all the zoned cells and find out if they
472                 // belong to a processor.  If so, add it to the zone as
473                 // necessary
475                 forAll (cz, zoneI)
476                 {
477                     const labelList& zoneCells = cz[zoneI];
479                     labelList procZoneCells(zoneCells.size());
480                     label nZoneCells = 0;
482                     forAll (zoneCells, cellI)
483                     {
484                         const label localIndex = cellLookup[zoneCells[cellI]];
486                         if (localIndex >= 0)
487                         {
488                             procZoneCells[nZoneCells] = localIndex;
489                             nZoneCells++;
490                         }
491                     }
493                     // Add the zone
494                     procZoneCells.setSize(nZoneCells);
496                     procCz[zoneI] = new cellZone
497                     (
498                         cz[zoneI].name(),
499                         procZoneCells,
500                         zoneI,
501                         procMesh.cellZones()
502                     );
503                 }
504             }
506             // Add zones
507             procMesh.addZones(procPz, procFz, procCz);
508         }
510         // Set the precision of the points data to 10
511         IOstream::defaultPrecision(10);
513         procMesh.write();
515         Info<< endl
516             << "Processor " << procI << nl
517             << "    Number of cells = " << procMesh.nCells()
518             << endl;
520         label nBoundaryFaces = 0;
521         label nProcPatches = 0;
522         label nProcFaces = 0;
524         forAll (procMesh.boundaryMesh(), patchi)
525         {
526             if
527             (
528                 procMesh.boundaryMesh()[patchi].type()
529              == processorPolyPatch::typeName
530             )
531             {
532                 const processorPolyPatch& ppp =
533                 refCast<const processorPolyPatch>
534                 (
535                     procMesh.boundaryMesh()[patchi]
536                 );
538                 Info<< "    Number of faces shared with processor "
539                     << ppp.neighbProcNo() << " = " << ppp.size() << endl;
541                 nProcPatches++;
542                 nProcFaces += ppp.size();
543             }
544             else
545             {
546                 nBoundaryFaces += procMesh.boundaryMesh()[patchi].size();
547             }
548         }
550         Info<< "    Number of processor patches = " << nProcPatches << nl
551             << "    Number of processor faces = " << nProcFaces << nl
552             << "    Number of boundary faces = " << nBoundaryFaces << endl;
554         totProcFaces += nProcFaces;
555         maxProcPatches = max(maxProcPatches, nProcPatches);
556         maxProcFaces = max(maxProcFaces, nProcFaces);
558         // create and write the addressing information
559         labelIOList pointProcAddressing
560         (
561             IOobject
562             (
563                 "pointProcAddressing",
564                 procMesh.facesInstance(),
565                 procMesh.meshSubDir,
566                 procMesh,
567                 IOobject::NO_READ,
568                 IOobject::NO_WRITE
569             ),
570             procPointAddressing_[procI]
571         );
572         pointProcAddressing.write();
574         labelIOList faceProcAddressing
575         (
576             IOobject
577             (
578                 "faceProcAddressing",
579                 procMesh.facesInstance(),
580                 procMesh.meshSubDir,
581                 procMesh,
582                 IOobject::NO_READ,
583                 IOobject::NO_WRITE
584             ),
585             procFaceAddressing_[procI]
586         );
587         faceProcAddressing.write();
589         labelIOList cellProcAddressing
590         (
591             IOobject
592             (
593                 "cellProcAddressing",
594                 procMesh.facesInstance(),
595                 procMesh.meshSubDir,
596                 procMesh,
597                 IOobject::NO_READ,
598                 IOobject::NO_WRITE
599             ),
600             procCellAddressing_[procI]
601         );
602         cellProcAddressing.write();
604         labelIOList boundaryProcAddressing
605         (
606             IOobject
607             (
608                 "boundaryProcAddressing",
609                 procMesh.facesInstance(),
610                 procMesh.meshSubDir,
611                 procMesh,
612                 IOobject::NO_READ,
613                 IOobject::NO_WRITE
614             ),
615             procBoundaryAddressing_[procI]
616         );
617         boundaryProcAddressing.write();
618     }
620     Info<< nl
621         << "Number of processor faces = " << totProcFaces/2 << nl
622         << "Max number of processor patches = " << maxProcPatches << nl
623         << "Max number of faces between processors = " << maxProcFaces
624         << endl;
626     return true;
630 // ************************************************************************* //