FreeFOAM The Cross-Platform CFD Toolkit
reconstructParMesh.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 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 Application
25  reconstructParMesh
26 
27 Description
28  Reconstructs a mesh using geometric information only.
29 
30  Writes point/face/cell procAddressing so afterwards reconstructPar can be
31  used to reconstruct fields.
32 
33 Note
34  Uses geometric matching tolerance (set with -mergeTol option)
35 
36  If the parallel case does not have correct procBoundaries use the
37  -fullMatch option which will check all boundary faces (bit slower).
38 
39 Usage
40 
41  - reconstructParMesh [OPTIONS]
42 
43  @param -fullMatch \n
44  Check all boundary faces.
45 
46  @param -mergeTol <number>\n
47  Relative merge distance.
48 
49  @param -region <name>\n
50  Only apply to named mesh region.
51 
52  @param -noZero \n
53  Ignore timestep 0.
54 
55  @param -constant \n
56  Include the constant directory.
57 
58  @param -time <time>\n
59  Apply only to specific time.
60 
61  @param -latestTime \n
62  Only apply to latest time step.
63 
64  @param -case <dir>\n
65  Case directory.
66 
67  @param -help \n
68  Display help message.
69 
70  @param -doc \n
71  Display Doxygen API documentation page for this application.
72 
73  @param -srcDoc \n
74  Display Doxygen source documentation page for this application.
75 
76 \*---------------------------------------------------------------------------*/
77 
78 #include <OpenFOAM/argList.H>
79 #include <OpenFOAM/Time.H>
80 #include <OpenFOAM/IOobjectList.H>
81 #include <OpenFOAM/labelIOList.H>
88 
89 using namespace Foam;
90 
91 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
92 
93 // Tolerance (as fraction of the bounding box). Needs to be fairly lax since
94 // usually meshes get written with limited precision (6 digits)
95 static const scalar defaultMergeTol = 1E-7;
96 
97 
98 static void renumber
99 (
100  const labelList& map,
101  labelList& elems
102 )
103 {
104  forAll(elems, i)
105  {
106  if (elems[i] >= 0)
107  {
108  elems[i] = map[elems[i]];
109  }
110  }
111 }
112 
113 
114 // Determine which faces are coupled. Uses geometric merge distance.
115 // Looks either at all boundaryFaces (fullMatch) or only at the
116 // procBoundaries for procI. Assumes that masterMesh contains already merged
117 // all the processors < procI.
118 autoPtr<faceCoupleInfo> determineCoupledFaces
119 (
120  const bool fullMatch,
121  const label procI,
122  const polyMesh& masterMesh,
123  const polyMesh& meshToAdd,
124  const scalar mergeDist
125 )
126 {
127  if (fullMatch || masterMesh.nCells() == 0)
128  {
130  (
131  new faceCoupleInfo
132  (
133  masterMesh,
134  meshToAdd,
135  mergeDist, // absolute merging distance
136  true // matching faces identical
137  )
138  );
139  }
140  else
141  {
142  // Pick up all patches on masterMesh ending in "toDDD" where DDD is
143  // the processor number procI.
144 
145  const polyBoundaryMesh& masterPatches = masterMesh.boundaryMesh();
146 
147  const string toProcString("to" + name(procI));
148 
149  DynamicList<label> masterFaces
150  (
151  masterMesh.nFaces()
152  - masterMesh.nInternalFaces()
153  );
154 
155  forAll(masterPatches, patchI)
156  {
157  const polyPatch& pp = masterPatches[patchI];
158 
159  if
160  (
161  isA<processorPolyPatch>(pp)
162  && (
163  pp.name().rfind(toProcString)
164  == (pp.name().size()-toProcString.size())
165  )
166  )
167  {
168  label meshFaceI = pp.start();
169  forAll(pp, i)
170  {
171  masterFaces.append(meshFaceI++);
172  }
173  }
174  }
175  masterFaces.shrink();
176 
177 
178  // Pick up all patches on meshToAdd ending in "procBoundaryDDDtoYYY"
179  // where DDD is the processor number procI and YYY is < procI.
180 
181  const polyBoundaryMesh& addPatches = meshToAdd.boundaryMesh();
182 
183  DynamicList<label> addFaces
184  (
185  meshToAdd.nFaces()
186  - meshToAdd.nInternalFaces()
187  );
188 
189  forAll(addPatches, patchI)
190  {
191  const polyPatch& pp = addPatches[patchI];
192 
193  if (isA<processorPolyPatch>(pp))
194  {
195  bool isConnected = false;
196 
197  for (label mergedProcI = 0; mergedProcI < procI; mergedProcI++)
198  {
199  const string fromProcString
200  (
201  "procBoundary"
202  + name(procI)
203  + "to"
204  + name(mergedProcI)
205  );
206 
207  if (pp.name() == fromProcString)
208  {
209  isConnected = true;
210  break;
211  }
212  }
213 
214  if (isConnected)
215  {
216  label meshFaceI = pp.start();
217  forAll(pp, i)
218  {
219  addFaces.append(meshFaceI++);
220  }
221  }
222  }
223  }
224  addFaces.shrink();
225 
227  (
228  new faceCoupleInfo
229  (
230  masterMesh,
231  masterFaces,
232  meshToAdd,
233  addFaces,
234  mergeDist, // absolute merging distance
235  true, // matching faces identical?
236  false, // if perfectmatch are faces already ordered
237  // (e.g. processor patches)
238  false // are faces each on separate patch?
239  )
240  );
241  }
242 }
243 
244 
245 autoPtr<mapPolyMesh> mergeSharedPoints
246 (
247  const scalar mergeDist,
248  polyMesh& mesh,
249  labelListList& pointProcAddressing
250 )
251 {
252  // Find out which sets of points get merged and create a map from
253  // mesh point to unique point.
254  Map<label> pointToMaster
255  (
257  (
258  mesh,
259  mergeDist
260  )
261  );
262 
263  Info<< "mergeSharedPoints : detected " << pointToMaster.size()
264  << " points that are to be merged." << endl;
265 
266  if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
267  {
268  return autoPtr<mapPolyMesh>(NULL);
269  }
270 
271  polyTopoChange meshMod(mesh);
272 
273  fvMeshAdder::mergePoints(mesh, pointToMaster, meshMod);
274 
275  // Change the mesh (no inflation). Note: parallel comms allowed.
276  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
277 
278  // Update fields. No inflation, parallel sync.
279  mesh.updateMesh(map);
280 
281  // pointProcAddressing give indices into the master mesh so adapt them
282  // for changed point numbering.
283 
284  // Adapt constructMaps for merged points.
285  forAll(pointProcAddressing, procI)
286  {
287  labelList& constructMap = pointProcAddressing[procI];
288 
289  forAll(constructMap, i)
290  {
291  label oldPointI = constructMap[i];
292 
293  // New label of point after changeMesh.
294  label newPointI = map().reversePointMap()[oldPointI];
295 
296  if (newPointI < -1)
297  {
298  constructMap[i] = -newPointI-2;
299  }
300  else if (newPointI >= 0)
301  {
302  constructMap[i] = newPointI;
303  }
304  else
305  {
306  FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
307  << "Problem. oldPointI:" << oldPointI
308  << " newPointI:" << newPointI << abort(FatalError);
309  }
310  }
311  }
312 
313  return map;
314 }
315 
316 
317 int main(int argc, char *argv[])
318 {
320  argList::validOptions.insert("mergeTol", "relative merge distance");
321  argList::validOptions.insert("fullMatch", "");
322 
323 # include <OpenFOAM/addTimeOptions.H>
324 # include <OpenFOAM/addRegionOption.H>
325 # include <OpenFOAM/setRootCase.H>
326 # include <OpenFOAM/createTime.H>
327 
328  Info<< "This is an experimental tool which tries to merge"
329  << " individual processor" << nl
330  << "meshes back into one master mesh. Use it if the original"
331  << " master mesh has" << nl
332  << "been deleted or if the processor meshes have been modified"
333  << " (topology change)." << nl
334  << "This tool will write the resulting mesh to a new time step"
335  << " and construct" << nl
336  << "xxxxProcAddressing files in the processor meshes so"
337  << " reconstructPar can be" << nl
338  << "used to regenerate the fields on the master mesh." << nl
339  << nl
340  << "Not well tested & use at your own risk!" << nl
341  << endl;
342 
343 
345  fileName regionPrefix = "";
346  if (args.optionFound("region"))
347  {
348  regionName = args.option("region");
349  regionPrefix = regionName;
350  Info<< "Operating on region " << regionName << nl << endl;
351  }
352 
353  scalar mergeTol = defaultMergeTol;
354  args.optionReadIfPresent("mergeTol", mergeTol);
355 
356  scalar writeTol = Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
357 
358  Info<< "Merge tolerance : " << mergeTol << nl
359  << "Write tolerance : " << writeTol << endl;
360 
361  if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
362  {
364  << "Your current settings specify ASCII writing with "
365  << IOstream::defaultPrecision() << " digits precision." << endl
366  << "Your merging tolerance (" << mergeTol << ") is finer than this."
367  << endl
368  << "Please change your writeFormat to binary"
369  << " or increase the writePrecision" << endl
370  << "or adjust the merge tolerance (-mergeTol)."
371  << exit(FatalError);
372  }
373 
374 
375  const bool fullMatch = args.optionFound("fullMatch");
376 
377  if (fullMatch)
378  {
379  Info<< "Doing geometric matching on all boundary faces." << nl << endl;
380  }
381  else
382  {
383  Info<< "Doing geometric matching on correct procBoundaries only."
384  << nl << "This assumes a correct decomposition." << endl;
385  }
386 
387 
388 
389  int nProcs = 0;
390 
391  while
392  (
393  isDir
394  (
395  args.rootPath()
396  / args.caseName()
397  / fileName(word("processor") + name(nProcs))
398  )
399  )
400  {
401  nProcs++;
402  }
403 
404  Info<< "Found " << nProcs << " processor directories" << nl << endl;
405 
406 
407  // Read all databases.
408  PtrList<Time> databases(nProcs);
409 
410  forAll (databases, procI)
411  {
412  Info<< "Reading database "
413  << args.caseName()/fileName(word("processor") + name(procI))
414  << endl;
415 
416  databases.set
417  (
418  procI,
419  new Time
420  (
422  args.rootPath(),
423  args.caseName()/fileName(word("processor") + name(procI))
424  )
425  );
426 
427  Time& procTime = databases[procI];
428 
429  instantList Times = procTime.times();
430 
431  // set startTime and endTime depending on -time and -latestTime options
432 # include <OpenFOAM/checkTimeOptions.H>
433 
434  procTime.setTime(Times[startTime], startTime);
435 
436  if (procI > 0 && databases[procI-1].value() != procTime.value())
437  {
439  << "Time not equal on processors." << nl
440  << "Processor:" << procI-1
441  << " time:" << databases[procI-1].value() << nl
442  << "Processor:" << procI
443  << " time:" << procTime.value()
444  << exit(FatalError);
445  }
446  }
447 
448  // Set master time
449  Info<< "Setting master time to " << databases[0].timeName() << nl << endl;
450  runTime.setTime(databases[0]);
451 
452 
453  // Read point on individual processors to determine merge tolerance
454  // (otherwise single cell domains might give problems)
455 
457 
458  for (label procI = 0; procI < nProcs; procI++)
459  {
460  fileName pointsInstance
461  (
462  databases[procI].findInstance
463  (
464  regionPrefix/polyMesh::meshSubDir,
465  "points"
466  )
467  );
468 
469  if (pointsInstance != databases[procI].timeName())
470  {
472  << "Your time was specified as " << databases[procI].timeName()
473  << " but there is no polyMesh/points in that time." << endl
474  << "(there is a points file in " << pointsInstance
475  << ")" << endl
476  << "Please rerun with the correct time specified"
477  << " (through the -constant, -time or -latestTime option)."
478  << endl << exit(FatalError);
479  }
480 
481  Info<< "Reading points from "
482  << databases[procI].caseName()
483  << " for time = " << databases[procI].timeName()
484  << nl << endl;
485 
487  (
488  IOobject
489  (
490  "points",
491  databases[procI].findInstance
492  (
493  regionPrefix/polyMesh::meshSubDir,
494  "points"
495  ),
496  regionPrefix/polyMesh::meshSubDir,
497  databases[procI],
500  false
501  )
502  );
503 
504  boundBox domainBb(points, false);
505 
506  bb.min() = min(bb.min(), domainBb.min());
507  bb.max() = max(bb.max(), domainBb.max());
508  }
509  const scalar mergeDist = mergeTol * bb.mag();
510 
511  Info<< "Overall mesh bounding box : " << bb << nl
512  << "Relative tolerance : " << mergeTol << nl
513  << "Absolute matching distance : " << mergeDist << nl
514  << endl;
515 
516 
517  // Addressing from processor to reconstructed case
518  labelListList cellProcAddressing(nProcs);
519  labelListList faceProcAddressing(nProcs);
520  labelListList pointProcAddressing(nProcs);
521  labelListList boundaryProcAddressing(nProcs);
522 
523  // Internal faces on the final reconstructed mesh
524  label masterInternalFaces;
525  // Owner addressing on the final reconstructed mesh
526  labelList masterOwner;
527 
528  {
529  // Construct empty mesh.
530  Info<< "Constructing empty mesh to add to." << nl << endl;
531  polyMesh masterMesh
532  (
533  IOobject
534  (
535  regionName,
536  runTime.timeName(),
537  runTime,
539  ),
540  xferCopy(pointField()),
541  xferCopy(faceList()),
542  xferCopy(cellList())
543  );
544 
545  for (label procI = 0; procI < nProcs; procI++)
546  {
547  Info<< "Reading mesh to add from "
548  << databases[procI].caseName()
549  << " for time = " << databases[procI].timeName()
550  << nl << endl;
551 
552  polyMesh meshToAdd
553  (
554  IOobject
555  (
556  regionName,
557  databases[procI].timeName(),
558  databases[procI]
559  )
560  );
561 
562  // Initialize its addressing
563  cellProcAddressing[procI] = identity(meshToAdd.nCells());
564  faceProcAddressing[procI] = identity(meshToAdd.nFaces());
565  pointProcAddressing[procI] = identity(meshToAdd.nPoints());
566  boundaryProcAddressing[procI] =
567  identity(meshToAdd.boundaryMesh().size());
568 
569 
570  // Find geometrically shared points/faces.
571  autoPtr<faceCoupleInfo> couples = determineCoupledFaces
572  (
573  fullMatch,
574  procI,
575  masterMesh,
576  meshToAdd,
577  mergeDist
578  );
579 
580 
581  // Add elements to mesh
582  Info<< "Adding to master mesh" << nl << endl;
583 
585  (
586  masterMesh,
587  meshToAdd,
588  couples
589  );
590 
591  // Update all addressing so xxProcAddressing points to correct item
592  // in masterMesh.
593 
594  // Processors that were already in masterMesh
595  for (label mergedI = 0; mergedI < procI; mergedI++)
596  {
597  renumber(map().oldCellMap(), cellProcAddressing[mergedI]);
598  renumber(map().oldFaceMap(), faceProcAddressing[mergedI]);
599  renumber(map().oldPointMap(), pointProcAddressing[mergedI]);
600  // Note: boundary is special since can contain -1.
601  renumber(map().oldPatchMap(), boundaryProcAddressing[mergedI]);
602  }
603 
604  // Added processor
605  renumber(map().addedCellMap(), cellProcAddressing[procI]);
606  renumber(map().addedFaceMap(), faceProcAddressing[procI]);
607  renumber(map().addedPointMap(), pointProcAddressing[procI]);
608  renumber(map().addedPatchMap(), boundaryProcAddressing[procI]);
609 
610  Info<< endl;
611  }
612 
613  // See if any points on the mastermesh have become connected
614  // because of connections through processor meshes.
615  mergeSharedPoints(mergeDist, masterMesh, pointProcAddressing);
616 
617  // Save some properties on the reconstructed mesh
618  masterInternalFaces = masterMesh.nInternalFaces();
619  masterOwner = masterMesh.faceOwner();
620 
621 
622  Info<< "\nWriting merged mesh to "
623  << runTime.path()/runTime.timeName()
624  << nl << endl;
625 
626  if (!masterMesh.write())
627  {
629  << "Failed writing polyMesh."
630  << exit(FatalError);
631  }
632  }
633 
634 
635  // Write the addressing
636 
637  Info<< "Reconstructing the addressing from the processor meshes"
638  << " to the newly reconstructed mesh" << nl << endl;
639 
640  forAll(databases, procI)
641  {
642  Info<< "Reading processor " << procI << " mesh from "
643  << databases[procI].caseName() << endl;
644 
645  polyMesh procMesh
646  (
647  IOobject
648  (
649  regionName,
650  databases[procI].timeName(),
651  databases[procI]
652  )
653  );
654 
655 
656  // From processor point to reconstructed mesh point
657 
658  Info<< "Writing pointProcAddressing to "
659  << databases[procI].caseName()
660  /procMesh.facesInstance()
662  << endl;
663 
665  (
666  IOobject
667  (
668  "pointProcAddressing",
669  procMesh.facesInstance(),
671  procMesh,
674  false // do not register
675  ),
676  pointProcAddressing[procI]
677  ).write();
678 
679 
680  // From processor face to reconstructed mesh face
681 
682  Info<< "Writing faceProcAddressing to "
683  << databases[procI].caseName()
684  /procMesh.facesInstance()
686  << endl;
687 
688  labelIOList faceProcAddr
689  (
690  IOobject
691  (
692  "faceProcAddressing",
693  procMesh.facesInstance(),
695  procMesh,
698  false // do not register
699  ),
700  faceProcAddressing[procI]
701  );
702 
703  // Now add turning index to faceProcAddressing.
704  // See reconstrurPar for meaning of turning index.
705  forAll(faceProcAddr, procFaceI)
706  {
707  label masterFaceI = faceProcAddr[procFaceI];
708 
709  if
710  (
711  !procMesh.isInternalFace(procFaceI)
712  && masterFaceI < masterInternalFaces
713  )
714  {
715  // proc face is now external but used to be internal face.
716  // Check if we have owner or neighbour.
717 
718  label procOwn = procMesh.faceOwner()[procFaceI];
719  label masterOwn = masterOwner[masterFaceI];
720 
721  if (cellProcAddressing[procI][procOwn] == masterOwn)
722  {
723  // No turning. Offset by 1.
724  faceProcAddr[procFaceI]++;
725  }
726  else
727  {
728  // Turned face.
729  faceProcAddr[procFaceI] =
730  -1 - faceProcAddr[procFaceI];
731  }
732  }
733  else
734  {
735  // No turning. Offset by 1.
736  faceProcAddr[procFaceI]++;
737  }
738  }
739 
740  faceProcAddr.write();
741 
742 
743  // From processor cell to reconstructed mesh cell
744 
745  Info<< "Writing cellProcAddressing to "
746  << databases[procI].caseName()
747  /procMesh.facesInstance()
749  << endl;
750 
752  (
753  IOobject
754  (
755  "cellProcAddressing",
756  procMesh.facesInstance(),
758  procMesh,
761  false // do not register
762  ),
763  cellProcAddressing[procI]
764  ).write();
765 
766 
767 
768  // From processor patch to reconstructed mesh patch
769 
770  Info<< "Writing boundaryProcAddressing to "
771  << databases[procI].caseName()
772  /procMesh.facesInstance()
774  << endl;
775 
777  (
778  IOobject
779  (
780  "boundaryProcAddressing",
781  procMesh.facesInstance(),
783  procMesh,
786  false // do not register
787  ),
788  boundaryProcAddressing[procI]
789  ).write();
790 
791  Info<< endl;
792  }
793 
794  Info<< "End.\n" << endl;
795 
796  return 0;
797 }
798 
799 
800 // ************************ vim: set sw=4 sts=4 et: ************************ //