conformalVoronoiMeshIO.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) 2012-2019 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 \*---------------------------------------------------------------------------*/
25 
26 #include "conformalVoronoiMesh.H"
27 #include "IOstreams.H"
28 #include "OFstream.H"
29 #include "pointMesh.H"
30 #include "pointFields.H"
31 #include "ListOps.H"
32 #include "polyMeshFilter.H"
33 #include "polyTopoChange.H"
34 #include "PrintTable.H"
35 #include "pointMesh.H"
36 #include "indexedVertexOps.H"
37 #include "DelaunayMeshTools.H"
38 #include "syncTools.H"
39 #include "faceSet.H"
40 #include "OBJstream.H"
41 
42 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
43 
45 (
46  const string& description
47 ) const
48 {
49  timeCheck(time(), description, foamyHexMeshControls().timeChecks());
50 }
51 
52 
54 (
55  const Time& runTime,
56  const string& description,
57  const bool check
58 )
59 {
60  if (check)
61  {
62  Info<< nl << "--- [ cpuTime "
63  << runTime.elapsedCpuTime() << " s, "
64  << "delta " << runTime.cpuTimeIncrement()<< " s";
65 
66  if (description != word::null)
67  {
68  Info<< ", " << description << " ";
69  }
70  else
71  {
72  Info<< " ";
73  }
74 
75  Info<< "] --- " << endl;
76 
77  memInfo m;
78 
79  if (m.valid())
80  {
81  PrintTable<word, label> memoryTable
82  (
83  "Memory Usage (kB): "
84  + description
85  );
86 
87  memoryTable.add("mSize", m.size());
88  memoryTable.add("mPeak", m.peak());
89  memoryTable.add("mRss", m.rss());
90 
91  Info<< incrIndent;
92  memoryTable.print(Info, true, true);
93  Info<< decrIndent;
94  }
95  }
96 }
97 
98 
99 void Foam::conformalVoronoiMesh::writeMesh(const fileName& instance)
100 {
102 
103  // Per cell the Delaunay vertex
104  labelList cellToDelaunayVertex;
105  // Per patch, per face the Delaunay vertex
106  labelListList patchToDelaunayVertex;
107 
108  labelList dualPatchStarts;
109 
110  {
112  labelList boundaryPts;
113  faceList faces;
114  labelList owner;
115  labelList neighbour;
117  PtrList<dictionary> patchDicts;
118  pointField cellCentres;
119  PackedBoolList boundaryFacesToRemove;
120 
121  calcDualMesh
122  (
123  points,
124  boundaryPts,
125  faces,
126  owner,
127  neighbour,
128  patchNames,
129  patchDicts,
130  cellCentres,
131  cellToDelaunayVertex,
132  patchToDelaunayVertex,
133  boundaryFacesToRemove
134  );
135 
136  Info<< nl << "Writing polyMesh to " << instance << endl;
137 
138  writeMesh
139  (
141  instance,
142  points,
143  boundaryPts,
144  faces,
145  owner,
146  neighbour,
147  patchNames,
148  patchDicts,
149  cellCentres,
150  boundaryFacesToRemove
151  );
152 
153  dualPatchStarts.setSize(patchDicts.size());
154 
155  forAll(dualPatchStarts, patchi)
156  {
157  dualPatchStarts[patchi] =
158  readLabel(patchDicts[patchi].lookup("startFace"));
159  }
160  }
161 
162  if (foamyHexMeshControls().writeCellShapeControlMesh())
163  {
165  }
166 
167  if (foamyHexMeshControls().writeBackgroundMeshDecomposition())
168  {
169  Info<< nl << "Writing " << "backgroundMeshDecomposition" << endl;
170 
171  // Have to explicitly update the mesh instance.
172  const_cast<fvMesh&>(decomposition_().mesh()).setInstance
173  (
174  time().timeName()
175  );
176 
177  decomposition_().mesh().write();
178  }
179 
180  if (foamyHexMeshControls().writeTetDualMesh())
181  {
182  label celli = 0;
183  for
184  (
185  Finite_cells_iterator cit = finite_cells_begin();
186  cit != finite_cells_end();
187  ++cit
188  )
189  {
190  if
191  (
192  !cit->hasFarPoint()
193  && !is_infinite(cit)
194  )
195  {
196  cit->cellIndex() = celli++;
197  }
198  }
199 
200  Info<< nl << "Writing " << "tetDualMesh" << endl;
201 
203  labelList cellMap;
204  autoPtr<polyMesh> tetMesh =
205  createMesh("tetDualMesh", vertexMap, cellMap);
206 
207  tetMesh().write();
208 
209 // // Determine map from Delaunay vertex to Dual mesh
210 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
211 //
212 // // From all Delaunay vertices to cell (positive index)
213 // // or patch face (negative index)
214 // labelList vertexToDualAddressing(number_of_vertices(), 0);
215 //
216 // forAll(cellToDelaunayVertex, celli)
217 // {
218 // label vertI = cellToDelaunayVertex[celli];
219 //
220 // if (vertexToDualAddressing[vertI] != 0)
221 // {
222 // FatalErrorInFunction
223 // << "Delaunay vertex " << vertI
224 // << " from cell " << celli
225 // << " is already mapped to "
226 // << vertexToDualAddressing[vertI]
227 // << exit(FatalError);
228 // }
229 // vertexToDualAddressing[vertI] = celli+1;
230 // }
231 //
232 // forAll(patchToDelaunayVertex, patchi)
233 // {
234 // const labelList& patchVertices = patchToDelaunayVertex[patchi];
235 //
236 // forAll(patchVertices, i)
237 // {
238 // label vertI = patchVertices[i];
239 //
240 // if (vertexToDualAddressing[vertI] > 0)
241 // {
242 // FatalErrorInFunction
243 // << "Delaunay vertex " << vertI
244 // << " from patch " << patchi
245 // << " local index " << i
246 // << " is already mapped to cell "
247 // << vertexToDualAddressing[vertI]-1
248 // << exit(FatalError);
249 // }
250 //
251 // // Vertex might be used by multiple faces. Which one to
252 // // use? For now last one wins.
253 // label dualFacei = dualPatchStarts[patchi]+i;
254 // vertexToDualAddressing[vertI] = -dualFacei-1;
255 // }
256 // }
257 //
258 //
259 // // Calculate tet mesh addressing
260 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
261 //
262 // pointField points;
263 // labelList boundaryPts(number_of_finite_cells(), -1);
264 // // From tet point back to Delaunay vertex index
265 // labelList pointToDelaunayVertex;
266 // faceList faces;
267 // labelList owner;
268 // labelList neighbour;
269 // wordList patchTypes;
270 // wordList patchNames;
271 // PtrList<dictionary> patchDicts;
272 // pointField cellCentres;
273 //
274 // calcTetMesh
275 // (
276 // points,
277 // pointToDelaunayVertex,
278 // faces,
279 // owner,
280 // neighbour,
281 // patchTypes,
282 // patchNames,
283 // patchDicts
284 // );
285 //
286 //
287 //
288 // // Calculate map from tet points to dual mesh cells/patch faces
289 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
290 //
291 // labelIOList pointDualAddressing
292 // (
293 // IOobject
294 // (
295 // "pointDualAddressing",
296 // instance,
297 // "tetDualMesh"/polyMesh::meshSubDir,
298 // runTime_,
299 // IOobject::NO_READ,
300 // IOobject::AUTO_WRITE,
301 // false
302 // ),
303 // UIndirectList<label>
304 // (
305 // vertexToDualAddressing,
306 // pointToDelaunayVertex
307 // )()
308 // );
309 //
310 // label pointi = findIndex(pointDualAddressing, -1);
311 // if (pointi != -1)
312 // {
313 // WarningInFunction
314 // << "Delaunay vertex " << pointi
315 // << " does not have a corresponding dual cell." << endl;
316 // }
317 //
318 // Info<< "Writing map from tetDualMesh points to Voronoi mesh to "
319 // << pointDualAddressing.objectPath() << endl;
320 // pointDualAddressing.write();
321 //
322 //
323 //
324 // // Write tet points corresponding to the Voronoi cell/face centre
325 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
326 // {
327 // // Read Voronoi mesh
328 // fvMesh mesh
329 // (
330 // IOobject
331 // (
332 // Foam::polyMesh::defaultRegion,
333 // instance,
334 // runTime_,
335 // IOobject::MUST_READ
336 // )
337 // );
338 // pointIOField dualPoints
339 // (
340 // IOobject
341 // (
342 // "dualPoints",
343 // instance,
344 // "tetDualMesh"/polyMesh::meshSubDir,
345 // runTime_,
346 // IOobject::NO_READ,
347 // IOobject::AUTO_WRITE,
348 // false
349 // ),
350 // points
351 // );
352 //
353 // forAll(pointDualAddressing, pointi)
354 // {
355 // label index = pointDualAddressing[pointi];
356 //
357 // if (index > 0)
358 // {
359 // label celli = index-1;
360 // dualPoints[pointi] = mesh.cellCentres()[celli];
361 // }
362 // else if (index < 0)
363 // {
364 // label facei = -index-1;
365 // if (facei >= mesh.nInternalFaces())
366 // {
367 // dualPoints[pointi] = mesh.faceCentres()[facei];
368 // }
369 // }
370 // }
371 //
372 // Info<< "Writing tetDualMesh points mapped onto Voronoi mesh to "
373 // << dualPoints.objectPath() << endl
374 // << "Replace the polyMesh/points with these." << endl;
375 // dualPoints.write();
376 // }
377 //
378 //
379 // Info<< nl << "Writing tetDualMesh to " << instance << endl;
380 //
381 // PackedBoolList boundaryFacesToRemove;
382 // writeMesh
383 // (
384 // "tetDualMesh",
385 // instance,
386 // points,
387 // boundaryPts,
388 // faces,
389 // owner,
390 // neighbour,
391 // patchTypes,
392 // patchNames,
393 // patchDicts,
394 // cellCentres,
395 // boundaryFacesToRemove
396 // );
397  }
398 }
399 
400 
401 Foam::autoPtr<Foam::fvMesh> Foam::conformalVoronoiMesh::createDummyMesh
402 (
403  const IOobject& io,
404  const wordList& patchNames,
405  const PtrList<dictionary>& patchDicts
406 ) const
407 {
408  autoPtr<fvMesh> meshPtr
409  (
410  new fvMesh
411  (
412  io,
413  pointField(),
414  faceList(),
415  cellList()
416  )
417  );
418  fvMesh& mesh = meshPtr();
419 
420  List<polyPatch*> patches(patchDicts.size());
421 
422  forAll(patches, patchi)
423  {
424  if
425  (
426  patchDicts.set(patchi)
427  && (
428  word(patchDicts[patchi].lookup("type"))
429  == processorPolyPatch::typeName
430  )
431  )
432  {
433  patches[patchi] = new processorPolyPatch
434  (
435  0, // patchSizes[p],
436  0, // patchStarts[p],
437  patchi,
438  mesh.boundaryMesh(),
439  readLabel(patchDicts[patchi].lookup("myProcNo")),
440  readLabel(patchDicts[patchi].lookup("neighbProcNo")),
442  );
443  }
444  else
445  {
446  patches[patchi] = polyPatch::New
447  (
448  patchDicts[patchi].lookup("type"),
449  patchNames[patchi],
450  0, // patchSizes[p],
451  0, // patchStarts[p],
452  patchi,
453  mesh.boundaryMesh()
454  ).ptr();
455  }
456  }
457 
458  mesh.addFvPatches(patches);
459 
460  return meshPtr;
461 }
462 
463 
464 void Foam::conformalVoronoiMesh::checkProcessorPatchesMatch
465 (
466  const PtrList<dictionary>& patchDicts
467 ) const
468 {
469  // Check patch sizes
470  labelListList procPatchSizes
471  (
472  Pstream::nProcs(),
474  );
475 
476  forAll(patchDicts, patchi)
477  {
478  if
479  (
480  patchDicts.set(patchi)
481  && (
482  word(patchDicts[patchi].lookup("type"))
483  == processorPolyPatch::typeName
484  )
485  )
486  {
487  const label procNeighb =
488  readLabel(patchDicts[patchi].lookup("neighbProcNo"));
489 
490  procPatchSizes[Pstream::myProcNo()][procNeighb]
491  = readLabel(patchDicts[patchi].lookup("nFaces"));
492  }
493  }
494 
495  Pstream::gatherList(procPatchSizes);
496 
497  if (Pstream::master())
498  {
499  bool allMatch = true;
500 
501  forAll(procPatchSizes, proci)
502  {
503  const labelList& patchSizes = procPatchSizes[proci];
504 
505  forAll(patchSizes, patchi)
506  {
507  if (patchSizes[patchi] != procPatchSizes[patchi][proci])
508  {
509  allMatch = false;
510 
511  Info<< indent << "Patches " << proci << " and " << patchi
512  << " have different sizes: " << patchSizes[patchi]
513  << " and " << procPatchSizes[patchi][proci] << endl;
514  }
515  }
516  }
517 
518  if (allMatch)
519  {
520  Info<< indent << "All processor patches have matching numbers of "
521  << "faces" << endl;
522  }
523  }
524 }
525 
526 
527 void Foam::conformalVoronoiMesh::reorderPoints
528 (
529  pointField& points,
530  labelList& boundaryPts,
531  faceList& faces,
532  const label nInternalFaces
533 ) const
534 {
535  Info<< incrIndent << indent << "Reordering points into internal/external"
536  << endl;
537 
538  labelList oldToNew(points.size(), label(0));
539 
540  // Find points that are internal
541  for (label fI = nInternalFaces; fI < faces.size(); ++fI)
542  {
543  const face& f = faces[fI];
544 
545  forAll(f, fpI)
546  {
547  oldToNew[f[fpI]] = 1;
548  }
549  }
550 
551  const label nInternalPoints = points.size() - sum(oldToNew);
552 
553  label countInternal = 0;
554  label countExternal = nInternalPoints;
555 
556  forAll(points, pI)
557  {
558  if (oldToNew[pI] == 0)
559  {
560  oldToNew[pI] = countInternal++;
561  }
562  else
563  {
564  oldToNew[pI] = countExternal++;
565  }
566  }
567 
568  Info<< indent
569  << "Number of internal points: " << countInternal << nl
570  << indent << "Number of external points: " << countExternal
571  << decrIndent << endl;
572 
573  inplaceReorder(oldToNew, points);
574  inplaceReorder(oldToNew, boundaryPts);
575 
576  forAll(faces, fI)
577  {
578  face& f = faces[fI];
579 
580  forAll(f, fpI)
581  {
582  f[fpI] = oldToNew[f[fpI]];
583  }
584  }
585 }
586 
587 
588 void Foam::conformalVoronoiMesh::reorderProcessorPatches
589 (
590  const word& meshName,
591  const fileName& instance,
592  const pointField& points,
593  faceList& faces,
594  labelList& owner,
595  labelList& neighbour,
596  const wordList& patchNames,
597  const PtrList<dictionary>& patchDicts
598 ) const
599 {
600  Info<< incrIndent << indent << "Reordering processor patches" << endl;
601 
602  Info<< incrIndent;
603  checkProcessorPatchesMatch(patchDicts);
604 
605  // Create dummy mesh with correct proc boundaries to do sorting
606  autoPtr<fvMesh> sortMeshPtr
607  (
608  createDummyMesh
609  (
610  IOobject
611  (
612  meshName,
613  instance,
614  runTime_,
615  IOobject::NO_READ,
617  false
618  ),
619  patchNames,
620  patchDicts
621  )
622  );
623  const fvMesh& sortMesh = sortMeshPtr();
624 
625  // Rotation on new faces.
626  labelList rotation(faces.size(), label(0));
627  labelList faceMap(faces.size(), label(-1));
628 
629  PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
630 
631  // Send ordering
632  forAll(sortMesh.boundaryMesh(), patchi)
633  {
634  const polyPatch& pp = sortMesh.boundaryMesh()[patchi];
635 
636  if (isA<processorPolyPatch>(pp))
637  {
638  refCast<const processorPolyPatch>(pp).initOrder
639  (
640  pBufs,
642  (
643  SubList<face>
644  (
645  faces,
646  readLabel(patchDicts[patchi].lookup("nFaces")),
647  readLabel(patchDicts[patchi].lookup("startFace"))
648  ),
649  points
650  )
651  );
652  }
653  }
654 
655  pBufs.finishedSends();
656 
657  Info<< incrIndent << indent << "Face ordering initialised..." << endl;
658 
659  // Receive and calculate ordering
660  bool anyChanged = false;
661 
662  forAll(sortMesh.boundaryMesh(), patchi)
663  {
664  const polyPatch& pp = sortMesh.boundaryMesh()[patchi];
665 
666  if (isA<processorPolyPatch>(pp))
667  {
668  const label nPatchFaces =
669  readLabel(patchDicts[patchi].lookup("nFaces"));
670  const label patchStartFace =
671  readLabel(patchDicts[patchi].lookup("startFace"));
672 
673  labelList patchFaceMap(nPatchFaces, label(-1));
674  labelList patchFaceRotation(nPatchFaces, label(0));
675 
676  bool changed = refCast<const processorPolyPatch>(pp).order
677  (
678  pBufs,
680  (
681  SubList<face>
682  (
683  faces,
684  nPatchFaces,
685  patchStartFace
686  ),
687  points
688  ),
689  patchFaceMap,
690  patchFaceRotation
691  );
692 
693  if (changed)
694  {
695  // Merge patch face reordering into mesh face reordering table
696  forAll(patchFaceRotation, patchFacei)
697  {
698  rotation[patchFacei + patchStartFace]
699  = patchFaceRotation[patchFacei];
700  }
701 
702  forAll(patchFaceMap, patchFacei)
703  {
704  if (patchFaceMap[patchFacei] != patchFacei)
705  {
706  faceMap[patchFacei + patchStartFace]
707  = patchFaceMap[patchFacei] + patchStartFace;
708  }
709  }
710 
711  anyChanged = true;
712  }
713  }
714  }
715 
716  Info<< incrIndent << indent << "Faces matched." << endl;
717 
718  reduce(anyChanged, orOp<bool>());
719 
720  if (anyChanged)
721  {
722  label nReorderedFaces = 0;
723 
724  forAll(faceMap, facei)
725  {
726  if (faceMap[facei] != -1)
727  {
728  nReorderedFaces++;
729  }
730  }
731 
732  if (nReorderedFaces > 0)
733  {
734  inplaceReorder(faceMap, faces);
735  inplaceReorder(faceMap, owner);
736  inplaceReorder(faceMap, neighbour);
737  }
738 
739  // Rotate faces (rotation is already in new face indices).
740  label nRotated = 0;
741 
742  forAll(rotation, facei)
743  {
744  if (rotation[facei] != 0)
745  {
746  faces[facei] = rotateList(faces[facei], rotation[facei]);
747  nRotated++;
748  }
749  }
750 
751  Info<< indent << returnReduce(nReorderedFaces, sumOp<label>())
752  << " faces have been reordered" << nl
753  << indent << returnReduce(nRotated, sumOp<label>())
754  << " faces have been rotated"
755  << decrIndent << decrIndent
756  << decrIndent << decrIndent << endl;
757  }
758 }
759 
760 
762 (
763  const word& meshName,
764  const fileName& instance,
765  pointField& points,
766  labelList& boundaryPts,
767  faceList& faces,
768  labelList& owner,
769  labelList& neighbour,
770  const wordList& patchNames,
771  const PtrList<dictionary>& patchDicts,
772  const pointField& cellCentres,
773  PackedBoolList& boundaryFacesToRemove
774 ) const
775 {
776  if (foamyHexMeshControls().objOutput())
777  {
779  (
780  time().path()/word(meshName + ".obj"),
781  points,
782  faces
783  );
784  }
785 
786  const label nInternalFaces = readLabel(patchDicts[0].lookup("startFace"));
787 
788  reorderPoints(points, boundaryPts, faces, nInternalFaces);
789 
790  if (Pstream::parRun())
791  {
792  reorderProcessorPatches
793  (
794  meshName,
795  instance,
796  points,
797  faces,
798  owner,
799  neighbour,
800  patchNames,
801  patchDicts
802  );
803  }
804 
805  Info<< incrIndent;
806  Info<< indent << "Constructing mesh" << endl;
807 
808  timeCheck("Before fvMesh construction");
809 
810  fvMesh mesh
811  (
812  IOobject
813  (
814  meshName,
815  instance,
816  runTime_,
817  IOobject::NO_READ,
819  ),
820  Foam::move(points),
821  Foam::move(faces),
822  Foam::move(owner),
823  Foam::move(neighbour)
824  );
825 
826  Info<< indent << "Adding patches to mesh" << endl;
827 
828  List<polyPatch*> patches(patchNames.size());
829 
830  label nValidPatches = 0;
831 
832  forAll(patches, p)
833  {
834  label totalPatchSize = readLabel(patchDicts[p].lookup("nFaces"));
835 
836  if
837  (
838  patchDicts.set(p)
839  && (
840  word(patchDicts[p].lookup("type"))
841  == processorPolyPatch::typeName
842  )
843  )
844  {
845  const_cast<dictionary&>(patchDicts[p]).set
846  (
847  "transform",
848  "coincidentFullMatch"
849  );
850 
851  // Do not create empty processor patches
852  if (totalPatchSize > 0)
853  {
854  patches[nValidPatches] = new processorPolyPatch
855  (
856  patchNames[p],
857  patchDicts[p],
858  nValidPatches,
859  mesh.boundaryMesh(),
860  processorPolyPatch::typeName
861  );
862 
863  nValidPatches++;
864  }
865  }
866  else
867  {
868  // Check that the patch is not empty on every processor
869  reduce(totalPatchSize, sumOp<label>());
870 
871  if (totalPatchSize > 0)
872  {
873  patches[nValidPatches] = polyPatch::New
874  (
875  patchNames[p],
876  patchDicts[p],
877  nValidPatches,
878  mesh.boundaryMesh()
879  ).ptr();
880 
881  nValidPatches++;
882  }
883  }
884  }
885 
886  patches.setSize(nValidPatches);
887 
888  mesh.addFvPatches(patches);
889 
890 
891 
892  // Add zones to the mesh
893  addZones(mesh, cellCentres);
894 
895 
896 
897  Info<< indent << "Add pointZones" << endl;
898  {
899  label sz = mesh.pointZones().size();
900 
901  DynamicList<label> bPts(boundaryPts.size());
902 
904  {
905  forAll(boundaryPts, ptI)
906  {
907  const label& bPtType = boundaryPts[ptI];
908 
909  if (bPtType == typeI)
910  {
911  bPts.append(ptI);
912  }
913  }
914 
915 // syncTools::syncPointList(mesh, bPts, maxEqOp<label>(), -1);
916 
917  Info<< incrIndent << indent
918  << "Adding " << bPts.size()
919  << " points of type " << dualMeshPointTypeNames_.words()[typeI]
920  << decrIndent << endl;
921 
922  mesh.pointZones().append
923  (
924  new pointZone
925  (
927  bPts,
928  sz + typeI,
929  mesh.pointZones()
930  )
931  );
932 
933  bPts.clear();
934  }
935  }
936 
937 
938 
939  // Add indirectPatchFaces to a face zone
940  Info<< indent << "Adding indirect patch faces set" << endl;
941 
943  (
944  mesh,
945  boundaryFacesToRemove,
946  orEqOp<unsigned int>()
947  );
948 
949  labelList addr(boundaryFacesToRemove.count());
950  label count = 0;
951 
952  forAll(boundaryFacesToRemove, facei)
953  {
954  if (boundaryFacesToRemove[facei])
955  {
956  addr[count++] = facei;
957  }
958  }
959 
960  addr.setSize(count);
961 
962  faceSet indirectPatchFaces
963  (
964  mesh,
965  "indirectPatchFaces",
966  addr,
968  );
969 
970  indirectPatchFaces.sync(mesh);
971 
972 
973  Info<< decrIndent;
974 
975  timeCheck("Before fvMesh filtering");
976 
977  autoPtr<polyMeshFilter> meshFilter;
978 
979  label nInitialBadFaces = 0;
980 
981  if (foamyHexMeshControls().filterEdges())
982  {
983  Info<< nl << "Filtering edges on polyMesh" << nl << endl;
984 
985  meshFilter.reset(new polyMeshFilter(mesh, boundaryPts));
986 
987  // Filter small edges only. This reduces the number of faces so that
988  // the face filtering is sped up.
989  nInitialBadFaces = meshFilter().filterEdges(0);
990  {
991  const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
992 
993  polyTopoChange meshMod(newMesh());
994 
995  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
996 
997  polyMeshFilter::copySets(newMesh(), mesh);
998  }
999  }
1000 
1001  if (foamyHexMeshControls().filterFaces())
1002  {
1003  labelIOList boundaryPtsIO
1004  (
1005  IOobject
1006  (
1007  "pointPriority",
1008  instance,
1009  time(),
1010  IOobject::NO_READ,
1012  ),
1013  labelList(mesh.nPoints(), labelMin)
1014  );
1015 
1016  forAll(mesh.points(), ptI)
1017  {
1018  boundaryPtsIO[ptI] = mesh.pointZones().whichZone(ptI);
1019  }
1020 
1021 
1022  Info<< nl << "Filtering faces on polyMesh" << nl << endl;
1023 
1024  meshFilter.reset(new polyMeshFilter(mesh, boundaryPtsIO));
1025 
1026  meshFilter().filter(nInitialBadFaces);
1027  {
1028  const autoPtr<fvMesh>& newMesh = meshFilter().filteredMesh();
1029 
1030  polyTopoChange meshMod(newMesh());
1031 
1032  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false);
1033 
1034  polyMeshFilter::copySets(newMesh(), mesh);
1035  }
1036  }
1037 
1038  timeCheck("After fvMesh filtering");
1039 
1040  mesh.setInstance(instance);
1041 
1042  if (!mesh.write())
1043  {
1045  << "Failed writing polyMesh."
1046  << exit(FatalError);
1047  }
1048  else
1049  {
1050  Info<< nl << "Written filtered mesh to "
1051  << mesh.polyMesh::instance() << nl
1052  << endl;
1053  }
1054 
1055  {
1056  pointScalarField boundaryPtsScalarField
1057  (
1058  IOobject
1059  (
1060  "boundaryPoints_collapsed",
1061  instance,
1062  time(),
1063  IOobject::NO_READ,
1065  ),
1066  pointMesh::New(mesh),
1067  scalar(labelMin)
1068  );
1069 
1070  labelIOList boundaryPtsIO
1071  (
1072  IOobject
1073  (
1074  "pointPriority",
1075  instance,
1076  time(),
1077  IOobject::NO_READ,
1079  ),
1080  labelList(mesh.nPoints(), labelMin)
1081  );
1082 
1083  forAll(mesh.points(), ptI)
1084  {
1085  boundaryPtsScalarField[ptI] = mesh.pointZones().whichZone(ptI);
1086  boundaryPtsIO[ptI] = mesh.pointZones().whichZone(ptI);
1087  }
1088 
1089  boundaryPtsScalarField.write();
1090  boundaryPtsIO.write();
1091  }
1092 
1093 // writeCellSizes(mesh);
1094 
1095 // writeCellAlignments(mesh);
1096 
1097 // writeCellCentres(mesh);
1098 
1100 }
1101 
1102 
1104 (
1105  const fvMesh& mesh
1106 ) const
1107 {
1108  {
1109  timeCheck("Start writeCellSizes");
1110 
1111  Info<< nl << "Create targetCellSize volScalarField" << endl;
1112 
1113  volScalarField targetCellSize
1114  (
1115  IOobject
1116  (
1117  "targetCellSize",
1118  mesh.polyMesh::instance(),
1119  mesh,
1122  ),
1123  mesh,
1125  zeroGradientFvPatchScalarField::typeName
1126  );
1127 
1128  scalarField& cellSize = targetCellSize.primitiveFieldRef();
1129 
1130  const vectorField& C = mesh.cellCentres();
1131 
1132  forAll(cellSize, i)
1133  {
1134  cellSize[i] = cellShapeControls().cellSize(C[i]);
1135  }
1136 
1137  // Info<< nl << "Create targetCellVolume volScalarField" << endl;
1138 
1139  // volScalarField targetCellVolume
1140  // (
1141  // IOobject
1142  // (
1143  // "targetCellVolume",
1144  // mesh.polyMesh::instance(),
1145  // mesh,
1146  // IOobject::NO_READ,
1147  // IOobject::AUTO_WRITE
1148  // ),
1149  // mesh,
1150  // dimensionedScalar(dimLength, 0),
1151  // zeroGradientFvPatchScalarField::typeName
1152  // );
1153 
1154  // targetCellVolume.primitiveFieldRef() = pow3(cellSize);
1155 
1156  // Info<< nl << "Create actualCellVolume volScalarField" << endl;
1157 
1158  // volScalarField actualCellVolume
1159  // (
1160  // IOobject
1161  // (
1162  // "actualCellVolume",
1163  // mesh.polyMesh::instance(),
1164  // mesh,
1165  // IOobject::NO_READ,
1166  // IOobject::AUTO_WRITE
1167  // ),
1168  // mesh,
1169  // dimensionedScalar(dimVolume, 0),
1170  // zeroGradientFvPatchScalarField::typeName
1171  // );
1172 
1173  // actualCellVolume.primitiveFieldRef() = mesh.cellVolumes();
1174 
1175  // Info<< nl << "Create equivalentCellSize volScalarField" << endl;
1176 
1177  // volScalarField equivalentCellSize
1178  // (
1179  // IOobject
1180  // (
1181  // "equivalentCellSize",
1182  // mesh.polyMesh::instance(),
1183  // mesh,
1184  // IOobject::NO_READ,
1185  // IOobject::AUTO_WRITE
1186  // ),
1187  // mesh,
1188  // dimensionedScalar(dimLength, 0),
1189  // zeroGradientFvPatchScalarField::typeName
1190  // );
1191 
1192  // equivalentCellSize.primitiveFieldRef() = pow
1193  // (
1194  // actualCellVolume.primitiveField(),
1195  // 1.0/3.0
1196  // );
1197 
1198  targetCellSize.correctBoundaryConditions();
1199  // targetCellVolume.correctBoundaryConditions();
1200  // actualCellVolume.correctBoundaryConditions();
1201  // equivalentCellSize.correctBoundaryConditions();
1202 
1203  targetCellSize.write();
1204  // targetCellVolume.write();
1205  // actualCellVolume.write();
1206  // equivalentCellSize.write();
1207  }
1208 
1209  // {
1210  // polyMesh tetMesh
1211  // (
1212  // IOobject
1213  // (
1214  // "tetDualMesh",
1215  // runTime_.constant(),
1216  // runTime_,
1217  // IOobject::MUST_READ
1218  // )
1219  // );
1220 
1221  // pointMesh ptMesh(tetMesh);
1222 
1223  // pointScalarField ptTargetCellSize
1224  // (
1225  // IOobject
1226  // (
1227  // "ptTargetCellSize",
1228  // runTime_.timeName(),
1229  // tetMesh,
1230  // IOobject::NO_READ,
1231  // IOobject::AUTO_WRITE
1232  // ),
1233  // ptMesh,
1234  // dimensionedScalar(dimLength, 0),
1235  // pointPatchVectorField::calculatedType()
1236  // );
1237 
1238  // scalarField& cellSize = ptTargetCellSize.primitiveFieldRef();
1239 
1240  // const vectorField& P = tetMesh.points();
1241 
1242  // forAll(cellSize, i)
1243  // {
1244  // cellSize[i] = cellShapeControls().cellSize(P[i]);
1245  // }
1246 
1247  // ptTargetCellSize.write();
1248  // }
1249 }
1250 
1251 
1253 (
1254  const fvMesh& mesh
1255 ) const
1256 {
1257 // Info<< nl << "Create cellAlignments volTensorField" << endl;
1258 //
1259 // volTensorField cellAlignments
1260 // (
1261 // IOobject
1262 // (
1263 // "cellAlignments",
1264 // mesh.polyMesh::instance(),
1265 // mesh,
1266 // IOobject::NO_READ,
1267 // IOobject::AUTO_WRITE
1268 // ),
1269 // mesh,
1270 // tensor::I,
1271 // zeroGradientFvPatchTensorField::typeName
1272 // );
1273 //
1274 // tensorField& cellAlignment = cellAlignments.primitiveFieldRef();
1275 //
1276 // const vectorField& C = mesh.cellCentres();
1277 //
1278 // vectorField xDir(cellAlignment.size());
1279 // vectorField yDir(cellAlignment.size());
1280 // vectorField zDir(cellAlignment.size());
1281 //
1282 // forAll(cellAlignment, i)
1283 // {
1284 // cellAlignment[i] = cellShapeControls().cellAlignment(C[i]);
1285 // xDir[i] = cellAlignment[i] & vector(1, 0, 0);
1286 // yDir[i] = cellAlignment[i] & vector(0, 1, 0);
1287 // zDir[i] = cellAlignment[i] & vector(0, 0, 1);
1288 // }
1289 //
1290 // OFstream xStr("xDir.obj");
1291 // OFstream yStr("yDir.obj");
1292 // OFstream zStr("zDir.obj");
1293 //
1294 // forAll(xDir, i)
1295 // {
1296 // meshTools::writeOBJ(xStr, C[i], C[i] + xDir[i]);
1297 // meshTools::writeOBJ(yStr, C[i], C[i] + yDir[i]);
1298 // meshTools::writeOBJ(zStr, C[i], C[i] + zDir[i]);
1299 // }
1300 //
1301 // cellAlignments.correctBoundaryConditions();
1302 //
1303 // cellAlignments.write();
1304 }
1305 
1306 
1308 (
1309  const fvMesh& mesh
1310 ) const
1311 {
1312  Info<< "Writing components of cellCentre positions to volScalarFields"
1313  << " ccx, ccy, ccz in " << runTime_.timeName() << endl;
1314 
1315  for (direction i=0; i<vector::nComponents; i++)
1316  {
1317  volScalarField cci
1318  (
1319  IOobject
1320  (
1321  "cc" + word(vector::componentNames[i]),
1322  runTime_.timeName(),
1323  mesh,
1326  ),
1327  mesh.C().component(i)
1328  );
1329 
1330  cci.write();
1331  }
1332 }
1333 
1334 
1336 (
1337  const polyMesh& mesh
1338 ) const
1339 {
1340  timeCheck("Start findRemainingProtrusionSet");
1341 
1342  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1343 
1344  labelHashSet protrudingBoundaryPoints;
1345 
1346  forAll(patches, patchi)
1347  {
1348  const polyPatch& patch = patches[patchi];
1349 
1350  forAll(patch.localPoints(), pLPI)
1351  {
1352  label meshPtI = patch.meshPoints()[pLPI];
1353 
1354  const Foam::point& pt = patch.localPoints()[pLPI];
1355 
1356  if
1357  (
1358  geometryToConformTo_.wellOutside
1359  (
1360  pt,
1361  sqr(targetCellSize(pt))
1362  )
1363  )
1364  {
1365  protrudingBoundaryPoints.insert(meshPtI);
1366  }
1367  }
1368  }
1369 
1370  cellSet protrudingCells
1371  (
1372  mesh,
1373  "foamyHexMesh_remainingProtrusions",
1374  mesh.nCells()/1000
1375  );
1376 
1377  forAllConstIter(labelHashSet, protrudingBoundaryPoints, iter)
1378  {
1379  const label pointi = iter.key();
1380  const labelList& pCells = mesh.pointCells()[pointi];
1381 
1382  forAll(pCells, pCI)
1383  {
1384  protrudingCells.insert(pCells[pCI]);
1385  }
1386  }
1387 
1388  label protrudingCellsSize = protrudingCells.size();
1389 
1390  reduce(protrudingCellsSize, sumOp<label>());
1391 
1392  if (foamyHexMeshControls().objOutput() && protrudingCellsSize > 0)
1393  {
1394  Info<< nl << "Found " << protrudingCellsSize
1395  << " cells protruding from the surface, writing cellSet "
1396  << protrudingCells.name()
1397  << endl;
1398 
1399  protrudingCells.write();
1400  }
1401 
1402  return Foam::move(protrudingCells);
1403 }
1404 
1405 
1406 void Foam::conformalVoronoiMesh::writePointPairs
1407 (
1408  const fileName& fName
1409 ) const
1410 {
1411  OBJstream os(fName);
1412 
1413  for
1414  (
1415  Delaunay::Finite_edges_iterator eit = finite_edges_begin();
1416  eit != finite_edges_end();
1417  ++eit
1418  )
1419  {
1420  Cell_handle c = eit->first;
1421  Vertex_handle vA = c->vertex(eit->second);
1422  Vertex_handle vB = c->vertex(eit->third);
1423 
1424  if (ptPairs_.isPointPair(vA, vB))
1425  {
1426  os.write
1427  (
1428  linePointRef(topoint(vA->point()), topoint(vB->point()))
1429  );
1430  }
1431  }
1432 }
1433 
1434 
1435 // ************************************************************************* //
GeometricField< scalar, pointPatchField, pointMesh > pointScalarField
Definition: pointFields.H:49
labelHashSet findRemainingProtrusionSet(const polyMesh &mesh) const
Find the cellSet of the boundary cells which have points that.
static wordList words()
The set of names as a list of words.
Definition: NamedEnum.C:105
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
Delaunay::Cell_handle Cell_handle
void inplaceReorder(const labelUList &oldToNew, ListType &)
Inplace reorder the elements of a list.
#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
Ostream & indent(Ostream &os)
Indent stream.
Definition: Ostream.H:226
Delaunay::Vertex_handle Vertex_handle
pointFromPoint topoint(const Point &P)
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
error FatalError
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
ListType rotateList(const ListType &list, const label n)
Rotate a list by n places. If n is positive rotate clockwise/right/down.
dimensionedSymmTensor sqr(const dimensionedVector &dv)
static fvMesh * meshPtr
Definition: globalFoam.H:52
const cellShapeControl & cellShapeControls() const
Return the cellShapeControl object.
uint8_t direction
Definition: direction.H:45
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:429
List< face > faceList
Definition: faceListFwd.H:43
HashTable< label, labelPair, FixedList< label, 2 >::Hash<> > labelTolabelPairHashTable
Definition: DelaunayMesh.H:90
static word defaultRegion
Return the default region name.
Definition: polyMesh.H:309
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
static bool master(const label communicator=0)
Am I the master process.
Definition: UPstream.H:423
static void timeCheck(const Time &runTime, const string &description=string::null, const bool check=true)
Write the elapsedCpuTime and memory usage, with an optional.
void writeCellCentres(const fvMesh &mesh) const
Calculate and write the cell centres.
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
patches[0]
void writeMesh(const fileName &instance)
Prepare data and call writeMesh for polyMesh and.
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
dimensioned< Type > sum(const DimensionedField< Type, GeoMesh > &df)
Various functions to operate on Lists.
static const pointMesh & New(const polyMesh &mesh)
Definition: MeshObject.C:44
Useful combination of include files which define Sin, Sout and Serr and the use of IO streams general...
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:208
void writeObjMesh(const fileName &fName, const pointField &points, const faceList &faces)
Write an OBJ mesh consisting of points and faces.
static const direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:99
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
stressControl lookup("compactNormalStress") >> compactNormalStress
void writeCellSizes(const fvMesh &mesh) const
Calculate and write a field of the target cell size,.
dynamicFvMesh & mesh
line< point, const point & > linePointRef
Line using referred points.
Definition: linePointRef.H:45
const pointField & points
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
wordList patchNames(nPatches)
static const word null
An empty word.
Definition: word.H:77
List< label > labelList
A List of labels.
Definition: labelList.H:56
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
label readLabel(Istream &is)
Definition: label.H:64
static const char nl
Definition: Ostream.H:265
Delaunay ::Finite_cells_iterator Finite_cells_iterator
Field< bool > wellOutside(const pointField &samplePts, const scalarField &testDistSqr) const
Check if point is outside surfaces to conform to by at least.
Ostream & decrIndent(Ostream &os)
Decrement the indent level.
Definition: Ostream.H:240
labelList f(nPoints)
PrimitivePatch< SubList< face >, const pointField & > primitivePatch
Addressing for a faceList slice.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
bool isPointPair(const Vertex_handle &vA, const Vertex_handle &vB) const
List< word > wordList
A List of words.
Definition: fileName.H:54
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:381
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:399
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
label patchi
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
messageStream Info
void writeCellAlignments(const fvMesh &mesh) const
static const NamedEnum< dualMeshPointType, 5 > dualMeshPointTypeNames_
cellShapeControlMesh & shapeControlMesh()
Field< vector > vectorField
Specialisation of Field<T> for vector.
static autoPtr< polyPatch > New(const word &patchType, const word &name, const label size, const label start, const label index, const polyBoundaryMesh &bm)
Return a pointer to a new patch created on freestore from.
Definition: polyPatchNew.C:32
virtual Ostream & write(const token &)=0
Write next token to stream.
const cvControls & foamyHexMeshControls() const
Return the foamyHexMeshControls object.
const Time & time() const
Return the Time object.
PtrList< dictionary > patchDicts
Definition: readKivaGrid.H:537
virtual bool write(const bool write=true) const
Write using setting from DB.
volScalarField & p
virtual void print(Ostream &) const
Print description of IOstream to Ostream.
Definition: IOstream.C:126
void writeInternalDelaunayVertices(const fileName &instance, const Triangulation &t)
Write the internal Delaunay vertices of the tessellation as a.
Ostream & incrIndent(Ostream &os)
Increment the indent level.
Definition: Ostream.H:233
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
List< cell > cellList
list of cells
Definition: cellList.H:42
static const label labelMin
Definition: label.H:61
IOList< label > labelIOList
Label container classes.
Definition: labelIOList.H:42
autoPtr< polyMesh > createMesh(const fileName &name, labelTolabelPairHashTable &vertexMap, labelList &cellMap, const bool writeDelaunayData=true) const
Create an fvMesh from the triangulation.
scalar cellSize(const point &pt) const
Return the cell size at the given location.
fileName path(UMean.rootPath()/UMean.caseName()/functionObjects::writeFile::outputPrefix/"graphs"/UMean.instance())