FreeFOAM The Cross-Platform CFD Toolkit
renumberMesh.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 | Copyright (C) 1991-2010 OpenCFD Ltd.
6  \\/ M anispulation |
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 Application
25  renumberMesh
26 
27 Description
28  Renumbers the cell list in order to reduce the bandwidth, reading and
29  renumbering all fields from all the time directories.
30 
31 Usage
32 
33  - renumberMesh [OPTIONS]
34 
35  @param -blockOrder \n
36  Order cells into regions and faces into region-internal/external.
37 
38  @param -writeMaps \n
39  Write renumber map (new to old)
40 
41  @param -orderPoints \n
42  Order points into internal/boundary points.
43 
44  @param -overwrite \n
45  Overwrite existing data.
46 
47  @param -noZero \n
48  Ignore timestep 0.
49 
50  @param -constant \n
51  Include the constant directory.
52 
53  @param -time <time>\n
54  Apply only to specific time.
55 
56  @param -latestTime \n
57  Only apply to latest time step.
58 
59  @param -case <dir>\n
60  Case directory.
61 
62  @param -parallel \n
63  Run in parallel.
64 
65  @param -help \n
66  Display help message.
67 
68  @param -doc \n
69  Display Doxygen API documentation page for this application.
70 
71  @param -srcDoc \n
72  Display Doxygen source documentation page for this application.
73 
74 \*---------------------------------------------------------------------------*/
75 
76 #include <OpenFOAM/argList.H>
77 #include <OpenFOAM/IOobjectList.H>
78 #include <finiteVolume/fvMesh.H>
80 #include <OpenFOAM/ReadFields.H>
81 #include <finiteVolume/volFields.H>
84 #include <meshTools/faceSet.H>
85 #include <OpenFOAM/SortableList.H>
89 
90 using namespace Foam;
91 
92 
93 // Calculate band of matrix
94 label getBand(const labelList& owner, const labelList& neighbour)
95 {
96  label band = 0;
97 
98  forAll(neighbour, faceI)
99  {
100  label diff = neighbour[faceI] - owner[faceI];
101 
102  if (diff > band)
103  {
104  band = diff;
105  }
106  }
107  return band;
108 }
109 
110 
111 // Return new to old cell numbering
112 labelList regionBandCompression
113 (
114  const fvMesh& mesh,
115  const labelList& cellToRegion
116 )
117 {
118  Pout<< "Determining cell order:" << endl;
119 
120  labelList cellOrder(cellToRegion.size());
121 
122  label nRegions = max(cellToRegion)+1;
123 
124  labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
125 
126  label cellI = 0;
127 
128  forAll(regionToCells, regionI)
129  {
130  Pout<< " region " << regionI << " starts at " << cellI << endl;
131 
132  // Per region do a reordering.
133  fvMeshSubset subsetter(mesh);
134  subsetter.setLargeCellSubset(cellToRegion, regionI);
135  const fvMesh& subMesh = subsetter.subMesh();
136  labelList subCellOrder(bandCompression(subMesh.cellCells()));
137 
138  const labelList& cellMap = subsetter.cellMap();
139 
140  forAll(subCellOrder, i)
141  {
142  cellOrder[cellI++] = cellMap[subCellOrder[i]];
143  }
144  }
145  Pout<< endl;
146 
147  return cellOrder;
148 }
149 
150 
151 // Determine face order such that inside region faces are sorted
152 // upper-triangular but inbetween region faces are handled like boundary faces.
153 labelList regionFaceOrder
154 (
155  const primitiveMesh& mesh,
156  const labelList& cellOrder, // new to old cell
157  const labelList& cellToRegion // old cell to region
158 )
159 {
160  Pout<< "Determining face order:" << endl;
161 
162  labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
163 
164  labelList oldToNewFace(mesh.nFaces(), -1);
165 
166  label newFaceI = 0;
167 
168  label prevRegion = -1;
169 
170  forAll (cellOrder, newCellI)
171  {
172  label oldCellI = cellOrder[newCellI];
173 
174  if (cellToRegion[oldCellI] != prevRegion)
175  {
176  prevRegion = cellToRegion[oldCellI];
177  Pout<< " region " << prevRegion << " internal faces start at "
178  << newFaceI << endl;
179  }
180 
181  const cell& cFaces = mesh.cells()[oldCellI];
182 
183  SortableList<label> nbr(cFaces.size());
184 
185  forAll(cFaces, i)
186  {
187  label faceI = cFaces[i];
188 
189  if (mesh.isInternalFace(faceI))
190  {
191  // Internal face. Get cell on other side.
192  label nbrCellI = reverseCellOrder[mesh.faceNeighbour()[faceI]];
193  if (nbrCellI == newCellI)
194  {
195  nbrCellI = reverseCellOrder[mesh.faceOwner()[faceI]];
196  }
197 
198  if (cellToRegion[oldCellI] != cellToRegion[cellOrder[nbrCellI]])
199  {
200  // Treat like external face. Do later.
201  nbr[i] = -1;
202  }
203  else if (newCellI < nbrCellI)
204  {
205  // CellI is master
206  nbr[i] = nbrCellI;
207  }
208  else
209  {
210  // nbrCell is master. Let it handle this face.
211  nbr[i] = -1;
212  }
213  }
214  else
215  {
216  // External face. Do later.
217  nbr[i] = -1;
218  }
219  }
220 
221  nbr.sort();
222 
223  forAll(nbr, i)
224  {
225  if (nbr[i] != -1)
226  {
227  oldToNewFace[cFaces[nbr.indices()[i]]] = newFaceI++;
228  }
229  }
230  }
231 
232  // Do region interfaces
233  label nRegions = max(cellToRegion)+1;
234  {
235  // Sort in increasing region
236  SortableList<label> sortKey(mesh.nFaces(), labelMax);
237 
238  for (label faceI = 0; faceI < mesh.nInternalFaces(); faceI++)
239  {
240  label ownRegion = cellToRegion[mesh.faceOwner()[faceI]];
241  label neiRegion = cellToRegion[mesh.faceNeighbour()[faceI]];
242 
243  if (ownRegion != neiRegion)
244  {
245  sortKey[faceI] =
246  min(ownRegion, neiRegion)*nRegions
247  +max(ownRegion, neiRegion);
248  }
249  }
250  sortKey.sort();
251 
252  // Extract.
253  label prevKey = -1;
254  forAll(sortKey, i)
255  {
256  label key = sortKey[i];
257 
258  if (key == labelMax)
259  {
260  break;
261  }
262 
263  if (prevKey != key)
264  {
265  Pout<< " faces inbetween region " << key/nRegions
266  << " and " << key%nRegions
267  << " start at " << newFaceI << endl;
268  prevKey = key;
269  }
270 
271  oldToNewFace[sortKey.indices()[i]] = newFaceI++;
272  }
273  }
274 
275  // Leave patch faces intact.
276  for (label faceI = newFaceI; faceI < mesh.nFaces(); faceI++)
277  {
278  oldToNewFace[faceI] = faceI;
279  }
280 
281 
282  // Check done all faces.
283  forAll(oldToNewFace, faceI)
284  {
285  if (oldToNewFace[faceI] == -1)
286  {
288  (
289  "polyDualMesh::getFaceOrder"
290  "(const labelList&, const labelList&, const label) const"
291  ) << "Did not determine new position"
292  << " for face " << faceI
293  << abort(FatalError);
294  }
295  }
296  Pout<< endl;
297 
298  return invert(mesh.nFaces(), oldToNewFace);
299 }
300 
301 
302 // cellOrder: old cell for every new cell
303 // faceOrder: old face for every new face. Ordering of boundary faces not
304 // changed.
305 autoPtr<mapPolyMesh> reorderMesh
306 (
307  polyMesh& mesh,
308  const labelList& cellOrder,
309  const labelList& faceOrder
310 )
311 {
312  labelList reverseCellOrder(invert(cellOrder.size(), cellOrder));
313  labelList reverseFaceOrder(invert(faceOrder.size(), faceOrder));
314 
315  faceList newFaces(reorder(reverseFaceOrder, mesh.faces()));
316  labelList newOwner
317  (
318  renumber
319  (
320  reverseCellOrder,
321  reorder(reverseFaceOrder, mesh.faceOwner())
322  )
323  );
324  labelList newNeighbour
325  (
326  renumber
327  (
328  reverseCellOrder,
329  reorder(reverseFaceOrder, mesh.faceNeighbour())
330  )
331  );
332 
333  // Check if any faces need swapping.
334  forAll(newNeighbour, faceI)
335  {
336  label own = newOwner[faceI];
337  label nei = newNeighbour[faceI];
338 
339  if (nei < own)
340  {
341  newFaces[faceI] = newFaces[faceI].reverseFace();
342  Swap(newOwner[faceI], newNeighbour[faceI]);
343  }
344  }
345 
346  const polyBoundaryMesh& patches = mesh.boundaryMesh();
347  labelList patchSizes(patches.size());
348  labelList patchStarts(patches.size());
349  labelList oldPatchNMeshPoints(patches.size());
350  labelListList patchPointMap(patches.size());
351  forAll(patches, patchI)
352  {
353  patchSizes[patchI] = patches[patchI].size();
354  patchStarts[patchI] = patches[patchI].start();
355  oldPatchNMeshPoints[patchI] = patches[patchI].nPoints();
356  patchPointMap[patchI] = identity(patches[patchI].nPoints());
357  }
358 
359  mesh.resetPrimitives
360  (
362  xferMove(newFaces),
363  xferMove(newOwner),
364  xferMove(newNeighbour),
365  patchSizes,
366  patchStarts
367  );
368 
369  return autoPtr<mapPolyMesh>
370  (
371  new mapPolyMesh
372  (
373  mesh, //const polyMesh& mesh,
374  mesh.nPoints(), // nOldPoints,
375  mesh.nFaces(), // nOldFaces,
376  mesh.nCells(), // nOldCells,
377  identity(mesh.nPoints()), // pointMap,
378  List<objectMap>(0), // pointsFromPoints,
379  faceOrder, // faceMap,
380  List<objectMap>(0), // facesFromPoints,
381  List<objectMap>(0), // facesFromEdges,
382  List<objectMap>(0), // facesFromFaces,
383  cellOrder, // cellMap,
384  List<objectMap>(0), // cellsFromPoints,
385  List<objectMap>(0), // cellsFromEdges,
386  List<objectMap>(0), // cellsFromFaces,
387  List<objectMap>(0), // cellsFromCells,
388  identity(mesh.nPoints()), // reversePointMap,
389  reverseFaceOrder, // reverseFaceMap,
390  reverseCellOrder, // reverseCellMap,
391  labelHashSet(0), // flipFaceFlux,
392  patchPointMap, // patchPointMap,
393  labelListList(0), // pointZoneMap,
394  labelListList(0), // faceZonePointMap,
395  labelListList(0), // faceZoneFaceMap,
396  labelListList(0), // cellZoneMap,
397  pointField(0), // preMotionPoints,
398  patchStarts, // oldPatchStarts,
399  oldPatchNMeshPoints // oldPatchNMeshPoints
400  )
401  );
402 }
403 
404 
405 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
406 
407 int main(int argc, char *argv[])
408 {
409  argList::validOptions.insert("blockOrder", "");
410  argList::validOptions.insert("orderPoints", "");
411  argList::validOptions.insert("writeMaps", "");
412  argList::validOptions.insert("overwrite", "");
413 
414 # include <OpenFOAM/addRegionOption.H>
415 # include <OpenFOAM/addTimeOptions.H>
416 
417 # include <OpenFOAM/setRootCase.H>
418 # include <OpenFOAM/createTime.H>
419  runTime.functionObjects().off();
420 
421  // Get times list
422  instantList Times = runTime.times();
423 
424  // set startTime and endTime depending on -time and -latestTime options
425 # include <OpenFOAM/checkTimeOptions.H>
426 
427  runTime.setTime(Times[startTime], startTime);
428 
429 # include <OpenFOAM/createNamedMesh.H>
430  const word oldInstance = mesh.pointsInstance();
431 
432  const bool blockOrder = args.optionFound("blockOrder");
433  if (blockOrder)
434  {
435  Info<< "Ordering cells into regions (using decomposition);"
436  << " ordering faces into region-internal and region-external." << nl
437  << endl;
438  }
439 
440  const bool orderPoints = args.optionFound("orderPoints");
441  if (orderPoints)
442  {
443  Info<< "Ordering points into internal and boundary points." << nl
444  << endl;
445  }
446 
447  const bool writeMaps = args.optionFound("writeMaps");
448 
449  if (writeMaps)
450  {
451  Info<< "Writing renumber maps (new to old) to polyMesh." << nl
452  << endl;
453  }
454 
455  bool overwrite = args.optionFound("overwrite");
456 
457  label band = getBand(mesh.faceOwner(), mesh.faceNeighbour());
458 
459  Info<< "Mesh size: " << returnReduce(mesh.nCells(), sumOp<label>()) << nl
460  << "Band before renumbering: "
461  << returnReduce(band, maxOp<label>()) << nl << endl;
462 
463 
464  // Read parallel reconstruct maps
465  labelIOList cellProcAddressing
466  (
467  IOobject
468  (
469  "cellProcAddressing",
470  mesh.facesInstance(),
472  mesh,
474  ),
475  labelList(0)
476  );
477 
478  labelIOList faceProcAddressing
479  (
480  IOobject
481  (
482  "faceProcAddressing",
483  mesh.facesInstance(),
485  mesh,
487  ),
488  labelList(0)
489  );
490  labelIOList pointProcAddressing
491  (
492  IOobject
493  (
494  "pointProcAddressing",
495  mesh.pointsInstance(),
497  mesh,
499  ),
500  labelList(0)
501  );
502  labelIOList boundaryProcAddressing
503  (
504  IOobject
505  (
506  "boundaryProcAddressing",
507  mesh.pointsInstance(),
509  mesh,
511  ),
512  labelList(0)
513  );
514 
515 
516  // Read objects in time directory
517  IOobjectList objects(mesh, runTime.timeName());
518 
519  // Read vol fields.
520 
522  ReadFields(mesh, objects, vsFlds);
523 
525  ReadFields(mesh, objects, vvFlds);
526 
528  ReadFields(mesh, objects, vstFlds);
529 
530  PtrList<volSymmTensorField> vsymtFlds;
531  ReadFields(mesh, objects, vsymtFlds);
532 
534  ReadFields(mesh, objects, vtFlds);
535 
536  // Read surface fields.
537 
539  ReadFields(mesh, objects, ssFlds);
540 
542  ReadFields(mesh, objects, svFlds);
543 
545  ReadFields(mesh, objects, sstFlds);
546 
548  ReadFields(mesh, objects, ssymtFlds);
549 
551  ReadFields(mesh, objects, stFlds);
552 
553 
555 
556  if (blockOrder)
557  {
558  // Renumbering in two phases. Should be done in one so mapping of
559  // fields is done correctly!
560 
561  // Read decomposePar dictionary
562  IOdictionary decomposeDict
563  (
564  IOobject
565  (
566  "decomposeParDict",
567  runTime.system(),
568  mesh,
571  )
572  );
574  (
575  decomposeDict,
576  mesh
577  );
578 
579  labelList cellToRegion(decomposePtr().decompose(mesh.cellCentres()));
580 
581  // For debugging: write out region
582  {
583  volScalarField cellDist
584  (
585  IOobject
586  (
587  "cellDist",
588  runTime.timeName(),
589  mesh,
592  false
593  ),
594  mesh,
595  dimensionedScalar("cellDist", dimless, 0),
596  zeroGradientFvPatchScalarField::typeName
597  );
598 
599  forAll(cellToRegion, cellI)
600  {
601  cellDist[cellI] = cellToRegion[cellI];
602  }
603 
604  cellDist.write();
605 
606  Info<< nl << "Written decomposition as volScalarField to "
607  << cellDist.name() << " for use in postprocessing."
608  << nl << endl;
609  }
610 
611  // Use block based renumbering.
612  //labelList cellOrder(bandCompression(mesh.cellCells()));
613  labelList cellOrder(regionBandCompression(mesh, cellToRegion));
614 
615  // Determine new to old face order with new cell numbering
616  labelList faceOrder
617  (
618  regionFaceOrder
619  (
620  mesh,
621  cellOrder,
622  cellToRegion
623  )
624  );
625 
626  if (!overwrite)
627  {
628  runTime++;
629  }
630 
631  // Change the mesh.
632  map = reorderMesh(mesh, cellOrder, faceOrder);
633  }
634  else
635  {
636  // Use built-in renumbering.
637 
638  polyTopoChange meshMod(mesh);
639 
640  if (!overwrite)
641  {
642  runTime++;
643  }
644 
645  // Change the mesh.
646  map = meshMod.changeMesh
647  (
648  mesh,
649  false, // inflate
650  true, // parallel sync
651  true, // cell ordering
652  orderPoints // point ordering
653  );
654  }
655 
656  // Update fields
657  mesh.updateMesh(map);
658 
659  // Update proc maps
660  if (cellProcAddressing.headerOk())
661  {
662  Info<< "Renumbering processor cell decomposition map "
663  << cellProcAddressing.name() << endl;
664 
665  cellProcAddressing = labelList
666  (
667  UIndirectList<label>(cellProcAddressing, map().cellMap())
668  );
669  }
670  if (faceProcAddressing.headerOk())
671  {
672  Info<< "Renumbering processor face decomposition map "
673  << faceProcAddressing.name() << endl;
674 
675  faceProcAddressing = labelList
676  (
677  UIndirectList<label>(faceProcAddressing, map().faceMap())
678  );
679  }
680  if (pointProcAddressing.headerOk())
681  {
682  Info<< "Renumbering processor point decomposition map "
683  << pointProcAddressing.name() << endl;
684 
685  pointProcAddressing = labelList
686  (
687  UIndirectList<label>(pointProcAddressing, map().pointMap())
688  );
689  }
690 
691 
692  // Move mesh (since morphing might not do this)
693  if (map().hasMotionPoints())
694  {
695  mesh.movePoints(map().preMotionPoints());
696  }
697 
698 
699  band = getBand(mesh.faceOwner(), mesh.faceNeighbour());
700 
701  Info<< "Band after renumbering: "
702  << returnReduce(band, maxOp<label>()) << nl << endl;
703 
704 
705  if (orderPoints)
706  {
707  // Force edge calculation (since only reason that points would need to
708  // be sorted)
709  (void)mesh.edges();
710 
711  label nTotPoints = returnReduce
712  (
713  mesh.nPoints(),
714  sumOp<label>()
715  );
716  label nTotIntPoints = returnReduce
717  (
718  mesh.nInternalPoints(),
719  sumOp<label>()
720  );
721 
722  label nTotEdges = returnReduce
723  (
724  mesh.nEdges(),
725  sumOp<label>()
726  );
727  label nTotIntEdges = returnReduce
728  (
729  mesh.nInternalEdges(),
730  sumOp<label>()
731  );
732  label nTotInt0Edges = returnReduce
733  (
734  mesh.nInternal0Edges(),
735  sumOp<label>()
736  );
737  label nTotInt1Edges = returnReduce
738  (
739  mesh.nInternal1Edges(),
740  sumOp<label>()
741  );
742 
743  Info<< "Points:" << nl
744  << " total : " << nTotPoints << nl
745  << " internal: " << nTotIntPoints << nl
746  << " boundary: " << nTotPoints-nTotIntPoints << nl
747  << "Edges:" << nl
748  << " total : " << nTotEdges << nl
749  << " internal: " << nTotIntEdges << nl
750  << " internal using 0 boundary points: "
751  << nTotInt0Edges << nl
752  << " internal using 1 boundary points: "
753  << nTotInt1Edges-nTotInt0Edges << nl
754  << " internal using 2 boundary points: "
755  << nTotIntEdges-nTotInt1Edges << nl
756  << " boundary: " << nTotEdges-nTotIntEdges << nl
757  << endl;
758  }
759 
760 
761  if (overwrite)
762  {
763  mesh.setInstance(oldInstance);
764  }
765 
766  Info<< "Writing mesh to " << mesh.facesInstance() << endl;
767 
768  mesh.write();
769  if (cellProcAddressing.headerOk())
770  {
771  cellProcAddressing.instance() = mesh.facesInstance();
772  cellProcAddressing.write();
773  }
774  if (faceProcAddressing.headerOk())
775  {
776  faceProcAddressing.instance() = mesh.facesInstance();
777  faceProcAddressing.write();
778  }
779  if (pointProcAddressing.headerOk())
780  {
781  pointProcAddressing.instance() = mesh.facesInstance();
782  pointProcAddressing.write();
783  }
784  if (boundaryProcAddressing.headerOk())
785  {
786  boundaryProcAddressing.instance() = mesh.facesInstance();
787  boundaryProcAddressing.write();
788  }
789 
790 
791  if (writeMaps)
792  {
794  (
795  IOobject
796  (
797  "cellMap",
798  mesh.facesInstance(),
800  mesh,
803  false
804  ),
805  map().cellMap()
806  ).write();
808  (
809  IOobject
810  (
811  "faceMap",
812  mesh.facesInstance(),
814  mesh,
817  false
818  ),
819  map().faceMap()
820  ).write();
822  (
823  IOobject
824  (
825  "pointMap",
826  mesh.facesInstance(),
828  mesh,
831  false
832  ),
833  map().pointMap()
834  ).write();
835  }
836 
837  Info<< "\nEnd.\n" << endl;
838 
839  return 0;
840 }
841 
842 
843 // ************************ vim: set sw=4 sts=4 et: ************************ //