meshRefinement.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2018 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 "meshRefinement.H"
27 #include "volMesh.H"
28 #include "volFields.H"
29 #include "surfaceMesh.H"
30 #include "syncTools.H"
31 #include "Time.H"
32 #include "refinementHistory.H"
33 #include "refinementSurfaces.H"
34 #include "refinementFeatures.H"
35 #include "decompositionMethod.H"
36 #include "regionSplit.H"
37 #include "fvMeshDistribute.H"
38 #include "indirectPrimitivePatch.H"
39 #include "polyTopoChange.H"
40 #include "removeCells.H"
41 #include "mapDistributePolyMesh.H"
42 #include "localPointRegion.H"
43 #include "pointMesh.H"
44 #include "pointFields.H"
45 #include "slipPointPatchFields.H"
49 #include "processorPointPatch.H"
50 #include "globalIndex.H"
51 #include "meshTools.H"
52 #include "OFstream.H"
53 #include "geomDecomp.H"
54 #include "Random.H"
55 #include "searchableSurfaces.H"
56 #include "treeBoundBox.H"
57 #include "fvMeshTools.H"
58 
59 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
60 
61 namespace Foam
62 {
63  defineTypeNameAndDebug(meshRefinement, 0);
64 
65  template<>
66  const char* Foam::NamedEnum
67  <
69  5
70  >::names[] =
71  {
72  "mesh",
73  "intersections",
74  "featureSeeds",
75  "attraction",
76  "layerInfo"
77  };
78 
79  template<>
80  const char* Foam::NamedEnum
81  <
83  1
84  >::names[] =
85  {
86  "layerInfo"
87  };
88 
89  template<>
90  const char* Foam::NamedEnum
91  <
93  5
94  >::names[] =
95  {
96  "mesh",
97  "noRefinement",
98  "scalarLevels",
99  "layerSets",
100  "layerFields"
101  };
102 
103 }
104 
107 
110 
113 
114 
115 Foam::meshRefinement::writeType Foam::meshRefinement::writeLevel_;
116 
117 Foam::meshRefinement::outputType Foam::meshRefinement::outputLevel_;
118 
119 
120 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
121 
122 void Foam::meshRefinement::calcNeighbourData
123 (
124  labelList& neiLevel,
125  pointField& neiCc
126 ) const
127 {
128  const labelList& cellLevel = meshCutter_.cellLevel();
129  const pointField& cellCentres = mesh_.cellCentres();
130 
131  label nBoundaryFaces = mesh_.nFaces() - mesh_.nInternalFaces();
132 
133  if (neiLevel.size() != nBoundaryFaces || neiCc.size() != nBoundaryFaces)
134  {
136  << nBoundaryFaces << " neiLevel:" << neiLevel.size()
137  << abort(FatalError);
138  }
139 
140  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
141 
142  labelHashSet addedPatchIDSet(meshedPatches());
143 
144  forAll(patches, patchi)
145  {
146  const polyPatch& pp = patches[patchi];
147 
148  const labelUList& faceCells = pp.faceCells();
149  const vectorField::subField faceCentres = pp.faceCentres();
150  const vectorField::subField faceAreas = pp.faceAreas();
151 
152  label bFacei = pp.start()-mesh_.nInternalFaces();
153 
154  if (pp.coupled())
155  {
156  forAll(faceCells, i)
157  {
158  neiLevel[bFacei] = cellLevel[faceCells[i]];
159  neiCc[bFacei] = cellCentres[faceCells[i]];
160  bFacei++;
161  }
162  }
163  else if (addedPatchIDSet.found(patchi))
164  {
165  // Face was introduced from cell-cell intersection. Try to
166  // reconstruct other side cell(centre). Three possibilities:
167  // - cells same size.
168  // - preserved cell smaller. Not handled.
169  // - preserved cell larger.
170  forAll(faceCells, i)
171  {
172  // Extrapolate the face centre.
173  vector fn = faceAreas[i];
174  fn /= mag(fn)+vSmall;
175 
176  label own = faceCells[i];
177  label ownLevel = cellLevel[own];
178  label faceLevel = meshCutter_.faceLevel(pp.start()+i);
179 
180  // Normal distance from face centre to cell centre
181  scalar d = ((faceCentres[i] - cellCentres[own]) & fn);
182  if (faceLevel > ownLevel)
183  {
184  // Other cell more refined. Adjust normal distance
185  d *= 0.5;
186  }
187  neiLevel[bFacei] = faceLevel;
188  // Calculate other cell centre by extrapolation
189  neiCc[bFacei] = faceCentres[i] + d*fn;
190  bFacei++;
191  }
192  }
193  else
194  {
195  forAll(faceCells, i)
196  {
197  neiLevel[bFacei] = cellLevel[faceCells[i]];
198  neiCc[bFacei] = faceCentres[i];
199  bFacei++;
200  }
201  }
202  }
203 
204  // Swap coupled boundaries. Apply separation to cc since is coordinate.
206  syncTools::swapBoundaryFaceList(mesh_, neiLevel);
207 }
208 
209 
210 void Foam::meshRefinement::updateIntersections(const labelList& changedFaces)
211 {
212  const pointField& cellCentres = mesh_.cellCentres();
213 
214  // Stats on edges to test. Count proc faces only once.
215  PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
216 
217  {
218  label nMasterFaces = 0;
219  forAll(isMasterFace, facei)
220  {
221  if (isMasterFace.get(facei) == 1)
222  {
223  nMasterFaces++;
224  }
225  }
226  reduce(nMasterFaces, sumOp<label>());
227 
228  label nChangedFaces = 0;
229  forAll(changedFaces, i)
230  {
231  if (isMasterFace.get(changedFaces[i]) == 1)
232  {
233  nChangedFaces++;
234  }
235  }
236  reduce(nChangedFaces, sumOp<label>());
237 
238  Info<< "Edge intersection testing:" << nl
239  << " Number of edges : " << nMasterFaces << nl
240  << " Number of edges to retest : " << nChangedFaces
241  << endl;
242  }
243 
244 
245  // Get boundary face centre and level. Coupled aware.
246  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
247  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
248  calcNeighbourData(neiLevel, neiCc);
249 
250  // Collect segments we want to test for
251  pointField start(changedFaces.size());
252  pointField end(changedFaces.size());
253 
254  forAll(changedFaces, i)
255  {
256  label facei = changedFaces[i];
257  label own = mesh_.faceOwner()[facei];
258 
259  start[i] = cellCentres[own];
260  if (mesh_.isInternalFace(facei))
261  {
262  end[i] = cellCentres[mesh_.faceNeighbour()[facei]];
263  }
264  else
265  {
266  end[i] = neiCc[facei-mesh_.nInternalFaces()];
267  }
268  }
269 
270  // Extend segments a bit
271  {
272  const vectorField smallVec(rootSmall*(end-start));
273  start -= smallVec;
274  end += smallVec;
275  }
276 
277 
278  // Do tests in one go
279  labelList surfaceHit;
280  {
281  labelList surfaceLevel;
282  surfaces_.findHigherIntersection
283  (
284  start,
285  end,
286  labelList(start.size(), -1), // accept any intersection
287  surfaceHit,
288  surfaceLevel
289  );
290  }
291 
292  // Keep just surface hit
293  forAll(surfaceHit, i)
294  {
295  surfaceIndex_[changedFaces[i]] = surfaceHit[i];
296  }
297 
298  // Make sure both sides have same information. This should be
299  // case in general since same vectors but just to make sure.
300  syncTools::syncFaceList(mesh_, surfaceIndex_, maxEqOp<label>());
301 
302  label nHits = countHits();
303  label nTotHits = returnReduce(nHits, sumOp<label>());
304 
305  Info<< " Number of intersected edges : " << nTotHits << endl;
306 
307  // Set files to same time as mesh
308  setInstance(mesh_.facesInstance());
309 }
310 
311 
313 (
314  const string& msg,
315  const polyMesh& mesh,
316  const List<scalar>& fld
317 )
318 {
319  if (fld.size() != mesh.nPoints())
320  {
322  << msg << endl
323  << "fld size:" << fld.size() << " mesh points:" << mesh.nPoints()
324  << abort(FatalError);
325  }
326 
327  Pout<< "Checking field " << msg << endl;
328  scalarField minFld(fld);
330  (
331  mesh,
332  minFld,
333  minEqOp<scalar>(),
334  great
335  );
336  scalarField maxFld(fld);
338  (
339  mesh,
340  maxFld,
341  maxEqOp<scalar>(),
342  -great
343  );
344  forAll(minFld, pointi)
345  {
346  const scalar& minVal = minFld[pointi];
347  const scalar& maxVal = maxFld[pointi];
348  if (mag(minVal-maxVal) > small)
349  {
350  Pout<< msg << " at:" << mesh.points()[pointi] << nl
351  << " minFld:" << minVal << nl
352  << " maxFld:" << maxVal << nl
353  << endl;
354  }
355  }
356 }
357 
358 
360 (
361  const string& msg,
362  const polyMesh& mesh,
363  const List<point>& fld
364 )
365 {
366  if (fld.size() != mesh.nPoints())
367  {
369  << msg << endl
370  << "fld size:" << fld.size() << " mesh points:" << mesh.nPoints()
371  << abort(FatalError);
372  }
373 
374  Pout<< "Checking field " << msg << endl;
375  pointField minFld(fld);
377  (
378  mesh,
379  minFld,
381  point(great, great, great)
382  );
383  pointField maxFld(fld);
385  (
386  mesh,
387  maxFld,
390  );
391  forAll(minFld, pointi)
392  {
393  const point& minVal = minFld[pointi];
394  const point& maxVal = maxFld[pointi];
395  if (mag(minVal-maxVal) > small)
396  {
397  Pout<< msg << " at:" << mesh.points()[pointi] << nl
398  << " minFld:" << minVal << nl
399  << " maxFld:" << maxVal << nl
400  << endl;
401  }
402  }
403 }
404 
405 
407 {
408  Pout<< "meshRefinement::checkData() : Checking refinement structure."
409  << endl;
410  meshCutter_.checkMesh();
411 
412  Pout<< "meshRefinement::checkData() : Checking refinement levels."
413  << endl;
414  meshCutter_.checkRefinementLevels(1, labelList(0));
415 
416 
417  label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
418 
419  Pout<< "meshRefinement::checkData() : Checking synchronization."
420  << endl;
421 
422  // Check face centres
423  {
424  // Boundary face centres
425  pointField::subList boundaryFc
426  (
427  mesh_.faceCentres(),
428  nBnd,
429  mesh_.nInternalFaces()
430  );
431 
432  // Get neighbouring face centres
433  pointField neiBoundaryFc(boundaryFc);
435  (
436  mesh_,
437  neiBoundaryFc,
438  eqOp<point>()
439  );
440 
441  // Compare
442  testSyncBoundaryFaceList
443  (
444  mergeDistance_,
445  "testing faceCentres : ",
446  boundaryFc,
447  neiBoundaryFc
448  );
449  }
450  // Check meshRefinement
451  {
452  // Get boundary face centre and level. Coupled aware.
453  labelList neiLevel(nBnd);
454  pointField neiCc(nBnd);
455  calcNeighbourData(neiLevel, neiCc);
456 
457  // Collect segments we want to test for
458  pointField start(mesh_.nFaces());
459  pointField end(mesh_.nFaces());
460 
461  forAll(start, facei)
462  {
463  start[facei] = mesh_.cellCentres()[mesh_.faceOwner()[facei]];
464 
465  if (mesh_.isInternalFace(facei))
466  {
467  end[facei] = mesh_.cellCentres()[mesh_.faceNeighbour()[facei]];
468  }
469  else
470  {
471  end[facei] = neiCc[facei-mesh_.nInternalFaces()];
472  }
473  }
474 
475  // Extend segments a bit
476  {
477  const vectorField smallVec(rootSmall*(end-start));
478  start -= smallVec;
479  end += smallVec;
480  }
481 
482 
483  // Do tests in one go
484  labelList surfaceHit;
485  {
486  labelList surfaceLevel;
487  surfaces_.findHigherIntersection
488  (
489  start,
490  end,
491  labelList(start.size(), -1), // accept any intersection
492  surfaceHit,
493  surfaceLevel
494  );
495  }
496  // Get the coupled hit
497  labelList neiHit
498  (
500  (
501  surfaceHit,
502  nBnd,
503  mesh_.nInternalFaces()
504  )
505  );
506  syncTools::swapBoundaryFaceList(mesh_, neiHit);
507 
508  // Check
509  forAll(surfaceHit, facei)
510  {
511  if (surfaceIndex_[facei] != surfaceHit[facei])
512  {
513  if (mesh_.isInternalFace(facei))
514  {
516  << "Internal face:" << facei
517  << " fc:" << mesh_.faceCentres()[facei]
518  << " cached surfaceIndex_:" << surfaceIndex_[facei]
519  << " current:" << surfaceHit[facei]
520  << " ownCc:"
521  << mesh_.cellCentres()[mesh_.faceOwner()[facei]]
522  << " neiCc:"
523  << mesh_.cellCentres()[mesh_.faceNeighbour()[facei]]
524  << endl;
525  }
526  else if
527  (
528  surfaceIndex_[facei]
529  != neiHit[facei-mesh_.nInternalFaces()]
530  )
531  {
533  << "Boundary face:" << facei
534  << " fc:" << mesh_.faceCentres()[facei]
535  << " cached surfaceIndex_:" << surfaceIndex_[facei]
536  << " current:" << surfaceHit[facei]
537  << " ownCc:"
538  << mesh_.cellCentres()[mesh_.faceOwner()[facei]]
539  << " end:" << end[facei]
540  << endl;
541  }
542  }
543  }
544  }
545  {
546  labelList::subList boundarySurface
547  (
548  surfaceIndex_,
549  mesh_.nFaces()-mesh_.nInternalFaces(),
550  mesh_.nInternalFaces()
551  );
552 
553  labelList neiBoundarySurface(boundarySurface);
555  (
556  mesh_,
557  neiBoundarySurface
558  );
559 
560  // Compare
561  testSyncBoundaryFaceList
562  (
563  0, // tolerance
564  "testing surfaceIndex() : ",
565  boundarySurface,
566  neiBoundarySurface
567  );
568  }
569 
570 
571  // Find duplicate faces
572  Pout<< "meshRefinement::checkData() : Counting duplicate faces."
573  << endl;
574 
575  labelList duplicateFace
576  (
578  (
579  mesh_,
580  identity(mesh_.nFaces()-mesh_.nInternalFaces())
581  + mesh_.nInternalFaces()
582  )
583  );
584 
585  // Count
586  {
587  label nDup = 0;
588 
589  forAll(duplicateFace, i)
590  {
591  if (duplicateFace[i] != -1)
592  {
593  nDup++;
594  }
595  }
596  nDup /= 2; // will have counted both faces of duplicate
597  Pout<< "meshRefinement::checkData() : Found " << nDup
598  << " duplicate pairs of faces." << endl;
599  }
600 }
601 
602 
604 {
605  meshCutter_.setInstance(inst);
606  surfaceIndex_.instance() = inst;
607 }
608 
609 
610 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::doRemoveCells
611 (
612  const labelList& cellsToRemove,
613  const labelList& exposedFaces,
614  const labelList& exposedPatchIDs,
615  removeCells& cellRemover
616 )
617 {
618  polyTopoChange meshMod(mesh_);
619 
620  // Arbitrary: put exposed faces into last patch.
621  cellRemover.setRefinement
622  (
623  cellsToRemove,
624  exposedFaces,
625  exposedPatchIDs,
626  meshMod
627  );
628 
629  // Change the mesh (no inflation)
630  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
631 
632  // Update fields
633  mesh_.updateMesh(map);
634 
635  // Move mesh (since morphing might not do this)
636  if (map().hasMotionPoints())
637  {
638  mesh_.movePoints(map().preMotionPoints());
639  }
640  else
641  {
642  // Delete mesh volumes. No other way to do this?
643  mesh_.clearOut();
644  }
645 
646  // Reset the instance for if in overwrite mode
647  mesh_.setInstance(timeName());
648  setInstance(mesh_.facesInstance());
649 
650  // Update local mesh data
651  cellRemover.updateMesh(map);
652 
653  // Update intersections. Recalculate intersections for exposed faces.
654  labelList newExposedFaces = renumber
655  (
656  map().reverseFaceMap(),
657  exposedFaces
658  );
659 
660  // Pout<< "removeCells : updating intersections for "
661  // << newExposedFaces.size() << " newly exposed faces." << endl;
662 
663  updateMesh(map, newExposedFaces);
664 
665  return map;
666 }
667 
668 
670 (
671  const labelList& splitFaces,
672  const labelPairList& splits
673 )
674 {
675  polyTopoChange meshMod(mesh_);
676 
677  forAll(splitFaces, i)
678  {
679  label facei = splitFaces[i];
680  const face& f = mesh_.faces()[facei];
681 
682  // Split as start and end index in face
683  const labelPair& split = splits[i];
684 
685  label nVerts = split[1]-split[0];
686  if (nVerts < 0)
687  {
688  nVerts += f.size();
689  }
690  nVerts += 1;
691 
692 
693  // Split into f0, f1
694  face f0(nVerts);
695 
696  label fp = split[0];
697  forAll(f0, i)
698  {
699  f0[i] = f[fp];
700  fp = f.fcIndex(fp);
701  }
702 
703  face f1(f.size()-f0.size()+2);
704  fp = split[1];
705  forAll(f1, i)
706  {
707  f1[i] = f[fp];
708  fp = f.fcIndex(fp);
709  }
710 
711 
712  // Determine face properties
713  label own = mesh_.faceOwner()[facei];
714  label nei = -1;
715  label patchi = -1;
716  if (facei >= mesh_.nInternalFaces())
717  {
718  patchi = mesh_.boundaryMesh().whichPatch(facei);
719  }
720  else
721  {
722  nei = mesh_.faceNeighbour()[facei];
723  }
724 
725  label zoneI = mesh_.faceZones().whichZone(facei);
726  bool zoneFlip = false;
727  if (zoneI != -1)
728  {
729  const faceZone& fz = mesh_.faceZones()[zoneI];
730  zoneFlip = fz.flipMap()[fz.whichFace(facei)];
731  }
732 
733 
734  if (debug)
735  {
736  Pout<< "face:" << facei << " verts:" << f
737  << " split into f0:" << f0
738  << " f1:" << f1 << endl;
739  }
740 
741  // Change/add faces
742  meshMod.modifyFace
743  (
744  f0, // modified face
745  facei, // label of face
746  own, // owner
747  nei, // neighbour
748  false, // face flip
749  patchi, // patch for face
750  zoneI, // zone for face
751  zoneFlip // face flip in zone
752  );
753 
754  meshMod.addFace
755  (
756  f1, // modified face
757  own, // owner
758  nei, // neighbour
759  -1, // master point
760  -1, // master edge
761  facei, // master face
762  false, // face flip
763  patchi, // patch for face
764  zoneI, // zone for face
765  zoneFlip // face flip in zone
766  );
767  }
768 
769 
770 
771  // Change the mesh (no inflation)
772  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false, true);
773 
774  // Update fields
775  mesh_.updateMesh(map);
776 
777  // Move mesh (since morphing might not do this)
778  if (map().hasMotionPoints())
779  {
780  mesh_.movePoints(map().preMotionPoints());
781  }
782  else
783  {
784  // Delete mesh volumes. No other way to do this?
785  mesh_.clearOut();
786  }
787 
788  // Reset the instance for if in overwrite mode
789  mesh_.setInstance(timeName());
790  setInstance(mesh_.facesInstance());
791 
792  // Update local mesh data
793  const labelList& oldToNew = map().reverseFaceMap();
794  labelList newSplitFaces(renumber(oldToNew, splitFaces));
795  // Add added faces (every splitFaces becomes two faces)
796  label sz = newSplitFaces.size();
797  newSplitFaces.setSize(2*sz);
798  forAll(map().faceMap(), facei)
799  {
800  label oldFacei = map().faceMap()[facei];
801  if (oldToNew[oldFacei] != facei)
802  {
803  // Added face
804  newSplitFaces[sz++] = facei;
805  }
806  }
807 
808  updateMesh(map, newSplitFaces);
809 
810  return map;
811 }
812 
813 
816 //void Foam::meshRefinement::getCoupledRegionMaster
817 //(
818 // const globalIndex& globalCells,
819 // const boolList& blockedFace,
820 // const regionSplit& globalRegion,
821 // Map<label>& regionToMaster
822 //) const
823 //{
824 // const polyBoundaryMesh& patches = mesh_.boundaryMesh();
825 //
826 // forAll(patches, patchi)
827 // {
828 // const polyPatch& pp = patches[patchi];
829 //
830 // if (isA<processorPolyPatch>(pp))
831 // {
832 // forAll(pp, i)
833 // {
834 // label facei = pp.start()+i;
835 //
836 // if (!blockedFace[facei])
837 // {
838 // // Only if there is a connection to the neighbour
839 // // will there be a multi-domain region; if not through
840 // // this face then through another.
841 //
842 // label celli = mesh_.faceOwner()[facei];
843 // label globalCelli = globalCells.toGlobal(celli);
844 //
845 // Map<label>::iterator iter =
846 // regionToMaster.find(globalRegion[celli]);
847 //
848 // if (iter != regionToMaster.end())
849 // {
850 // label master = iter();
851 // iter() = min(master, globalCelli);
852 // }
853 // else
854 // {
855 // regionToMaster.insert
856 // (
857 // globalRegion[celli],
858 // globalCelli
859 // );
860 // }
861 // }
862 // }
863 // }
864 // }
865 //
866 // // Do reduction
867 // Pstream::mapCombineGather(regionToMaster, minEqOp<label>());
868 // Pstream::mapCombineScatter(regionToMaster);
869 //}
870 //
871 //
872 //void Foam::meshRefinement::calcLocalRegions
873 //(
874 // const globalIndex& globalCells,
875 // const labelList& globalRegion,
876 // const Map<label>& coupledRegionToMaster,
877 // const scalarField& cellWeights,
878 //
879 // Map<label>& globalToLocalRegion,
880 // pointField& localPoints,
881 // scalarField& localWeights
882 //) const
883 //{
884 // globalToLocalRegion.resize(globalRegion.size());
885 // DynamicList<point> localCc(globalRegion.size()/2);
886 // DynamicList<scalar> localWts(globalRegion.size()/2);
887 //
888 // forAll(globalRegion, celli)
889 // {
890 // Map<label>::const_iterator fndMaster =
891 // coupledRegionToMaster.find(globalRegion[celli]);
892 //
893 // if (fndMaster != coupledRegionToMaster.end())
894 // {
895 // // Multi-processor region.
896 // if (globalCells.toGlobal(celli) == fndMaster())
897 // {
898 // // I am master. Allocate region for me.
899 // globalToLocalRegion.insert
900 // (
901 // globalRegion[celli],
902 // localCc.size()
903 // );
904 // localCc.append(mesh_.cellCentres()[celli]);
905 // localWts.append(cellWeights[celli]);
906 // }
907 // }
908 // else
909 // {
910 // // Single processor region.
911 // if
912 // (
913 // globalToLocalRegion.insert
914 // (
915 // globalRegion[celli],
916 // localCc.size()
917 // )
918 // )
919 // {
920 // localCc.append(mesh_.cellCentres()[celli]);
921 // localWts.append(cellWeights[celli]);
922 // }
923 // }
924 // }
925 //
926 // localPoints.transfer(localCc);
927 // localWeights.transfer(localWts);
928 //
929 // if (localPoints.size() != globalToLocalRegion.size())
930 // {
931 // FatalErrorInFunction
932 // << "localPoints:" << localPoints.size()
933 // << " globalToLocalRegion:" << globalToLocalRegion.size()
934 // << exit(FatalError);
935 // }
936 //}
937 //
938 //
939 //Foam::label Foam::meshRefinement::getShiftedRegion
940 //(
941 // const globalIndex& indexer,
942 // const Map<label>& globalToLocalRegion,
943 // const Map<label>& coupledRegionToShifted,
944 // const label globalRegion
945 //)
946 //{
947 // Map<label>::const_iterator iter =
948 // globalToLocalRegion.find(globalRegion);
949 //
950 // if (iter != globalToLocalRegion.end())
951 // {
952 // // Region is 'owned' locally. Convert local region index into global.
953 // return indexer.toGlobal(iter());
954 // }
955 // else
956 // {
957 // return coupledRegionToShifted[globalRegion];
958 // }
959 //}
960 //
961 //
963 //void Foam::meshRefinement::addUnique(const label elem, labelList& lst)
964 //{
965 // if (findIndex(lst, elem) == -1)
966 // {
967 // label sz = lst.size();
968 // lst.setSize(sz+1);
969 // lst[sz] = elem;
970 // }
971 //}
972 //
973 //
974 //void Foam::meshRefinement::calcRegionRegions
975 //(
976 // const labelList& globalRegion,
977 // const Map<label>& globalToLocalRegion,
978 // const Map<label>& coupledRegionToMaster,
979 // labelListList& regionRegions
980 //) const
981 //{
982 // // Global region indexing since we now know the shifted regions.
983 // globalIndex shiftIndexer(globalToLocalRegion.size());
984 //
985 // // Redo the coupledRegionToMaster to be in shifted region indexing.
986 // Map<label> coupledRegionToShifted(coupledRegionToMaster.size());
987 // forAllConstIter(Map<label>, coupledRegionToMaster, iter)
988 // {
989 // label region = iter.key();
990 //
991 // Map<label>::const_iterator fndRegion = globalToLocalRegion.find
992 // (region);
993 //
994 // if (fndRegion != globalToLocalRegion.end())
995 // {
996 // // A local cell is master of this region. Get its shifted region.
997 // coupledRegionToShifted.insert
998 // (
999 // region,
1000 // shiftIndexer.toGlobal(fndRegion())
1001 // );
1002 // }
1003 // Pstream::mapCombineGather(coupledRegionToShifted, minEqOp<label>());
1004 // Pstream::mapCombineScatter(coupledRegionToShifted);
1005 // }
1006 //
1007 //
1008 // // Determine region-region connectivity.
1009 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1010 // // This is for all my regions (so my local ones or the ones I am master
1011 // // of) the neighbouring region indices.
1012 //
1013 //
1014 // // Transfer lists.
1015 // PtrList<HashSet<edge, Hash<edge>>> regionConnectivity
1016 // (Pstream::nProcs());
1017 // forAll(regionConnectivity, proci)
1018 // {
1019 // if (proci != Pstream::myProcNo())
1020 // {
1021 // regionConnectivity.set
1022 // (
1023 // proci,
1024 // new HashSet<edge, Hash<edge>>
1025 // (
1026 // coupledRegionToShifted.size()
1027 // / Pstream::nProcs()
1028 // )
1029 // );
1030 // }
1031 // }
1032 //
1033 //
1034 // // Connectivity. For all my local regions the connected regions.
1035 // regionRegions.setSize(globalToLocalRegion.size());
1036 //
1037 // // Add all local connectivity to regionRegions; add all non-local
1038 // // to the transferlists.
1039 // for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1040 // {
1041 // label ownRegion = globalRegion[mesh_.faceOwner()[facei]];
1042 // label neiRegion = globalRegion[mesh_.faceNeighbour()[facei]];
1043 //
1044 // if (ownRegion != neiRegion)
1045 // {
1046 // label shiftOwnRegion = getShiftedRegion
1047 // (
1048 // shiftIndexer,
1049 // globalToLocalRegion,
1050 // coupledRegionToShifted,
1051 // ownRegion
1052 // );
1053 // label shiftNeiRegion = getShiftedRegion
1054 // (
1055 // shiftIndexer,
1056 // globalToLocalRegion,
1057 // coupledRegionToShifted,
1058 // neiRegion
1059 // );
1060 //
1061 //
1062 // // Connection between two regions. Send to owner of region.
1063 // // - is ownRegion 'owned' by me
1064 // // - is neiRegion 'owned' by me
1065 //
1066 // if (shiftIndexer.isLocal(shiftOwnRegion))
1067 // {
1068 // label localI = shiftIndexer.toLocal(shiftOwnRegion);
1069 // addUnique(shiftNeiRegion, regionRegions[localI]);
1070 // }
1071 // else
1072 // {
1073 // label masterProc = shiftIndexer.whichProcID(shiftOwnRegion);
1074 // regionConnectivity[masterProc].insert
1075 // (
1076 // edge(shiftOwnRegion, shiftNeiRegion)
1077 // );
1078 // }
1079 //
1080 // if (shiftIndexer.isLocal(shiftNeiRegion))
1081 // {
1082 // label localI = shiftIndexer.toLocal(shiftNeiRegion);
1083 // addUnique(shiftOwnRegion, regionRegions[localI]);
1084 // }
1085 // else
1086 // {
1087 // label masterProc = shiftIndexer.whichProcID(shiftNeiRegion);
1088 // regionConnectivity[masterProc].insert
1089 // (
1090 // edge(shiftOwnRegion, shiftNeiRegion)
1091 // );
1092 // }
1093 // }
1094 // }
1095 //
1096 //
1097 // // Send
1098 // forAll(regionConnectivity, proci)
1099 // {
1100 // if (proci != Pstream::myProcNo())
1101 // {
1102 // OPstream str(Pstream::commsTypes::blocking, proci);
1103 // str << regionConnectivity[proci];
1104 // }
1105 // }
1106 // // Receive
1107 // forAll(regionConnectivity, proci)
1108 // {
1109 // if (proci != Pstream::myProcNo())
1110 // {
1111 // IPstream str(Pstream::commsTypes::blocking, proci);
1112 // str >> regionConnectivity[proci];
1113 // }
1114 // }
1115 //
1116 // // Add to addressing.
1117 // forAll(regionConnectivity, proci)
1118 // {
1119 // if (proci != Pstream::myProcNo())
1120 // {
1121 // for
1122 // (
1123 // HashSet<edge, Hash<edge>>::const_iterator iter =
1124 // regionConnectivity[proci].begin();
1125 // iter != regionConnectivity[proci].end();
1126 // ++iter
1127 // )
1128 // {
1129 // const edge& e = iter.key();
1130 //
1131 // bool someLocal = false;
1132 // if (shiftIndexer.isLocal(e[0]))
1133 // {
1134 // label localI = shiftIndexer.toLocal(e[0]);
1135 // addUnique(e[1], regionRegions[localI]);
1136 // someLocal = true;
1137 // }
1138 // if (shiftIndexer.isLocal(e[1]))
1139 // {
1140 // label localI = shiftIndexer.toLocal(e[1]);
1141 // addUnique(e[0], regionRegions[localI]);
1142 // someLocal = true;
1143 // }
1144 //
1145 // if (!someLocal)
1146 // {
1147 // FatalErrorInFunction
1148 // << "Received from processor " << proci
1149 // << " connection " << e
1150 // << " where none of the elements is local to me."
1151 // << abort(FatalError);
1152 // }
1153 // }
1154 // }
1155 // }
1156 //}
1157 
1158 
1159 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
1160 
1161 Foam::meshRefinement::meshRefinement
1163  fvMesh& mesh,
1164  const scalar mergeDistance,
1165  const bool overwrite,
1166  const refinementSurfaces& surfaces,
1167  const refinementFeatures& features,
1168  const shellSurfaces& shells
1169 )
1170 :
1171  mesh_(mesh),
1172  mergeDistance_(mergeDistance),
1173  overwrite_(overwrite),
1174  oldInstance_(mesh.pointsInstance()),
1175  surfaces_(surfaces),
1176  features_(features),
1177  shells_(shells),
1178  meshCutter_
1179  (
1180  mesh,
1181  false // do not try to read history.
1182  ),
1183  surfaceIndex_
1184  (
1185  IOobject
1186  (
1187  "surfaceIndex",
1188  mesh_.facesInstance(),
1190  mesh_,
1193  false
1194  ),
1195  labelList(mesh_.nFaces(), -1)
1196  ),
1197  userFaceData_(0)
1198 {
1199  // recalculate intersections for all faces
1200  updateIntersections(identity(mesh_.nFaces()));
1201 }
1202 
1203 
1204 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1205 
1207 {
1208  // Stats on edges to test. Count proc faces only once.
1209  PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
1210 
1211  label nHits = 0;
1212 
1213  forAll(surfaceIndex_, facei)
1214  {
1215  if (surfaceIndex_[facei] >= 0 && isMasterFace.get(facei) == 1)
1216  {
1217  nHits++;
1218  }
1219  }
1220  return nHits;
1221 }
1222 
1223 
1225 //Foam::labelList Foam::meshRefinement::decomposeCombineRegions
1226 //(
1227 // const scalarField& cellWeights,
1228 // const boolList& blockedFace,
1229 // const List<labelPair>& explicitConnections,
1230 // decompositionMethod& decomposer
1231 //) const
1232 //{
1233 // // Determine global regions, separated by blockedFaces
1234 // regionSplit globalRegion(mesh_, blockedFace, explicitConnections);
1235 //
1236 // // Now globalRegion has global region per cell. Problem is that
1237 // // the region might span multiple domains so we want to get
1238 // // a "region master" per domain. Note that multi-processor
1239 // // regions can only occur on cells on coupled patches.
1240 // // Note: since the number of regions does not change by this the
1241 // // process can be seen as just a shift of a region onto a single
1242 // // processor.
1243 //
1244 //
1245 // // Global cell numbering engine
1246 // globalIndex globalCells(mesh_.nCells());
1247 //
1248 //
1249 // // Determine per coupled region the master cell (lowest numbered cell
1250 // // on lowest numbered processor)
1251 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1252 // // (does not determine master for non-coupled=fully-local regions)
1253 //
1254 // Map<label> coupledRegionToMaster(mesh_.nFaces()-mesh_.nInternalFaces());
1255 // getCoupledRegionMaster
1256 // (
1257 // globalCells,
1258 // blockedFace,
1259 // globalRegion,
1260 // coupledRegionToMaster
1261 // );
1262 //
1263 // // Determine my regions
1264 // // ~~~~~~~~~~~~~~~~~~~~
1265 // // These are the fully local ones or the coupled ones of which I am master
1266 //
1267 // Map<label> globalToLocalRegion;
1268 // pointField localPoints;
1269 // scalarField localWeights;
1270 // calcLocalRegions
1271 // (
1272 // globalCells,
1273 // globalRegion,
1274 // coupledRegionToMaster,
1275 // cellWeights,
1276 //
1277 // globalToLocalRegion,
1278 // localPoints,
1279 // localWeights
1280 // );
1281 //
1282 //
1283 //
1284 // // Find distribution for regions
1285 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1286 //
1287 // labelList regionDistribution;
1288 //
1289 // if (isA<geomDecomp>(decomposer))
1290 // {
1291 // regionDistribution = decomposer.decompose(localPoints, localWeights);
1292 // }
1293 // else
1294 // {
1295 // labelListList regionRegions;
1296 // calcRegionRegions
1297 // (
1298 // globalRegion,
1299 // globalToLocalRegion,
1300 // coupledRegionToMaster,
1301 //
1302 // regionRegions
1303 // );
1304 //
1305 // regionDistribution = decomposer.decompose
1306 // (
1307 // regionRegions,
1308 // localPoints,
1309 // localWeights
1310 // );
1311 // }
1312 //
1313 //
1314 //
1315 // // Convert region-based decomposition back to cell-based one
1316 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1317 //
1318 // // Transfer destination processor back to all. Use global reduce for now.
1319 // Map<label> regionToDist(coupledRegionToMaster.size());
1320 // forAllConstIter(Map<label>, coupledRegionToMaster, iter)
1321 // {
1322 // label region = iter.key();
1323 //
1324 // Map<label>::const_iterator regionFnd = globalToLocalRegion.find
1325 // (region);
1326 //
1327 // if (regionFnd != globalToLocalRegion.end())
1328 // {
1329 // // Master cell is local. Store my distribution.
1330 // regionToDist.insert(iter.key(), regionDistribution[regionFnd()]);
1331 // }
1332 // else
1333 // {
1334 // // Master cell is not on this processor. Make sure it is
1335 // // overridden.
1336 // regionToDist.insert(iter.key(), labelMax);
1337 // }
1338 // }
1339 // Pstream::mapCombineGather(regionToDist, minEqOp<label>());
1340 // Pstream::mapCombineScatter(regionToDist);
1341 //
1342 //
1343 // // Determine destination for all cells
1344 // labelList distribution(mesh_.nCells());
1345 //
1346 // forAll(globalRegion, celli)
1347 // {
1348 // Map<label>::const_iterator fndRegion =
1349 // regionToDist.find(globalRegion[celli]);
1350 //
1351 // if (fndRegion != regionToDist.end())
1352 // {
1353 // distribution[celli] = fndRegion();
1354 // }
1355 // else
1356 // {
1357 // // region is local to the processor.
1358 // label localRegionI = globalToLocalRegion[globalRegion[celli]];
1359 //
1360 // distribution[celli] = regionDistribution[localRegionI];
1361 // }
1362 // }
1363 //
1364 // return distribution;
1365 //}
1366 
1367 
1370  const bool keepZoneFaces,
1371  const bool keepBaffles,
1372  const scalarField& cellWeights,
1373  decompositionMethod& decomposer,
1374  fvMeshDistribute& distributor
1375 )
1376 {
1378 
1379  if (Pstream::parRun())
1380  {
1381  // Wanted distribution
1383 
1384 
1385  // Faces where owner and neighbour are not 'connected' so can
1386  // go to different processors.
1387  boolList blockedFace;
1388  label nUnblocked = 0;
1389 
1390  // Faces that move as block onto single processor
1391  PtrList<labelList> specifiedProcessorFaces;
1392  labelList specifiedProcessor;
1393 
1394  // Pairs of baffles
1395  List<labelPair> couples;
1396 
1397  // Constraints from decomposeParDict
1398  decomposer.setConstraints
1399  (
1400  mesh_,
1401  blockedFace,
1402  specifiedProcessorFaces,
1403  specifiedProcessor,
1404  couples
1405  );
1406 
1407 
1408  if (keepZoneFaces || keepBaffles)
1409  {
1410  if (keepZoneFaces)
1411  {
1412  // Determine decomposition to keep/move surface zones
1413  // on one processor. The reason is that snapping will make these
1414  // into baffles, move and convert them back so if they were
1415  // proc boundaries after baffling&moving the points might be no
1416  // longer synchronised so recoupling will fail. To prevent this
1417  // keep owner&neighbour of such a surface zone on the same
1418  // processor.
1419 
1420  const PtrList<surfaceZonesInfo>& surfZones =
1421  surfaces().surfZones();
1422  const faceZoneMesh& fZones = mesh_.faceZones();
1423  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
1424 
1425  // Get faces whose owner and neighbour should stay together,
1426  // i.e. they are not 'blocked'.
1427 
1428  forAll(surfZones, surfI)
1429  {
1430  const word& fzName = surfZones[surfI].faceZoneName();
1431 
1432  if (fzName.size())
1433  {
1434  // Get zone
1435  const faceZone& fZone = fZones[fzName];
1436 
1437  forAll(fZone, i)
1438  {
1439  label facei = fZone[i];
1440  if (blockedFace[facei])
1441  {
1442  if
1443  (
1444  mesh_.isInternalFace(facei)
1445  || pbm[pbm.whichPatch(facei)].coupled()
1446  )
1447  {
1448  blockedFace[facei] = false;
1449  nUnblocked++;
1450  }
1451  }
1452  }
1453  }
1454  }
1455 
1456 
1457  // If the faceZones are not synchronised the blockedFace
1458  // might not be synchronised. If you are sure the faceZones
1459  // are synchronised remove below check.
1461  (
1462  mesh_,
1463  blockedFace,
1464  andEqOp<bool>() // combine operator
1465  );
1466  }
1467  reduce(nUnblocked, sumOp<label>());
1468 
1469  if (keepZoneFaces)
1470  {
1471  Info<< "Found " << nUnblocked
1472  << " zoned faces to keep together." << endl;
1473  }
1474 
1475 
1476  // Extend un-blockedFaces with any cyclics
1477  {
1478  boolList separatedCoupledFace(mesh_.nFaces(), false);
1479  selectSeparatedCoupledFaces(separatedCoupledFace);
1480 
1481  label nSeparated = 0;
1482  forAll(separatedCoupledFace, facei)
1483  {
1484  if (separatedCoupledFace[facei])
1485  {
1486  if (blockedFace[facei])
1487  {
1488  blockedFace[facei] = false;
1489  nSeparated++;
1490  }
1491  }
1492  }
1493  reduce(nSeparated, sumOp<label>());
1494  Info<< "Found " << nSeparated
1495  << " separated coupled faces to keep together." << endl;
1496 
1497  nUnblocked += nSeparated;
1498  }
1499 
1500 
1501  if (keepBaffles)
1502  {
1503  label nBnd = mesh_.nFaces()-mesh_.nInternalFaces();
1504 
1505  labelList coupledFace(mesh_.nFaces(), -1);
1506  {
1507  // Get boundary baffles that need to stay together
1508  List<labelPair> allCouples
1509  (
1511  );
1512 
1513  // Merge with any couples from
1514  // decompositionMethod::setConstraints
1515  forAll(couples, i)
1516  {
1517  const labelPair& baffle = couples[i];
1518  coupledFace[baffle.first()] = baffle.second();
1519  coupledFace[baffle.second()] = baffle.first();
1520  }
1521  forAll(allCouples, i)
1522  {
1523  const labelPair& baffle = allCouples[i];
1524  coupledFace[baffle.first()] = baffle.second();
1525  coupledFace[baffle.second()] = baffle.first();
1526  }
1527  }
1528 
1529  couples.setSize(nBnd);
1530  label nCpl = 0;
1531  forAll(coupledFace, facei)
1532  {
1533  if (coupledFace[facei] != -1 && facei < coupledFace[facei])
1534  {
1535  couples[nCpl++] = labelPair(facei, coupledFace[facei]);
1536  }
1537  }
1538  couples.setSize(nCpl);
1539  }
1540  label nCouples = returnReduce(couples.size(), sumOp<label>());
1541 
1542  if (keepBaffles)
1543  {
1544  Info<< "Found " << nCouples << " baffles to keep together."
1545  << endl;
1546  }
1547 
1548  // if (nUnblocked > 0 || nCouples > 0)
1549  //{
1550  // Info<< "Applying special decomposition to keep baffles"
1551  // << " and zoned faces together." << endl;
1552  //
1553  // distribution = decomposeCombineRegions
1554  // (
1555  // cellWeights,
1556  // blockedFace,
1557  // couples,
1558  // decomposer
1559  // );
1560  //
1561  // labelList nProcCells(distributor.countCells(distribution));
1562  // Pstream::listCombineGather(nProcCells, plusEqOp<label>());
1563  // Pstream::listCombineScatter(nProcCells);
1564  //
1565  // Info<< "Calculated decomposition:" << endl;
1566  // forAll(nProcCells, proci)
1567  // {
1568  // Info<< " " << proci << '\t' << nProcCells[proci]
1569  // << endl;
1570  // }
1571  // Info<< endl;
1572  //}
1573  // else
1574  //{
1575  // // Normal decomposition
1576  // distribution = decomposer.decompose
1577  // (
1578  // mesh_,
1579  // mesh_.cellCentres(),
1580  // cellWeights
1581  // );
1582  //}
1583  }
1584  // else
1585  //{
1586  // // Normal decomposition
1587  // distribution = decomposer.decompose
1588  // (
1589  // mesh_,
1590  // cellWeights
1591  // );
1592  //}
1593 
1594 
1595  // Make sure blockedFace not set on couples
1596  forAll(couples, i)
1597  {
1598  const labelPair& baffle = couples[i];
1599  blockedFace[baffle.first()] = false;
1600  blockedFace[baffle.second()] = false;
1601  }
1602 
1603  distribution = decomposer.decompose
1604  (
1605  mesh_,
1606  cellWeights,
1607  blockedFace,
1608  specifiedProcessorFaces,
1609  specifiedProcessor,
1610  couples // explicit connections
1611  );
1612 
1613  if (debug)
1614  {
1615  labelList nProcCells(distributor.countCells(distribution));
1616  Pout<< "Wanted distribution:" << nProcCells << endl;
1617 
1619  Pstream::listCombineScatter(nProcCells);
1620 
1621  Pout<< "Wanted resulting decomposition:" << endl;
1622  forAll(nProcCells, proci)
1623  {
1624  Pout<< " " << proci << '\t' << nProcCells[proci] << endl;
1625  }
1626  Pout<< endl;
1627  }
1628  // Do actual sending/receiving of mesh
1629  map = distributor.distribute(distribution);
1630 
1631  // Update numbering of meshRefiner
1632  distribute(map);
1633 
1634  // Set correct instance (for if overwrite)
1635  mesh_.setInstance(timeName());
1636  setInstance(mesh_.facesInstance());
1637  }
1638 
1639  return map;
1640 }
1641 
1642 
1644 {
1645  label nBoundaryFaces = 0;
1646 
1647  forAll(surfaceIndex_, facei)
1648  {
1649  if (surfaceIndex_[facei] != -1)
1650  {
1651  nBoundaryFaces++;
1652  }
1653  }
1654 
1655  labelList surfaceFaces(nBoundaryFaces);
1656  nBoundaryFaces = 0;
1657 
1658  forAll(surfaceIndex_, facei)
1659  {
1660  if (surfaceIndex_[facei] != -1)
1661  {
1662  surfaceFaces[nBoundaryFaces++] = facei;
1663  }
1664  }
1665  return surfaceFaces;
1666 }
1667 
1668 
1670 {
1671  const faceList& faces = mesh_.faces();
1672 
1673  // Mark all points on faces that will become baffles
1674  PackedBoolList isBoundaryPoint(mesh_.nPoints());
1675  label nBoundaryPoints = 0;
1676 
1677  forAll(surfaceIndex_, facei)
1678  {
1679  if (surfaceIndex_[facei] != -1)
1680  {
1681  const face& f = faces[facei];
1682 
1683  forAll(f, fp)
1684  {
1685  if (isBoundaryPoint.set(f[fp], 1u))
1686  {
1687  nBoundaryPoints++;
1688  }
1689  }
1690  }
1691  }
1692 
1694  // labelList adaptPatchIDs(meshedPatches());
1695  // forAll(adaptPatchIDs, i)
1696  //{
1697  // label patchi = adaptPatchIDs[i];
1698  //
1699  // if (patchi != -1)
1700  // {
1701  // const polyPatch& pp = mesh_.boundaryMesh()[patchi];
1702  //
1703  // label facei = pp.start();
1704  //
1705  // forAll(pp, i)
1706  // {
1707  // const face& f = faces[facei];
1708  //
1709  // forAll(f, fp)
1710  // {
1711  // if (isBoundaryPoint.set(f[fp], 1u))
1712  // nBoundaryPoints++;
1713  // }
1714  // }
1715  // facei++;
1716  // }
1717  // }
1718  //}
1719 
1720 
1721  // Pack
1722  labelList boundaryPoints(nBoundaryPoints);
1723  nBoundaryPoints = 0;
1724  forAll(isBoundaryPoint, pointi)
1725  {
1726  if (isBoundaryPoint.get(pointi) == 1u)
1727  {
1728  boundaryPoints[nBoundaryPoints++] = pointi;
1729  }
1730  }
1731 
1732  return boundaryPoints;
1733 }
1734 
1735 
1738  const polyMesh& mesh,
1739  const labelList& patchIDs
1740 )
1741 {
1742  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1743 
1744  // Count faces.
1745  label nFaces = 0;
1746 
1747  forAll(patchIDs, i)
1748  {
1749  const polyPatch& pp = patches[patchIDs[i]];
1750 
1751  nFaces += pp.size();
1752  }
1753 
1754  // Collect faces.
1755  labelList addressing(nFaces);
1756  nFaces = 0;
1757 
1758  forAll(patchIDs, i)
1759  {
1760  const polyPatch& pp = patches[patchIDs[i]];
1761 
1762  label meshFacei = pp.start();
1763 
1764  forAll(pp, i)
1765  {
1766  addressing[nFaces++] = meshFacei++;
1767  }
1768  }
1769 
1771  (
1773  (
1774  IndirectList<face>(mesh.faces(), addressing),
1775  mesh.points()
1776  )
1777  );
1778 }
1779 
1780 
1783  const pointMesh& pMesh,
1784  const labelList& adaptPatchIDs
1785 )
1786 {
1787  const polyMesh& mesh = pMesh();
1788 
1789  // Construct displacement field.
1790  const pointBoundaryMesh& pointPatches = pMesh.boundary();
1791 
1792  wordList patchFieldTypes
1793  (
1794  pointPatches.size(),
1795  slipPointPatchVectorField::typeName
1796  );
1797 
1798  forAll(adaptPatchIDs, i)
1799  {
1800  patchFieldTypes[adaptPatchIDs[i]] =
1801  fixedValuePointPatchVectorField::typeName;
1802  }
1803 
1804  forAll(pointPatches, patchi)
1805  {
1806  if (isA<processorPointPatch>(pointPatches[patchi]))
1807  {
1808  patchFieldTypes[patchi] = calculatedPointPatchVectorField::typeName;
1809  }
1810  else if (isA<cyclicPointPatch>(pointPatches[patchi]))
1811  {
1812  patchFieldTypes[patchi] = cyclicSlipPointPatchVectorField::typeName;
1813  }
1814  }
1815 
1816  // Note: time().timeName() instead of meshRefinement::timeName() since
1817  // postprocessable field.
1819  (
1820  new pointVectorField
1821  (
1822  IOobject
1823  (
1824  "pointDisplacement",
1825  mesh.time().timeName(),
1826  mesh,
1829  ),
1830  pMesh,
1831  dimensionedVector("displacement", dimLength, Zero),
1832  patchFieldTypes
1833  )
1834  );
1835  return tfld;
1836 }
1837 
1838 
1840 {
1841  const faceZoneMesh& fZones = mesh.faceZones();
1842 
1843  // Check any zones are present anywhere and in same order
1844 
1845  {
1846  List<wordList> zoneNames(Pstream::nProcs());
1847  zoneNames[Pstream::myProcNo()] = fZones.names();
1848  Pstream::gatherList(zoneNames);
1849  Pstream::scatterList(zoneNames);
1850  // All have same data now. Check.
1851  forAll(zoneNames, proci)
1852  {
1853  if (proci != Pstream::myProcNo())
1854  {
1855  if (zoneNames[proci] != zoneNames[Pstream::myProcNo()])
1856  {
1858  << "faceZones are not synchronised on processors." << nl
1859  << "Processor " << proci << " has faceZones "
1860  << zoneNames[proci] << nl
1861  << "Processor " << Pstream::myProcNo()
1862  << " has faceZones "
1863  << zoneNames[Pstream::myProcNo()] << nl
1864  << exit(FatalError);
1865  }
1866  }
1867  }
1868  }
1869 
1870  // Check that coupled faces are present on both sides.
1871 
1872  labelList faceToZone(mesh.nFaces()-mesh.nInternalFaces(), -1);
1873 
1874  forAll(fZones, zoneI)
1875  {
1876  const faceZone& fZone = fZones[zoneI];
1877 
1878  forAll(fZone, i)
1879  {
1880  label bFacei = fZone[i]-mesh.nInternalFaces();
1881 
1882  if (bFacei >= 0)
1883  {
1884  if (faceToZone[bFacei] == -1)
1885  {
1886  faceToZone[bFacei] = zoneI;
1887  }
1888  else if (faceToZone[bFacei] == zoneI)
1889  {
1891  << "Face " << fZone[i] << " in zone "
1892  << fZone.name()
1893  << " is twice in zone!"
1894  << abort(FatalError);
1895  }
1896  else
1897  {
1899  << "Face " << fZone[i] << " in zone "
1900  << fZone.name()
1901  << " is also in zone "
1902  << fZones[faceToZone[bFacei]].name()
1903  << abort(FatalError);
1904  }
1905  }
1906  }
1907  }
1908 
1909  labelList neiFaceToZone(faceToZone);
1910  syncTools::swapBoundaryFaceList(mesh, neiFaceToZone);
1911 
1912  forAll(faceToZone, i)
1913  {
1914  if (faceToZone[i] != neiFaceToZone[i])
1915  {
1917  << "Face " << mesh.nInternalFaces()+i
1918  << " is in zone " << faceToZone[i]
1919  << ", its coupled face is in zone " << neiFaceToZone[i]
1920  << abort(FatalError);
1921  }
1922  }
1923 }
1924 
1925 
1928  const polyMesh& mesh,
1929  const PackedBoolList& isMasterEdge,
1930  const labelList& meshPoints,
1931  const edgeList& edges,
1932  scalarField& edgeWeights,
1933  scalarField& invSumWeight
1934 )
1935 {
1936  const pointField& pts = mesh.points();
1937 
1938  // Calculate edgeWeights and inverse sum of edge weights
1939  edgeWeights.setSize(isMasterEdge.size());
1940  invSumWeight.setSize(meshPoints.size());
1941 
1942  forAll(edges, edgeI)
1943  {
1944  const edge& e = edges[edgeI];
1945  scalar eMag = max
1946  (
1947  small,
1948  mag
1949  (
1950  pts[meshPoints[e[1]]]
1951  - pts[meshPoints[e[0]]]
1952  )
1953  );
1954  edgeWeights[edgeI] = 1.0/eMag;
1955  }
1956 
1957  // Sum per point all edge weights
1958  weightedSum
1959  (
1960  mesh,
1961  isMasterEdge,
1962  meshPoints,
1963  edges,
1964  edgeWeights,
1965  scalarField(meshPoints.size(), 1.0), // data
1966  invSumWeight
1967  );
1968 
1969  // Inplace invert
1970  forAll(invSumWeight, pointi)
1971  {
1972  scalar w = invSumWeight[pointi];
1973 
1974  if (w > 0.0)
1975  {
1976  invSumWeight[pointi] = 1.0/w;
1977  }
1978  }
1979 }
1980 
1981 
1984  fvMesh& mesh,
1985  const label insertPatchi,
1986  const word& patchName,
1987  const dictionary& patchDict
1988 )
1989 {
1990  // Clear local fields and e.g. polyMesh parallelInfo.
1991  mesh.clearOut();
1992 
1993  polyBoundaryMesh& polyPatches =
1994  const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
1995  fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
1996 
1997  label patchi = polyPatches.size();
1998 
1999  // Add polyPatch at the end
2000  polyPatches.setSize(patchi+1);
2001  polyPatches.set
2002  (
2003  patchi,
2005  (
2006  patchName,
2007  patchDict,
2008  insertPatchi,
2009  polyPatches
2010  )
2011  );
2012  fvPatches.setSize(patchi+1);
2013  fvPatches.set
2014  (
2015  patchi,
2016  fvPatch::New
2017  (
2018  polyPatches[patchi], // point to newly added polyPatch
2019  mesh.boundary()
2020  )
2021  );
2022 
2023  addPatchFields<volScalarField>
2024  (
2025  mesh,
2027  );
2028  addPatchFields<volVectorField>
2029  (
2030  mesh,
2032  );
2033  addPatchFields<volSphericalTensorField>
2034  (
2035  mesh,
2037  );
2038  addPatchFields<volSymmTensorField>
2039  (
2040  mesh,
2042  );
2043  addPatchFields<volTensorField>
2044  (
2045  mesh,
2047  );
2048 
2049  // Surface fields
2050 
2051  addPatchFields<surfaceScalarField>
2052  (
2053  mesh,
2055  );
2056  addPatchFields<surfaceVectorField>
2057  (
2058  mesh,
2060  );
2061  addPatchFields<surfaceSphericalTensorField>
2062  (
2063  mesh,
2065  );
2066  addPatchFields<surfaceSymmTensorField>
2067  (
2068  mesh,
2070  );
2071  addPatchFields<surfaceTensorField>
2072  (
2073  mesh,
2075  );
2076  return patchi;
2077 }
2078 
2079 
2082  fvMesh& mesh,
2083  const word& patchName,
2084  const dictionary& patchInfo
2085 )
2086 {
2087  polyBoundaryMesh& polyPatches =
2088  const_cast<polyBoundaryMesh&>(mesh.boundaryMesh());
2089  fvBoundaryMesh& fvPatches = const_cast<fvBoundaryMesh&>(mesh.boundary());
2090 
2091  const label patchi = polyPatches.findPatchID(patchName);
2092  if (patchi != -1)
2093  {
2094  // Already there
2095  return patchi;
2096  }
2097 
2098 
2099  label insertPatchi = polyPatches.size();
2100  label startFacei = mesh.nFaces();
2101 
2102  forAll(polyPatches, patchi)
2103  {
2104  const polyPatch& pp = polyPatches[patchi];
2105 
2106  if (isA<processorPolyPatch>(pp))
2107  {
2108  insertPatchi = patchi;
2109  startFacei = pp.start();
2110  break;
2111  }
2112  }
2113 
2114  dictionary patchDict(patchInfo);
2115  patchDict.set("nFaces", 0);
2116  patchDict.set("startFace", startFacei);
2117 
2118  // Below is all quite a hack. Feel free to change once there is a better
2119  // mechanism to insert and reorder patches.
2120 
2121  label addedPatchi = appendPatch(mesh, insertPatchi, patchName, patchDict);
2122 
2123 
2124  // Create reordering list
2125  // patches before insert position stay as is
2126  labelList oldToNew(addedPatchi+1);
2127  for (label i = 0; i < insertPatchi; i++)
2128  {
2129  oldToNew[i] = i;
2130  }
2131  // patches after insert position move one up
2132  for (label i = insertPatchi; i < addedPatchi; i++)
2133  {
2134  oldToNew[i] = i+1;
2135  }
2136  // appended patch gets moved to insert position
2137  oldToNew[addedPatchi] = insertPatchi;
2138 
2139  // Shuffle into place
2140  polyPatches.reorder(oldToNew, true);
2141  fvPatches.reorder(oldToNew);
2142 
2143  reorderPatchFields<volScalarField>(mesh, oldToNew);
2144  reorderPatchFields<volVectorField>(mesh, oldToNew);
2145  reorderPatchFields<volSphericalTensorField>(mesh, oldToNew);
2146  reorderPatchFields<volSymmTensorField>(mesh, oldToNew);
2147  reorderPatchFields<volTensorField>(mesh, oldToNew);
2148  reorderPatchFields<surfaceScalarField>(mesh, oldToNew);
2149  reorderPatchFields<surfaceVectorField>(mesh, oldToNew);
2150  reorderPatchFields<surfaceSphericalTensorField>(mesh, oldToNew);
2151  reorderPatchFields<surfaceSymmTensorField>(mesh, oldToNew);
2152  reorderPatchFields<surfaceTensorField>(mesh, oldToNew);
2153 
2154  return insertPatchi;
2155 }
2156 
2157 
2160  const word& name,
2161  const dictionary& patchInfo
2162 )
2163 {
2164  label meshedI = findIndex(meshedPatches_, name);
2165 
2166  if (meshedI != -1)
2167  {
2168  // Already there. Get corresponding polypatch
2169  return mesh_.boundaryMesh().findPatchID(name);
2170  }
2171  else
2172  {
2173  // Add patch
2174  label patchi = addPatch(mesh_, name, patchInfo);
2175 
2176 // dictionary patchDict(patchInfo);
2177 // patchDict.set("nFaces", 0);
2178 // patchDict.set("startFace", 0);
2179 // autoPtr<polyPatch> ppPtr
2180 // (
2181 // polyPatch::New
2182 // (
2183 // name,
2184 // patchDict,
2185 // 0,
2186 // mesh_.boundaryMesh()
2187 // )
2188 // );
2189 // label patchi = fvMeshTools::addPatch
2190 // (
2191 // mesh_,
2192 // ppPtr(),
2193 // dictionary(), // optional field values
2194 // calculatedFvPatchField<scalar>::typeName,
2195 // true
2196 // );
2197 
2198  // Store
2199  meshedPatches_.append(name);
2200 
2201  return patchi;
2202  }
2203 }
2204 
2205 
2207 {
2208  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2209 
2210  DynamicList<label> patchIDs(meshedPatches_.size());
2211  forAll(meshedPatches_, i)
2212  {
2213  label patchi = patches.findPatchID(meshedPatches_[i]);
2214 
2215  if (patchi == -1)
2216  {
2218  << "Problem : did not find patch " << meshedPatches_[i]
2219  << endl << "Valid patches are " << patches.names()
2220  << abort(FatalError);
2221  }
2222  if (!polyPatch::constraintType(patches[patchi].type()))
2223  {
2224  patchIDs.append(patchi);
2225  }
2226  }
2227 
2228  return patchIDs;
2229 }
2230 
2231 
2233 {
2234  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2235 
2236  forAll(patches, patchi)
2237  {
2238  // Check all coupled. Avoid using .coupled() so we also pick up AMI.
2239  if (isA<coupledPolyPatch>(patches[patchi]))
2240  {
2241  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
2242  (
2243  patches[patchi]
2244  );
2245 
2246  if (cpp.separated() || !cpp.parallel())
2247  {
2248  forAll(cpp, i)
2249  {
2250  selected[cpp.start()+i] = true;
2251  }
2252  }
2253  }
2254  }
2255 }
2256 
2257 
2260  const polyMesh& mesh,
2261  const labelList& cellToRegion,
2262  const vector& perturbVec,
2263  const point& p
2264 )
2265 {
2266  label regionI = -1;
2267  label celli = mesh.findCell(p);
2268  if (celli != -1)
2269  {
2270  regionI = cellToRegion[celli];
2271  }
2272  reduce(regionI, maxOp<label>());
2273 
2274  if (regionI == -1)
2275  {
2276  // See if we can perturb a bit
2277  celli = mesh.findCell(p+perturbVec);
2278  if (celli != -1)
2279  {
2280  regionI = cellToRegion[celli];
2281  }
2282  reduce(regionI, maxOp<label>());
2283  }
2284  return regionI;
2285 }
2286 
2287 
2290  const labelList& globalToMasterPatch,
2291  const labelList& globalToSlavePatch,
2292  const point& keepPoint
2293 )
2294 {
2295  // Force calculation of face decomposition (used in findCell)
2296  (void)mesh_.tetBasePtIs();
2297 
2298  // Determine connected regions. regionSplit is the labelList with the
2299  // region per cell.
2300 
2301  boolList blockedFace(mesh_.nFaces(), false);
2302  selectSeparatedCoupledFaces(blockedFace);
2303 
2304  regionSplit cellRegion(mesh_, blockedFace);
2305 
2306  label regionI = findRegion
2307  (
2308  mesh_,
2309  cellRegion,
2310  mergeDistance_*vector(1,1,1), // note:1,1,1 should really be normalised
2311  keepPoint
2312  );
2313 
2314  if (regionI == -1)
2315  {
2317  << "Point " << keepPoint
2318  << " is not inside the mesh." << nl
2319  << "Bounding box of the mesh:" << mesh_.bounds()
2320  << exit(FatalError);
2321  }
2322 
2323  // Subset
2324  // ~~~~~~
2325 
2326  // Get cells to remove
2327  DynamicList<label> cellsToRemove(mesh_.nCells());
2328  forAll(cellRegion, celli)
2329  {
2330  if (cellRegion[celli] != regionI)
2331  {
2332  cellsToRemove.append(celli);
2333  }
2334  }
2335  cellsToRemove.shrink();
2336 
2337  label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
2338  reduce(nCellsToKeep, sumOp<label>());
2339 
2340  Info<< "Keeping all cells in region " << regionI
2341  << " containing point " << keepPoint << endl
2342  << "Selected for keeping : "
2343  << nCellsToKeep
2344  << " cells." << endl;
2345 
2346 
2347  // Remove cells
2348  removeCells cellRemover(mesh_);
2349 
2350  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
2351  labelList exposedPatch;
2352 
2353  label nExposedFaces = returnReduce(exposedFaces.size(), sumOp<label>());
2354  if (nExposedFaces)
2355  {
2356  // FatalErrorInFunction
2357  // << "Removing non-reachable cells should only expose"
2358  // << " boundary faces" << nl
2359  // << "ExposedFaces:" << exposedFaces << abort(FatalError);
2360 
2361  // Patch for exposed faces for lack of anything sensible.
2362  label defaultPatch = 0;
2363  if (globalToMasterPatch.size())
2364  {
2365  defaultPatch = globalToMasterPatch[0];
2366  }
2367 
2369  << "Removing non-reachable cells exposes "
2370  << nExposedFaces << " internal or coupled faces." << endl
2371  << " These get put into patch " << defaultPatch << endl;
2372  exposedPatch.setSize(exposedFaces.size(), defaultPatch);
2373  }
2374 
2375  return doRemoveCells
2376  (
2377  cellsToRemove,
2378  exposedFaces,
2379  exposedPatch,
2380  cellRemover
2381  );
2382 }
2383 
2384 
2386 {
2387  // mesh_ already distributed; distribute my member data
2388 
2389  // surfaceQueries_ ok.
2390 
2391  // refinement
2392  meshCutter_.distribute(map);
2393 
2394  // surfaceIndex is face data.
2395  map.distributeFaceData(surfaceIndex_);
2396 
2397  // maintainedFaces are indices of faces.
2398  forAll(userFaceData_, i)
2399  {
2400  map.distributeFaceData(userFaceData_[i].second());
2401  }
2402 
2403  // Redistribute surface and any fields on it.
2404  {
2405  // Get local mesh bounding box. Single box for now.
2407  treeBoundBox& bb = meshBb[0];
2408  bb = treeBoundBox(mesh_.points()).extend(1e-4);
2409 
2410  // Distribute all geometry (so refinementSurfaces and shellSurfaces)
2411  searchableSurfaces& geometry =
2412  const_cast<searchableSurfaces&>(surfaces_.geometry());
2413 
2414  forAll(geometry, i)
2415  {
2417  autoPtr<mapDistribute> pointMap;
2418  geometry[i].distribute
2419  (
2420  meshBb,
2421  false, // do not keep outside triangles
2422  faceMap,
2423  pointMap
2424  );
2425 
2426  if (faceMap.valid())
2427  {
2428  // (ab)use the instance() to signal current modification time
2429  geometry[i].instance() = geometry[i].time().timeName();
2430  }
2431 
2432  faceMap.clear();
2433  pointMap.clear();
2434  }
2435  }
2436 }
2437 
2438 
2441  const mapPolyMesh& map,
2442  const labelList& changedFaces
2443 )
2444 {
2445  Map<label> dummyMap(0);
2446 
2447  updateMesh(map, changedFaces, dummyMap, dummyMap, dummyMap);
2448 }
2449 
2450 
2453  const labelList& pointsToStore,
2454  const labelList& facesToStore,
2455  const labelList& cellsToStore
2456 )
2457 {
2458  // For now only meshCutter has storable/retrievable data.
2459  meshCutter_.storeData
2460  (
2461  pointsToStore,
2462  facesToStore,
2463  cellsToStore
2464  );
2465 }
2466 
2467 
2470  const mapPolyMesh& map,
2471  const labelList& changedFaces,
2472  const Map<label>& pointsToRestore,
2473  const Map<label>& facesToRestore,
2474  const Map<label>& cellsToRestore
2475 )
2476 {
2477  // For now only meshCutter has storable/retrievable data.
2478 
2479  // Update numbering of cells/vertices.
2480  meshCutter_.updateMesh
2481  (
2482  map,
2483  pointsToRestore,
2484  facesToRestore,
2485  cellsToRestore
2486  );
2487 
2488  // Update surfaceIndex
2489  updateList(map.faceMap(), label(-1), surfaceIndex_);
2490 
2491  // Update cached intersection information
2492  updateIntersections(changedFaces);
2493 
2494  // Update maintained faces
2495  forAll(userFaceData_, i)
2496  {
2497  labelList& data = userFaceData_[i].second();
2498 
2499  if (userFaceData_[i].first() == KEEPALL)
2500  {
2501  // extend list with face-from-face data
2502  updateList(map.faceMap(), label(-1), data);
2503  }
2504  else if (userFaceData_[i].first() == MASTERONLY)
2505  {
2506  // keep master only
2507  labelList newFaceData(map.faceMap().size(), -1);
2508 
2509  forAll(newFaceData, facei)
2510  {
2511  label oldFacei = map.faceMap()[facei];
2512 
2513  if (oldFacei >= 0 && map.reverseFaceMap()[oldFacei] == facei)
2514  {
2515  newFaceData[facei] = data[oldFacei];
2516  }
2517  }
2518  data.transfer(newFaceData);
2519  }
2520  else
2521  {
2522  // remove any face that has been refined i.e. referenced more than
2523  // once.
2524 
2525  // 1. Determine all old faces that get referenced more than once.
2526  // These get marked with -1 in reverseFaceMap
2527  labelList reverseFaceMap(map.reverseFaceMap());
2528 
2529  forAll(map.faceMap(), facei)
2530  {
2531  label oldFacei = map.faceMap()[facei];
2532 
2533  if (oldFacei >= 0)
2534  {
2535  if (reverseFaceMap[oldFacei] != facei)
2536  {
2537  // facei is slave face. Mark old face.
2538  reverseFaceMap[oldFacei] = -1;
2539  }
2540  }
2541  }
2542 
2543  // 2. Map only faces with intact reverseFaceMap
2544  labelList newFaceData(map.faceMap().size(), -1);
2545  forAll(newFaceData, facei)
2546  {
2547  label oldFacei = map.faceMap()[facei];
2548 
2549  if (oldFacei >= 0)
2550  {
2551  if (reverseFaceMap[oldFacei] == facei)
2552  {
2553  newFaceData[facei] = data[oldFacei];
2554  }
2555  }
2556  }
2557  data.transfer(newFaceData);
2558  }
2559  }
2560 }
2561 
2562 
2564 {
2565  bool writeOk = mesh_.write();
2566 
2567  // Make sure that any distributed surfaces (so ones which probably have
2568  // been changed) get written as well.
2569  // Note: should ideally have some 'modified' flag to say whether it
2570  // has been changed or not.
2571  searchableSurfaces& geometry =
2572  const_cast<searchableSurfaces&>(surfaces_.geometry());
2573 
2574  forAll(geometry, i)
2575  {
2576  searchableSurface& s = geometry[i];
2577 
2578  // Check if instance() of surface is not constant or system.
2579  // Is good hint that surface is distributed.
2580  if
2581  (
2582  s.instance() != s.time().system()
2583  && s.instance() != s.time().caseSystem()
2584  && s.instance() != s.time().constant()
2585  && s.instance() != s.time().caseConstant()
2586  )
2587  {
2588  // Make sure it gets written to current time, not constant.
2589  s.instance() = s.time().timeName();
2590  writeOk = writeOk && s.write();
2591  }
2592  }
2593 
2594  return writeOk;
2595 }
2596 
2597 
2600  const polyMesh& mesh,
2601  const labelList& meshPoints
2602 )
2603 {
2604  const globalIndex globalPoints(meshPoints.size());
2605 
2606  labelList myPoints(meshPoints.size());
2607  forAll(meshPoints, pointi)
2608  {
2609  myPoints[pointi] = globalPoints.toGlobal(pointi);
2610  }
2611 
2613  (
2614  mesh,
2615  meshPoints,
2616  myPoints,
2617  minEqOp<label>(),
2618  labelMax
2619  );
2620 
2621 
2622  PackedBoolList isPatchMasterPoint(meshPoints.size());
2623  forAll(meshPoints, pointi)
2624  {
2625  if (myPoints[pointi] == globalPoints.toGlobal(pointi))
2626  {
2627  isPatchMasterPoint[pointi] = true;
2628  }
2629  }
2630 
2631  return isPatchMasterPoint;
2632 }
2633 
2634 
2637  const polyMesh& mesh,
2638  const labelList& meshEdges
2639 )
2640 {
2641  const globalIndex globalEdges(meshEdges.size());
2642 
2643  labelList myEdges(meshEdges.size());
2644  forAll(meshEdges, edgeI)
2645  {
2646  myEdges[edgeI] = globalEdges.toGlobal(edgeI);
2647  }
2648 
2650  (
2651  mesh,
2652  meshEdges,
2653  myEdges,
2654  minEqOp<label>(),
2655  labelMax
2656  );
2657 
2658 
2659  PackedBoolList isMasterEdge(meshEdges.size());
2660  forAll(meshEdges, edgeI)
2661  {
2662  if (myEdges[edgeI] == globalEdges.toGlobal(edgeI))
2663  {
2664  isMasterEdge[edgeI] = true;
2665  }
2666  }
2667 
2668  return isMasterEdge;
2669 }
2670 
2671 
2672 void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
2673 const
2674 {
2675  const globalMeshData& pData = mesh_.globalData();
2676 
2677  if (debug)
2678  {
2679  Pout<< msg.c_str()
2680  << " : cells(local):" << mesh_.nCells()
2681  << " faces(local):" << mesh_.nFaces()
2682  << " points(local):" << mesh_.nPoints()
2683  << endl;
2684  }
2685 
2686  {
2687  PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
2688  label nMasterFaces = 0;
2689  forAll(isMasterFace, i)
2690  {
2691  if (isMasterFace[i])
2692  {
2693  nMasterFaces++;
2694  }
2695  }
2696 
2697  PackedBoolList isMeshMasterPoint(syncTools::getMasterPoints(mesh_));
2698  label nMasterPoints = 0;
2699  forAll(isMeshMasterPoint, i)
2700  {
2701  if (isMeshMasterPoint[i])
2702  {
2703  nMasterPoints++;
2704  }
2705  }
2706 
2707  Info<< msg.c_str()
2708  << " : cells:" << pData.nTotalCells()
2709  << " faces:" << returnReduce(nMasterFaces, sumOp<label>())
2710  << " points:" << returnReduce(nMasterPoints, sumOp<label>())
2711  << endl;
2712  }
2713 
2714 
2715  // if (debug)
2716  {
2717  const labelList& cellLevel = meshCutter_.cellLevel();
2718 
2719  labelList nCells(gMax(cellLevel)+1, 0);
2720 
2721  forAll(cellLevel, celli)
2722  {
2723  nCells[cellLevel[celli]]++;
2724  }
2725 
2728 
2729  Info<< "Cells per refinement level:" << endl;
2730  forAll(nCells, levelI)
2731  {
2732  Info<< " " << levelI << '\t' << nCells[levelI]
2733  << endl;
2734  }
2735  }
2736 }
2737 
2738 
2740 {
2741  if (overwrite_ && mesh_.time().timeIndex() == 0)
2742  {
2743  return oldInstance_;
2744  }
2745  else
2746  {
2747  return mesh_.time().timeName();
2748  }
2749 }
2750 
2751 
2753 {
2754  // Note: use time().timeName(), not meshRefinement::timeName()
2755  // so as to dump the fields to 0, not to constant.
2756  {
2757  volScalarField volRefLevel
2758  (
2759  IOobject
2760  (
2761  "cellLevel",
2762  mesh_.time().timeName(),
2763  mesh_,
2766  false
2767  ),
2768  mesh_,
2769  dimensionedScalar("zero", dimless, 0)
2770  );
2771 
2772  const labelList& cellLevel = meshCutter_.cellLevel();
2773 
2774  forAll(volRefLevel, celli)
2775  {
2776  volRefLevel[celli] = cellLevel[celli];
2777  }
2778 
2779  volRefLevel.write();
2780  }
2781 
2782  // Dump pointLevel
2783  {
2784  const pointMesh& pMesh = pointMesh::New(mesh_);
2785 
2786  pointScalarField pointRefLevel
2787  (
2788  IOobject
2789  (
2790  "pointLevel",
2791  mesh_.time().timeName(),
2792  mesh_,
2795  false
2796  ),
2797  pMesh,
2798  dimensionedScalar("zero", dimless, 0)
2799  );
2800 
2801  const labelList& pointLevel = meshCutter_.pointLevel();
2802 
2803  forAll(pointRefLevel, pointi)
2804  {
2805  pointRefLevel[pointi] = pointLevel[pointi];
2806  }
2807 
2808  pointRefLevel.write();
2809  }
2810 }
2811 
2812 
2814 {
2815  {
2816  const pointField& cellCentres = mesh_.cellCentres();
2817 
2818  OFstream str(prefix + "_edges.obj");
2819  label vertI = 0;
2820  Pout<< "meshRefinement::dumpIntersections :"
2821  << " Writing cellcentre-cellcentre intersections to file "
2822  << str.name() << endl;
2823 
2824 
2825  // Redo all intersections
2826  // ~~~~~~~~~~~~~~~~~~~~~~
2827 
2828  // Get boundary face centre and level. Coupled aware.
2829  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2830  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2831  calcNeighbourData(neiLevel, neiCc);
2832 
2833  labelList intersectionFaces(intersectedFaces());
2834 
2835  // Collect segments we want to test for
2836  pointField start(intersectionFaces.size());
2837  pointField end(intersectionFaces.size());
2838 
2839  forAll(intersectionFaces, i)
2840  {
2841  label facei = intersectionFaces[i];
2842  start[i] = cellCentres[mesh_.faceOwner()[facei]];
2843 
2844  if (mesh_.isInternalFace(facei))
2845  {
2846  end[i] = cellCentres[mesh_.faceNeighbour()[facei]];
2847  }
2848  else
2849  {
2850  end[i] = neiCc[facei-mesh_.nInternalFaces()];
2851  }
2852  }
2853 
2854  // Extend segments a bit
2855  {
2856  const vectorField smallVec(rootSmall*(end-start));
2857  start -= smallVec;
2858  end += smallVec;
2859  }
2860 
2861 
2862  // Do tests in one go
2863  labelList surfaceHit;
2864  List<pointIndexHit> surfaceHitInfo;
2865  surfaces_.findAnyIntersection
2866  (
2867  start,
2868  end,
2869  surfaceHit,
2870  surfaceHitInfo
2871  );
2872 
2873  forAll(intersectionFaces, i)
2874  {
2875  if (surfaceHit[i] != -1)
2876  {
2877  meshTools::writeOBJ(str, start[i]);
2878  vertI++;
2879  meshTools::writeOBJ(str, surfaceHitInfo[i].hitPoint());
2880  vertI++;
2881  meshTools::writeOBJ(str, end[i]);
2882  vertI++;
2883  str << "l " << vertI-2 << ' ' << vertI-1 << nl
2884  << "l " << vertI-1 << ' ' << vertI << nl;
2885  }
2886  }
2887  }
2888 
2889  Pout<< endl;
2890 }
2891 
2892 
2895  const debugType debugFlags,
2896  const writeType writeFlags,
2897  const fileName& prefix
2898 ) const
2899 {
2900  if (writeFlags & WRITEMESH)
2901  {
2902  write();
2903  }
2904 
2905  if (writeFlags && !(writeFlags & NOWRITEREFINEMENT))
2906  {
2907  meshCutter_.write();
2908  surfaceIndex_.write();
2909  }
2910 
2911  if (writeFlags & WRITELEVELS)
2912  {
2913  dumpRefinementLevel();
2914  }
2915 
2916  if (debugFlags & OBJINTERSECTIONS && prefix.size())
2917  {
2918  dumpIntersections(prefix);
2919  }
2920 }
2921 
2922 
2924 {
2925  return writeLevel_;
2926 }
2927 
2928 
2930 {
2931  writeLevel_ = flags;
2932 }
2933 
2934 
2936 {
2937  return outputLevel_;
2938 }
2939 
2940 
2942 {
2943  outputLevel_ = flags;
2944 }
2945 
2946 
2947 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:424
ListType renumber(const labelUList &oldToNew, const ListType &)
Renumber the values (not the indices) of a list.
virtual labelList decompose(const pointField &points, const scalarField &pointWeights)
Return for every coordinate the wanted processor number.
This class separates the mesh into distinct unconnected regions, each of which is then given a label ...
Definition: regionSplit.H:119
void clearOut()
Clear all geometry and addressing.
Definition: fvMesh.C:230
static void scatterList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Scatter data. Reverse of gatherList.
void printMeshInfo(const bool, const string &) const
Print some mesh stats.
static void listCombineScatter(const List< commsStruct > &comms, List< T > &Value, const int tag, const label comm)
Scatter data. Reverse of combineGather.
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
const word & name() const
Return name.
Definition: IOobject.H:297
IOoutputType
Enumeration for what to output.
A class for handling file names.
Definition: fileName.H:69
bool set(const label) const
Is element set.
Definition: PtrListI.H:65
virtual bool separated() const
Are the planes separated.
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:466
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Foam::pointBoundaryMesh.
unsigned int get(const label) const
Get value at index I.
Definition: PackedListI.H:961
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
void setRefinement(const labelList &cellsToRemove, const labelList &facesToExpose, const labelList &patchIDs, polyTopoChange &) const
Play commands into polyTopoChange to remove cells.
Definition: removeCells.C:179
label countHits() const
Count number of intersections (local)
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:137
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
void reorder(const labelUList &)
Reorders elements. Ordering does not have to be done in.
Definition: PtrList.C:197
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
label nInternalFaces() const
static void swapBoundaryFacePositions(const polyMesh &mesh, UList< point > &l)
Swap coupled positions.
Definition: syncTools.H:446
IOdebugType
Enumeration for what to debug.
label nFaces() const
static label findRegion(const polyMesh &, const labelList &cellRegion, const vector &perturbVec, const point &p)
Find region point is in. Uses optional perturbation to re-test.
Various mesh related information for a parallel run. Upon construction, constructs all info using par...
Output to file stream.
Definition: OFstream.H:82
static word meshSubDir
Return the mesh sub-directory name (usually "polyMesh")
Definition: polyMesh.H:312
static const NamedEnum< IOoutputType, 1 > IOoutputTypeNames
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:163
const boolList & flipMap() const
Return face flip map.
Definition: faceZone.H:253
label fcIndex(const label i) const
Return the forward circular index, i.e. the next index.
Definition: UListI.H:58
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:427
Given list of cells to remove insert all the topology changes.
Definition: removeCells.H:59
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
label nTotalCells() const
Return total number of cells in decomposed mesh.
labelList intersectedFaces() const
Get faces with intersection.
static writeType writeLevel()
Get/set write level.
label findPatchID(const word &patchName) const
Find patch index given a name.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
Class containing mesh-to-mesh mapping information after a mesh distribution where we send parts of me...
The coupledPolyPatch is an abstract base class for patches that couple regions of the computational d...
void writeOBJ(Ostream &os, const point &pt)
Write obj representation of point.
Definition: meshTools.C:203
label whichFace(const label globalCellID) const
Helper function to re-direct to zone::localID(...)
Definition: faceZone.C:306
static word timeName(const scalar, const int precision=precision_)
Return time name of given scalar time.
Definition: Time.C:626
Base class of (analytical or triangulated) surface. Encapsulates all the search routines. WIP.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
patches[0]
Pre-declare related SubField type.
Definition: Field.H:61
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
labelList intersectedPoints() const
Get points on surfaces with intersection and boundary faces.
const fileName & name() const
Return the name of the stream.
Definition: OFstream.H:120
Initialise the NamedEnum HashTable from the static list of names.
Definition: NamedEnum.H:51
static void listCombineGather(const List< commsStruct > &comms, List< T > &Value, const CombineOp &cop, const int tag, const label comm)
static PackedBoolList getMasterFaces(const polyMesh &)
Get per face whether it is uncoupled or a master of a.
Definition: syncTools.C:153
scalar f1
Definition: createFields.H:28
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
Pair< int > faceMap(const label facePi, const face &faceP, const label faceNi, const face &faceN)
autoPtr< mapPolyMesh > splitMeshRegions(const labelList &globalToMasterPatch, const labelList &globalToSlavePatch, const point &keepPoint)
Split mesh. Keep part containing point.
static labelPairList findDuplicateFacePairs(const polyMesh &)
Helper routine to find all baffles (two boundary faces.
static const pointMesh & New(const polyMesh &mesh)
Definition: MeshObject.C:44
IOwriteType
Enumeration for what to write.
Container for data on surfaces used for surface-driven refinement. Contains all the data about the le...
static label addPatch(fvMesh &, const word &name, const dictionary &)
Helper:add patch to mesh. Update all registered fields.
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1003
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:48
fileName caseSystem() const
Return system name for the case.
Definition: TimePaths.C:110
static bool constraintType(const word &pt)
Return true if the given type is a constraint type.
Definition: polyPatch.C:246
Encapsulates queries for features.
static label appendPatch(fvMesh &, const label insertPatchi, const word &, const dictionary &)
Helper:append patch to end of mesh.
void clear()
Delete object (if the pointer is valid) and set pointer to.
Definition: autoPtrI.H:126
void distribute(const mapDistributePolyMesh &)
Update local numbering for mesh redistribution.
virtual bool parallel() const
Are the cyclic planes parallel.
void write(Ostream &, const label, const dictionary &)
Write with dictionary lookup.
Calculates a unique integer (label so might not have enough room - 2G max) for processor + local inde...
Definition: globalIndex.H:63
A List obtained as a section of another List.
Definition: SubList.H:53
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){ gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
const fileName & pointsInstance() const
Return the current instance directory for points.
Definition: polyMesh.C:795
dynamicFvMesh & mesh
An ordered pair of two objects of type <T> with first() and second() elements.
Definition: contiguous.H:49
virtual bool coupled() const
Return true if this patch is geometrically coupled (i.e. faces and.
Definition: polyPatch.H:310
void distributeFaceData(List< T > &lst) const
Distribute list of face data.
Definition: ops.H:70
label addFace(const face &f, const label own, const label nei, const label masterPointID, const label masterEdgeID, const label masterFaceID, const bool flipFaceFlux, const label patchID, const label zoneID, const bool zoneFlip)
Add face to cells. Return new face label.
static void calculateEdgeWeights(const polyMesh &mesh, const PackedBoolList &isMasterEdge, const labelList &meshPoints, const edgeList &edges, scalarField &edgeWeights, scalarField &invSumWeight)
Helper: calculate edge weights (1/length)
const labelUList & faceCells() const
Return face-cell addressing.
Definition: polyPatch.C:315
An edge is a list of two point labels. The functionality it provides supports the discretisation on a...
Definition: edge.H:58
bool found(const Key &) const
Return true if hashedEntry is found in table.
Definition: HashTable.C:113
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
autoPtr< mapPolyMesh > splitFaces(const labelList &splitFaces, const labelPairList &splits)
Split faces into two.
A class for handling words, derived from string.
Definition: word.H:59
static PackedBoolList getMasterEdges(const polyMesh &mesh, const labelList &meshEdges)
Determine master edge for subset of edges. If coupled.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
wordList names() const
Return a list of patch names.
void append(const T &)
Append an element at the end of the list.
Definition: ListI.H:184
void setConstraints(const polyMesh &mesh, boolList &blockedFace, PtrList< labelList > &specifiedProcessorFaces, labelList &specifiedProcessor, List< labelPair > &explicitConnections)
Helper: extract constraints:
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
This boundary condition is not designed to be evaluated; it is assmued that the value is assigned via...
static const NamedEnum< IOdebugType, 5 > IOdebugTypeNames
static void checkCoupledFaceZones(const polyMesh &)
Helper function: check that face zones are synced.
Container for searchableSurfaces.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
static const label labelMax
Definition: label.H:62
word timeName
Definition: getTimeIndex.H:3
List< label > labelList
A List of labels.
Definition: labelList.H:56
static const zero Zero
Definition: zero.H:97
Abstract base class for decomposition.
Encapsulates queries for volume refinement (&#39;refine all cells within shell&#39;).
Definition: shellSurfaces.H:52
const word & name() const
Return name.
Definition: zone.H:147
autoPtr< mapDistributePolyMesh > balance(const bool keepZoneFaces, const bool keepBaffles, const scalarField &cellWeights, decompositionMethod &decomposer, fvMeshDistribute &distributor)
Redecompose according to cell count.
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1028
Accumulating histogram of values. Specified bin resolution automatic generation of bins...
Definition: distribution.H:58
void updateMesh(const mapPolyMesh &)
Force recalculation of locally stored data on topological change.
Definition: removeCells.H:115
const PtrList< surfaceZonesInfo > & surfZones() const
bool valid() const
Return true if the autoPtr valid (ie, the pointer is set)
Definition: autoPtrI.H:83
errorManip< error > abort(error &err)
Definition: errorManip.H:131
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: HashTable.H:61
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
void setSize(const label)
Reset size of PtrList. If extending the PtrList, new entries are.
Definition: PtrList.C:131
const word & system() const
Return system name.
Definition: TimePaths.H:114
Foam::polyBoundaryMesh.
List< treeBoundBox > meshBb(1, treeBoundBox(boundBox(coarseMesh.points(), false)).extend(1e-3))
const vectorField::subField faceAreas() const
Return face normals.
Definition: polyPatch.C:290
fileName caseConstant() const
Return constant name for the case.
Definition: TimePaths.C:123
static const char nl
Definition: Ostream.H:265
void checkData()
Debugging: check that all faces still obey start()>end()
Type gMax(const FieldField< Field, Type > &f)
const Time & time() const
Return time.
defineTypeNameAndDebug(combustionModel, 0)
Database for solution data, solver performance and other reduced data.
Definition: data.H:52
label addMeshedPatch(const word &name, const dictionary &)
Add patch originating from meshing. Update meshedPatches_.
static PackedBoolList getMasterPoints(const polyMesh &)
Get per point whether it is uncoupled or a master of a.
Definition: syncTools.C:65
labelList meshedPatches() const
Get patchIDs for patches added in addMeshedPatch.
autoPtr< mapDistributePolyMesh > distribute(const labelList &dist)
Send cells to neighbours according to distribution.
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
static autoPtr< indirectPrimitivePatch > makePatch(const polyMesh &, const labelList &)
Create patch from set of patches.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
static outputType outputLevel()
Get/set output level.
const labelList & faceMap() const
Old face map.
Definition: mapPolyMesh.H:404
fileName::Type type(const fileName &, const bool followLink=true)
Return the file type: DIRECTORY or FILE.
Definition: POSIX.C:481
void updateMesh(const mapPolyMesh &, const labelList &changedFaces)
Update for external change to mesh. changedFaces are in new mesh.
const Type & second() const
Return second.
Definition: Pair.H:99
label size() const
Return the number of elements in the UPtrList.
Definition: UPtrListI.H:29
void setInstance(const fileName &)
Set instance of all local IOobjects.
static void syncFaceList(const polyMesh &mesh, UList< T > &l, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:381
Calculates points shared by more than two processor patches or cyclic patches.
Definition: globalPoints.H:100
void setSize(const label)
Reset size of List.
Definition: List.C:281
void modifyFace(const face &f, const label facei, const label own, const label nei, const bool flipFaceFlux, const label patchID, const label zoneID, const bool zoneFlip)
Modify vertices or cell of face.
A bit-packed bool list.
static bool & parRun()
Is this a parallel run?
Definition: UPstream.H:397
const fileName & instance() const
Definition: IOobject.H:392
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:409
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
void storeData(const labelList &pointsToStore, const labelList &facesToStore, const labelList &cellsToStore)
Signal points/face/cells for which to store data.
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
static void testSyncPointList(const string &msg, const polyMesh &mesh, const List< scalar > &fld)
const dimensionSet dimLength(0, 1, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:50
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const Time & time() const
Return time.
Definition: IOobject.C:367
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:63
Foam::fvBoundaryMesh.
prefixOSstream Pout(cout, "Pout")
Definition: IOstreams.H:53
static void syncBoundaryFacePositions(const polyMesh &mesh, UList< point > &l, const CombineOp &cop)
Synchronize locations on boundary faces only.
Definition: syncTools.H:363
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
Direct mesh changes based on v1.3 polyTopoChange syntax.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:300
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
label findCell(const point &p, const cellDecomposition=CELL_TETS) const
Find cell enclosing this location and return index.
Definition: polyMesh.C:1386
void set(entry *)
Assign a new entry, overwrite any existing entry.
Definition: dictionary.C:941
static labelList findDuplicateFaces(const primitiveMesh &, const labelList &)
Helper routine to find baffles (two boundary faces using the.
void selectSeparatedCoupledFaces(boolList &) const
Select coupled faces that are not collocated.
Standard boundBox + extra functionality for use in octree.
Definition: treeBoundBox.H:87
wordList names() const
Return a list of zone names.
Definition: ZoneMesh.C:256
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
label nPoints() const
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
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:98
label size() const
Number of entries.
Definition: PackedListI.H:711
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:74
labelList getExposedFaces(const labelList &cellsToRemove) const
Get labels of exposed faces.
Definition: removeCells.C:73
static autoPtr< fvPatch > New(const polyPatch &, const fvBoundaryMesh &)
Return a pointer to a new patch created on freestore from polyPatch.
Definition: fvPatchNew.C:33
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &l)
Swap coupled boundary face values.
Definition: syncTools.H:430
A subset of mesh faces organised as a primitive patch.
Definition: faceZone.H:64
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
bool write() const
Write mesh and all data.
static tmp< pointVectorField > makeDisplacementField(const pointMesh &pMesh, const labelList &adaptPatchIDs)
Helper function to make a pointVectorField with correct.
A class for managing temporary objects.
Definition: PtrList.H:53
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
static PackedBoolList getMasterPoints(const polyMesh &mesh, const labelList &meshPoints)
Determine master point for subset of points. If coupled.
virtual bool write(const bool valid=true) const
Write using setting from DB.
void dumpIntersections(const fileName &prefix) const
Debug: Write intersection information to OBJ format.
const pointBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: pointMesh.H:106
static labelList countCells(const labelList &)
Helper function: count cells per processor in wanted distribution.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:92
static void gatherList(const List< commsStruct > &comms, List< T > &Values, const int tag, const label comm)
Gather data but keep individual values separate.
const Type & first() const
Return first.
Definition: Pair.H:87
A List with indirect addressing.
Definition: IndirectList.H:102
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
void dumpRefinementLevel() const
Write refinement level as volScalarFields for postprocessing.
void transfer(List< T > &)
Transfer the contents of the argument List into this list.
Definition: List.C:342
void reorder(const labelUList &, const bool validBoundary)
Reorders patches. Ordering does not have to be done in.
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:495
Namespace for OpenFOAM.
static const NamedEnum< IOwriteType, 5 > IOwriteTypeNames
const fvBoundaryMesh & boundary() const
Return reference to boundary mesh.
Definition: fvMesh.C:545
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:284