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-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 "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 
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 
1629  mesh_.clearOut();
1630 
1631  // Do actual sending/receiving of mesh
1632  map = distributor.distribute(distribution);
1633 
1634  // Update numbering of meshRefiner
1635  distribute(map);
1636 
1637  // Set correct instance (for if overwrite)
1638  mesh_.setInstance(timeName());
1639  setInstance(mesh_.facesInstance());
1640  }
1641 
1642  return map;
1643 }
1644 
1645 
1647 {
1648  label nBoundaryFaces = 0;
1649 
1650  forAll(surfaceIndex_, facei)
1651  {
1652  if (surfaceIndex_[facei] != -1)
1653  {
1654  nBoundaryFaces++;
1655  }
1656  }
1657 
1658  labelList surfaceFaces(nBoundaryFaces);
1659  nBoundaryFaces = 0;
1660 
1661  forAll(surfaceIndex_, facei)
1662  {
1663  if (surfaceIndex_[facei] != -1)
1664  {
1665  surfaceFaces[nBoundaryFaces++] = facei;
1666  }
1667  }
1668  return surfaceFaces;
1669 }
1670 
1671 
1673 {
1674  const faceList& faces = mesh_.faces();
1675 
1676  // Mark all points on faces that will become baffles
1677  PackedBoolList isBoundaryPoint(mesh_.nPoints());
1678  label nBoundaryPoints = 0;
1679 
1680  forAll(surfaceIndex_, facei)
1681  {
1682  if (surfaceIndex_[facei] != -1)
1683  {
1684  const face& f = faces[facei];
1685 
1686  forAll(f, fp)
1687  {
1688  if (isBoundaryPoint.set(f[fp], 1u))
1689  {
1690  nBoundaryPoints++;
1691  }
1692  }
1693  }
1694  }
1695 
1697  // labelList adaptPatchIDs(meshedPatches());
1698  // forAll(adaptPatchIDs, i)
1699  //{
1700  // label patchi = adaptPatchIDs[i];
1701  //
1702  // if (patchi != -1)
1703  // {
1704  // const polyPatch& pp = mesh_.boundaryMesh()[patchi];
1705  //
1706  // label facei = pp.start();
1707  //
1708  // forAll(pp, i)
1709  // {
1710  // const face& f = faces[facei];
1711  //
1712  // forAll(f, fp)
1713  // {
1714  // if (isBoundaryPoint.set(f[fp], 1u))
1715  // nBoundaryPoints++;
1716  // }
1717  // }
1718  // facei++;
1719  // }
1720  // }
1721  //}
1722 
1723 
1724  // Pack
1725  labelList boundaryPoints(nBoundaryPoints);
1726  nBoundaryPoints = 0;
1727  forAll(isBoundaryPoint, pointi)
1728  {
1729  if (isBoundaryPoint.get(pointi) == 1u)
1730  {
1731  boundaryPoints[nBoundaryPoints++] = pointi;
1732  }
1733  }
1734 
1735  return boundaryPoints;
1736 }
1737 
1738 
1741  const polyMesh& mesh,
1742  const labelList& patchIDs
1743 )
1744 {
1745  const polyBoundaryMesh& patches = mesh.boundaryMesh();
1746 
1747  // Count faces.
1748  label nFaces = 0;
1749 
1750  forAll(patchIDs, i)
1751  {
1752  const polyPatch& pp = patches[patchIDs[i]];
1753 
1754  nFaces += pp.size();
1755  }
1756 
1757  // Collect faces.
1758  labelList addressing(nFaces);
1759  nFaces = 0;
1760 
1761  forAll(patchIDs, i)
1762  {
1763  const polyPatch& pp = patches[patchIDs[i]];
1764 
1765  label meshFacei = pp.start();
1766 
1767  forAll(pp, i)
1768  {
1769  addressing[nFaces++] = meshFacei++;
1770  }
1771  }
1772 
1774  (
1776  (
1777  IndirectList<face>(mesh.faces(), addressing),
1778  mesh.points()
1779  )
1780  );
1781 }
1782 
1783 
1786  const pointMesh& pMesh,
1787  const labelList& adaptPatchIDs
1788 )
1789 {
1790  // Construct displacement field.
1791  const pointBoundaryMesh& pointPatches = pMesh.boundary();
1792 
1793  wordList patchFieldTypes
1794  (
1795  pointPatches.size(),
1796  slipPointPatchVectorField::typeName
1797  );
1798 
1799  forAll(adaptPatchIDs, i)
1800  {
1801  patchFieldTypes[adaptPatchIDs[i]] =
1802  fixedValuePointPatchVectorField::typeName;
1803  }
1804 
1805  forAll(pointPatches, patchi)
1806  {
1807  if (isA<processorPointPatch>(pointPatches[patchi]))
1808  {
1809  patchFieldTypes[patchi] = calculatedPointPatchVectorField::typeName;
1810  }
1811  else if (isA<cyclicPointPatch>(pointPatches[patchi]))
1812  {
1813  patchFieldTypes[patchi] = cyclicSlipPointPatchVectorField::typeName;
1814  }
1815  }
1816 
1817  // Note: time().timeName() instead of meshRefinement::timeName() since
1818  // postprocessable field.
1820  (
1822  (
1823  "pointDisplacement",
1824  pMesh,
1826  patchFieldTypes
1827  )
1828  );
1829  return tfld;
1830 }
1831 
1832 
1834 {
1835  const faceZoneMesh& fZones = mesh.faceZones();
1836 
1837  // Check any zones are present anywhere and in same order
1838 
1839  {
1840  List<wordList> zoneNames(Pstream::nProcs());
1841  zoneNames[Pstream::myProcNo()] = fZones.names();
1842  Pstream::gatherList(zoneNames);
1843  Pstream::scatterList(zoneNames);
1844  // All have same data now. Check.
1845  forAll(zoneNames, proci)
1846  {
1847  if (proci != Pstream::myProcNo())
1848  {
1849  if (zoneNames[proci] != zoneNames[Pstream::myProcNo()])
1850  {
1852  << "faceZones are not synchronised on processors." << nl
1853  << "Processor " << proci << " has faceZones "
1854  << zoneNames[proci] << nl
1855  << "Processor " << Pstream::myProcNo()
1856  << " has faceZones "
1857  << zoneNames[Pstream::myProcNo()] << nl
1858  << exit(FatalError);
1859  }
1860  }
1861  }
1862  }
1863 
1864  // Check that coupled faces are present on both sides.
1865 
1866  labelList faceToZone(mesh.nFaces()-mesh.nInternalFaces(), -1);
1867 
1868  forAll(fZones, zoneI)
1869  {
1870  const faceZone& fZone = fZones[zoneI];
1871 
1872  forAll(fZone, i)
1873  {
1874  label bFacei = fZone[i]-mesh.nInternalFaces();
1875 
1876  if (bFacei >= 0)
1877  {
1878  if (faceToZone[bFacei] == -1)
1879  {
1880  faceToZone[bFacei] = zoneI;
1881  }
1882  else if (faceToZone[bFacei] == zoneI)
1883  {
1885  << "Face " << fZone[i] << " in zone "
1886  << fZone.name()
1887  << " is twice in zone!"
1888  << abort(FatalError);
1889  }
1890  else
1891  {
1893  << "Face " << fZone[i] << " in zone "
1894  << fZone.name()
1895  << " is also in zone "
1896  << fZones[faceToZone[bFacei]].name()
1897  << abort(FatalError);
1898  }
1899  }
1900  }
1901  }
1902 
1903  labelList neiFaceToZone(faceToZone);
1904  syncTools::swapBoundaryFaceList(mesh, neiFaceToZone);
1905 
1906  forAll(faceToZone, i)
1907  {
1908  if (faceToZone[i] != neiFaceToZone[i])
1909  {
1911  << "Face " << mesh.nInternalFaces()+i
1912  << " is in zone " << faceToZone[i]
1913  << ", its coupled face is in zone " << neiFaceToZone[i]
1914  << abort(FatalError);
1915  }
1916  }
1917 }
1918 
1919 
1922  const polyMesh& mesh,
1923  const PackedBoolList& isMasterEdge,
1924  const labelList& meshPoints,
1925  const edgeList& edges,
1926  scalarField& edgeWeights,
1927  scalarField& invSumWeight
1928 )
1929 {
1930  const pointField& pts = mesh.points();
1931 
1932  // Calculate edgeWeights and inverse sum of edge weights
1933  edgeWeights.setSize(isMasterEdge.size());
1934  invSumWeight.setSize(meshPoints.size());
1935 
1936  forAll(edges, edgeI)
1937  {
1938  const edge& e = edges[edgeI];
1939  scalar eMag = max
1940  (
1941  small,
1942  mag
1943  (
1944  pts[meshPoints[e[1]]]
1945  - pts[meshPoints[e[0]]]
1946  )
1947  );
1948  edgeWeights[edgeI] = 1.0/eMag;
1949  }
1950 
1951  // Sum per point all edge weights
1952  weightedSum
1953  (
1954  mesh,
1955  isMasterEdge,
1956  meshPoints,
1957  edges,
1958  edgeWeights,
1959  scalarField(meshPoints.size(), 1.0), // data
1960  invSumWeight
1961  );
1962 
1963  // Inplace invert
1964  forAll(invSumWeight, pointi)
1965  {
1966  scalar w = invSumWeight[pointi];
1967 
1968  if (w > 0.0)
1969  {
1970  invSumWeight[pointi] = 1.0/w;
1971  }
1972  }
1973 }
1974 
1975 
1978  const word& name,
1979  const dictionary& patchInfo
1980 )
1981 {
1982  const polyBoundaryMesh& pbm = mesh_.boundaryMesh();
1983 
1984  label meshedI = findIndex(meshedPatches_, name);
1985 
1986  if (meshedI != -1)
1987  {
1988  // Already there. Get corresponding polypatch
1989  return pbm.findPatchID(name);
1990  }
1991  else
1992  {
1993  // Add patch
1994 
1995  label patchi = pbm.findPatchID(name);
1996  if (patchi == -1)
1997  {
1998  patchi = pbm.size();
1999  forAll(pbm, i)
2000  {
2001  const polyPatch& pp = pbm[i];
2002 
2003  if (isA<processorPolyPatch>(pp))
2004  {
2005  patchi = i;
2006  break;
2007  }
2008  }
2009 
2010  dictionary patchDict(patchInfo);
2011  patchDict.set("nFaces", 0);
2012  patchDict.set("startFace", 0);
2013  autoPtr<polyPatch> ppPtr
2014  (
2016  (
2017  name,
2018  patchDict,
2019  0,
2020  pbm
2021  )
2022  );
2023  mesh_.addPatch
2024  (
2025  patchi,
2026  ppPtr(),
2027  dictionary(), // optional field values
2029  true // validBoundary
2030  );
2031  }
2032 
2033  // Store
2034  meshedPatches_.append(name);
2035 
2036  return patchi;
2037  }
2038 }
2039 
2040 
2042 {
2043  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2044 
2045  DynamicList<label> patchIDs(meshedPatches_.size());
2046  forAll(meshedPatches_, i)
2047  {
2048  label patchi = patches.findPatchID(meshedPatches_[i]);
2049 
2050  if (patchi == -1)
2051  {
2053  << "Problem : did not find patch " << meshedPatches_[i]
2054  << endl << "Valid patches are " << patches.names()
2055  << abort(FatalError);
2056  }
2057  if (!polyPatch::constraintType(patches[patchi].type()))
2058  {
2059  patchIDs.append(patchi);
2060  }
2061  }
2062 
2063  return move(patchIDs);
2064 }
2065 
2066 
2068 {
2069  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
2070 
2071  forAll(patches, patchi)
2072  {
2073  // Check all coupled. Avoid using .coupled() so we also pick up AMI.
2074  if (isA<coupledPolyPatch>(patches[patchi]))
2075  {
2076  const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
2077  (
2078  patches[patchi]
2079  );
2080 
2081  if (cpp.separated() || !cpp.parallel())
2082  {
2083  forAll(cpp, i)
2084  {
2085  selected[cpp.start()+i] = true;
2086  }
2087  }
2088  }
2089  }
2090 }
2091 
2092 
2095  const polyMesh& mesh,
2096  const labelList& cellToRegion,
2097  const vector& perturbVec,
2098  const point& p
2099 )
2100 {
2101  label regionI = -1;
2102  label celli = mesh.findCell(p);
2103  if (celli != -1)
2104  {
2105  regionI = cellToRegion[celli];
2106  }
2107  reduce(regionI, maxOp<label>());
2108 
2109  if (regionI == -1)
2110  {
2111  // See if we can perturb a bit
2112  celli = mesh.findCell(p+perturbVec);
2113  if (celli != -1)
2114  {
2115  regionI = cellToRegion[celli];
2116  }
2117  reduce(regionI, maxOp<label>());
2118  }
2119  return regionI;
2120 }
2121 
2122 
2125  const labelList& globalToMasterPatch,
2126  const labelList& globalToSlavePatch,
2127  const point& keepPoint
2128 )
2129 {
2130  // Force calculation of face decomposition (used in findCell)
2131  (void)mesh_.tetBasePtIs();
2132 
2133  // Determine connected regions. regionSplit is the labelList with the
2134  // region per cell.
2135 
2136  boolList blockedFace(mesh_.nFaces(), false);
2137  selectSeparatedCoupledFaces(blockedFace);
2138 
2139  regionSplit cellRegion(mesh_, blockedFace);
2140 
2141  label regionI = findRegion
2142  (
2143  mesh_,
2144  cellRegion,
2145  mergeDistance_*vector(1,1,1), // note:1,1,1 should really be normalised
2146  keepPoint
2147  );
2148 
2149  if (regionI == -1)
2150  {
2152  << "Point " << keepPoint
2153  << " is not inside the mesh." << nl
2154  << "Bounding box of the mesh:" << mesh_.bounds()
2155  << exit(FatalError);
2156  }
2157 
2158  // Subset
2159  // ~~~~~~
2160 
2161  // Get cells to remove
2162  DynamicList<label> cellsToRemove(mesh_.nCells());
2163  forAll(cellRegion, celli)
2164  {
2165  if (cellRegion[celli] != regionI)
2166  {
2167  cellsToRemove.append(celli);
2168  }
2169  }
2170  cellsToRemove.shrink();
2171 
2172  label nCellsToKeep = mesh_.nCells() - cellsToRemove.size();
2173  reduce(nCellsToKeep, sumOp<label>());
2174 
2175  Info<< "Keeping all cells in region " << regionI
2176  << " containing point " << keepPoint << endl
2177  << "Selected for keeping : "
2178  << nCellsToKeep
2179  << " cells." << endl;
2180 
2181 
2182  // Remove cells
2183  removeCells cellRemover(mesh_);
2184 
2185  labelList exposedFaces(cellRemover.getExposedFaces(cellsToRemove));
2186  labelList exposedPatch;
2187 
2188  label nExposedFaces = returnReduce(exposedFaces.size(), sumOp<label>());
2189  if (nExposedFaces)
2190  {
2191  // FatalErrorInFunction
2192  // << "Removing non-reachable cells should only expose"
2193  // << " boundary faces" << nl
2194  // << "ExposedFaces:" << exposedFaces << abort(FatalError);
2195 
2196  // Patch for exposed faces for lack of anything sensible.
2197  label defaultPatch = 0;
2198  if (globalToMasterPatch.size())
2199  {
2200  defaultPatch = globalToMasterPatch[0];
2201  }
2202 
2204  << "Removing non-reachable cells exposes "
2205  << nExposedFaces << " internal or coupled faces." << endl
2206  << " These get put into patch " << defaultPatch << endl;
2207  exposedPatch.setSize(exposedFaces.size(), defaultPatch);
2208  }
2209 
2210  return doRemoveCells
2211  (
2212  cellsToRemove,
2213  exposedFaces,
2214  exposedPatch,
2215  cellRemover
2216  );
2217 }
2218 
2219 
2221 {
2222  // mesh_ already distributed; distribute my member data
2223 
2224  // surfaceQueries_ ok.
2225 
2226  // refinement
2227  meshCutter_.distribute(map);
2228 
2229  // surfaceIndex is face data.
2230  map.distributeFaceData(surfaceIndex_);
2231 
2232  // maintainedFaces are indices of faces.
2233  forAll(userFaceData_, i)
2234  {
2235  map.distributeFaceData(userFaceData_[i].second());
2236  }
2237 
2238  // Redistribute surface and any fields on it.
2239  {
2240  // Get local mesh bounding box. Single box for now.
2242  treeBoundBox& bb = meshBb[0];
2243  bb = treeBoundBox(mesh_.points()).extend(1e-4);
2244 
2245  // Distribute all geometry (so refinementSurfaces and shellSurfaces)
2246  searchableSurfaces& geometry =
2247  const_cast<searchableSurfaces&>(surfaces_.geometry());
2248 
2249  forAll(geometry, i)
2250  {
2252  autoPtr<mapDistribute> pointMap;
2253  geometry[i].distribute
2254  (
2255  meshBb,
2256  false, // do not keep outside triangles
2257  faceMap,
2258  pointMap
2259  );
2260 
2261  if (faceMap.valid())
2262  {
2263  // (ab)use the instance() to signal current modification time
2264  geometry[i].instance() = geometry[i].time().timeName();
2265  }
2266 
2267  faceMap.clear();
2268  pointMap.clear();
2269  }
2270  }
2271 }
2272 
2273 
2276  const mapPolyMesh& map,
2277  const labelList& changedFaces
2278 )
2279 {
2280  Map<label> dummyMap(0);
2281 
2282  updateMesh(map, changedFaces, dummyMap, dummyMap, dummyMap);
2283 }
2284 
2285 
2288  const labelList& pointsToStore,
2289  const labelList& facesToStore,
2290  const labelList& cellsToStore
2291 )
2292 {
2293  // For now only meshCutter has storable/retrievable data.
2294  meshCutter_.storeData
2295  (
2296  pointsToStore,
2297  facesToStore,
2298  cellsToStore
2299  );
2300 }
2301 
2302 
2305  const mapPolyMesh& map,
2306  const labelList& changedFaces,
2307  const Map<label>& pointsToRestore,
2308  const Map<label>& facesToRestore,
2309  const Map<label>& cellsToRestore
2310 )
2311 {
2312  // For now only meshCutter has storable/retrievable data.
2313 
2314  // Update numbering of cells/vertices.
2315  meshCutter_.updateMesh
2316  (
2317  map,
2318  pointsToRestore,
2319  facesToRestore,
2320  cellsToRestore
2321  );
2322 
2323  // Update surfaceIndex
2324  updateList(map.faceMap(), label(-1), surfaceIndex_);
2325 
2326  // Update cached intersection information
2327  updateIntersections(changedFaces);
2328 
2329  // Update maintained faces
2330  forAll(userFaceData_, i)
2331  {
2332  labelList& data = userFaceData_[i].second();
2333 
2334  if (userFaceData_[i].first() == KEEPALL)
2335  {
2336  // extend list with face-from-face data
2337  updateList(map.faceMap(), label(-1), data);
2338  }
2339  else if (userFaceData_[i].first() == MASTERONLY)
2340  {
2341  // keep master only
2342  labelList newFaceData(map.faceMap().size(), -1);
2343 
2344  forAll(newFaceData, facei)
2345  {
2346  label oldFacei = map.faceMap()[facei];
2347 
2348  if (oldFacei >= 0 && map.reverseFaceMap()[oldFacei] == facei)
2349  {
2350  newFaceData[facei] = data[oldFacei];
2351  }
2352  }
2353  data.transfer(newFaceData);
2354  }
2355  else
2356  {
2357  // remove any face that has been refined i.e. referenced more than
2358  // once.
2359 
2360  // 1. Determine all old faces that get referenced more than once.
2361  // These get marked with -1 in reverseFaceMap
2362  labelList reverseFaceMap(map.reverseFaceMap());
2363 
2364  forAll(map.faceMap(), facei)
2365  {
2366  label oldFacei = map.faceMap()[facei];
2367 
2368  if (oldFacei >= 0)
2369  {
2370  if (reverseFaceMap[oldFacei] != facei)
2371  {
2372  // facei is slave face. Mark old face.
2373  reverseFaceMap[oldFacei] = -1;
2374  }
2375  }
2376  }
2377 
2378  // 2. Map only faces with intact reverseFaceMap
2379  labelList newFaceData(map.faceMap().size(), -1);
2380  forAll(newFaceData, facei)
2381  {
2382  label oldFacei = map.faceMap()[facei];
2383 
2384  if (oldFacei >= 0)
2385  {
2386  if (reverseFaceMap[oldFacei] == facei)
2387  {
2388  newFaceData[facei] = data[oldFacei];
2389  }
2390  }
2391  }
2392  data.transfer(newFaceData);
2393  }
2394  }
2395 }
2396 
2397 
2399 {
2400  bool writeOk = mesh_.write();
2401 
2402  // Make sure that any distributed surfaces (so ones which probably have
2403  // been changed) get written as well.
2404  // Note: should ideally have some 'modified' flag to say whether it
2405  // has been changed or not.
2406  searchableSurfaces& geometry =
2407  const_cast<searchableSurfaces&>(surfaces_.geometry());
2408 
2409  forAll(geometry, i)
2410  {
2411  searchableSurface& s = geometry[i];
2412 
2413  // Check if instance() of surface is not constant or system.
2414  // Is good hint that surface is distributed.
2415  if
2416  (
2417  s.instance() != s.time().system()
2418  && s.instance() != s.time().caseSystem()
2419  && s.instance() != s.time().constant()
2420  && s.instance() != s.time().caseConstant()
2421  )
2422  {
2423  // Make sure it gets written to current time, not constant.
2424  s.instance() = s.time().timeName();
2425  writeOk = writeOk && s.write();
2426  }
2427  }
2428 
2429  return writeOk;
2430 }
2431 
2432 
2435  const polyMesh& mesh,
2436  const labelList& meshPoints
2437 )
2438 {
2439  const globalIndex globalPoints(meshPoints.size());
2440 
2441  labelList myPoints(meshPoints.size());
2442  forAll(meshPoints, pointi)
2443  {
2444  myPoints[pointi] = globalPoints.toGlobal(pointi);
2445  }
2446 
2448  (
2449  mesh,
2450  meshPoints,
2451  myPoints,
2452  minEqOp<label>(),
2453  labelMax
2454  );
2455 
2456 
2457  PackedBoolList isPatchMasterPoint(meshPoints.size());
2458  forAll(meshPoints, pointi)
2459  {
2460  if (myPoints[pointi] == globalPoints.toGlobal(pointi))
2461  {
2462  isPatchMasterPoint[pointi] = true;
2463  }
2464  }
2465 
2466  return isPatchMasterPoint;
2467 }
2468 
2469 
2472  const polyMesh& mesh,
2473  const labelList& meshEdges
2474 )
2475 {
2476  const globalIndex globalEdges(meshEdges.size());
2477 
2478  labelList myEdges(meshEdges.size());
2479  forAll(meshEdges, edgeI)
2480  {
2481  myEdges[edgeI] = globalEdges.toGlobal(edgeI);
2482  }
2483 
2485  (
2486  mesh,
2487  meshEdges,
2488  myEdges,
2489  minEqOp<label>(),
2490  labelMax
2491  );
2492 
2493 
2494  PackedBoolList isMasterEdge(meshEdges.size());
2495  forAll(meshEdges, edgeI)
2496  {
2497  if (myEdges[edgeI] == globalEdges.toGlobal(edgeI))
2498  {
2499  isMasterEdge[edgeI] = true;
2500  }
2501  }
2502 
2503  return isMasterEdge;
2504 }
2505 
2506 
2507 void Foam::meshRefinement::printMeshInfo(const bool debug, const string& msg)
2508 const
2509 {
2510  const globalMeshData& pData = mesh_.globalData();
2511 
2512  if (debug)
2513  {
2514  Pout<< msg.c_str()
2515  << " : cells(local):" << mesh_.nCells()
2516  << " faces(local):" << mesh_.nFaces()
2517  << " points(local):" << mesh_.nPoints()
2518  << endl;
2519  }
2520 
2521  {
2522  PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
2523  label nMasterFaces = 0;
2524  forAll(isMasterFace, i)
2525  {
2526  if (isMasterFace[i])
2527  {
2528  nMasterFaces++;
2529  }
2530  }
2531 
2532  PackedBoolList isMeshMasterPoint(syncTools::getMasterPoints(mesh_));
2533  label nMasterPoints = 0;
2534  forAll(isMeshMasterPoint, i)
2535  {
2536  if (isMeshMasterPoint[i])
2537  {
2538  nMasterPoints++;
2539  }
2540  }
2541 
2542  Info<< msg.c_str()
2543  << " : cells:" << pData.nTotalCells()
2544  << " faces:" << returnReduce(nMasterFaces, sumOp<label>())
2545  << " points:" << returnReduce(nMasterPoints, sumOp<label>())
2546  << endl;
2547  }
2548 
2549 
2550  // if (debug)
2551  {
2552  const labelList& cellLevel = meshCutter_.cellLevel();
2553 
2554  labelList nCells(gMax(cellLevel)+1, 0);
2555 
2556  forAll(cellLevel, celli)
2557  {
2558  nCells[cellLevel[celli]]++;
2559  }
2560 
2563 
2564  Info<< "Cells per refinement level:" << endl;
2565  forAll(nCells, levelI)
2566  {
2567  Info<< " " << levelI << '\t' << nCells[levelI]
2568  << endl;
2569  }
2570  }
2571 }
2572 
2573 
2575 {
2576  if (overwrite_ && mesh_.time().timeIndex() == 0)
2577  {
2578  return oldInstance_;
2579  }
2580  else
2581  {
2582  return mesh_.time().timeName();
2583  }
2584 }
2585 
2586 
2588 {
2589  // Note: use time().timeName(), not meshRefinement::timeName()
2590  // so as to dump the fields to 0, not to constant.
2591  {
2592  volScalarField volRefLevel
2593  (
2594  IOobject
2595  (
2596  "cellLevel",
2597  mesh_.time().timeName(),
2598  mesh_,
2601  false
2602  ),
2603  mesh_,
2605  );
2606 
2607  const labelList& cellLevel = meshCutter_.cellLevel();
2608 
2609  forAll(volRefLevel, celli)
2610  {
2611  volRefLevel[celli] = cellLevel[celli];
2612  }
2613 
2614  volRefLevel.write();
2615  }
2616 
2617  // Dump pointLevel
2618  {
2619  const pointMesh& pMesh = pointMesh::New(mesh_);
2620 
2621  pointScalarField pointRefLevel
2622  (
2623  IOobject
2624  (
2625  "pointLevel",
2626  mesh_.time().timeName(),
2627  mesh_,
2630  false
2631  ),
2632  pMesh,
2634  );
2635 
2636  const labelList& pointLevel = meshCutter_.pointLevel();
2637 
2638  forAll(pointRefLevel, pointi)
2639  {
2640  pointRefLevel[pointi] = pointLevel[pointi];
2641  }
2642 
2643  pointRefLevel.write();
2644  }
2645 }
2646 
2647 
2649 {
2650  {
2651  const pointField& cellCentres = mesh_.cellCentres();
2652 
2653  OFstream str(prefix + "_edges.obj");
2654  label vertI = 0;
2655  Pout<< "meshRefinement::dumpIntersections :"
2656  << " Writing cellcentre-cellcentre intersections to file "
2657  << str.name() << endl;
2658 
2659 
2660  // Redo all intersections
2661  // ~~~~~~~~~~~~~~~~~~~~~~
2662 
2663  // Get boundary face centre and level. Coupled aware.
2664  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2665  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2666  calcNeighbourData(neiLevel, neiCc);
2667 
2668  labelList intersectionFaces(intersectedFaces());
2669 
2670  // Collect segments we want to test for
2671  pointField start(intersectionFaces.size());
2672  pointField end(intersectionFaces.size());
2673 
2674  forAll(intersectionFaces, i)
2675  {
2676  label facei = intersectionFaces[i];
2677  start[i] = cellCentres[mesh_.faceOwner()[facei]];
2678 
2679  if (mesh_.isInternalFace(facei))
2680  {
2681  end[i] = cellCentres[mesh_.faceNeighbour()[facei]];
2682  }
2683  else
2684  {
2685  end[i] = neiCc[facei-mesh_.nInternalFaces()];
2686  }
2687  }
2688 
2689  // Extend segments a bit
2690  {
2691  const vectorField smallVec(rootSmall*(end-start));
2692  start -= smallVec;
2693  end += smallVec;
2694  }
2695 
2696 
2697  // Do tests in one go
2698  labelList surfaceHit;
2699  List<pointIndexHit> surfaceHitInfo;
2700  surfaces_.findAnyIntersection
2701  (
2702  start,
2703  end,
2704  surfaceHit,
2705  surfaceHitInfo
2706  );
2707 
2708  forAll(intersectionFaces, i)
2709  {
2710  if (surfaceHit[i] != -1)
2711  {
2712  meshTools::writeOBJ(str, start[i]);
2713  vertI++;
2714  meshTools::writeOBJ(str, surfaceHitInfo[i].hitPoint());
2715  vertI++;
2716  meshTools::writeOBJ(str, end[i]);
2717  vertI++;
2718  str << "l " << vertI-2 << ' ' << vertI-1 << nl
2719  << "l " << vertI-1 << ' ' << vertI << nl;
2720  }
2721  }
2722  }
2723 
2724  Pout<< endl;
2725 }
2726 
2727 
2730  const debugType debugFlags,
2731  const writeType writeFlags,
2732  const fileName& prefix
2733 ) const
2734 {
2735  if (writeFlags & WRITEMESH)
2736  {
2737  write();
2738  }
2739 
2740  if (writeFlags && !(writeFlags & NOWRITEREFINEMENT))
2741  {
2742  meshCutter_.write();
2743  surfaceIndex_.write();
2744  }
2745 
2746  if (writeFlags & WRITELEVELS)
2747  {
2748  dumpRefinementLevel();
2749  }
2750 
2751  if (debugFlags & OBJINTERSECTIONS && prefix.size())
2752  {
2753  dumpIntersections(prefix);
2754  }
2755 }
2756 
2757 
2759 {
2760  return writeLevel_;
2761 }
2762 
2763 
2765 {
2766  writeLevel_ = flags;
2767 }
2768 
2769 
2771 {
2772  return outputLevel_;
2773 }
2774 
2775 
2777 {
2778  outputLevel_ = flags;
2779 }
2780 
2781 
2782 // ************************************************************************* //
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:434
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
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: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
const word & name() const
Return name.
Definition: IOobject.H:295
IOoutputType
Enumeration for what to output.
A class for handling file names.
Definition: fileName.H:79
virtual bool separated() const
Are the planes separated.
const faceZoneMesh & faceZones() const
Return face zone mesh.
Definition: polyMesh.H:476
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:954
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:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#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:248
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:429
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.
Abstract base class with a fat-interface to all derived classes covering all possible ways in which t...
Definition: fvPatchField.H:66
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:60
PrimitivePatch< IndirectList< face >, 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...
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1131
Mesh representing a set of points created from polyMesh.
Definition: pointMesh.H:48
static tmp< GeometricField< vector, pointPatchField, pointMesh > > New(const word &name, const Mesh &, const dimensionSet &, const word &patchFieldType=pointPatchField< vector >::calculatedType())
Return a temporary field constructed from name, mesh, dimensionSet.
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.
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:802
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:313
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:177
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
static const NamedEnum< IOdebugType, 5 > IOdebugTypeNames
meshRefinement(fvMesh &mesh, const scalar mergeDistance, const bool overwrite, const refinementSurfaces &, const refinementFeatures &, const shellSurfaces &)
Construct from components.
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:1156
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:60
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
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)
defineTypeNameAndDebug(combustionModel, 0)
Database for solution and other reduced data.
Definition: data.H:51
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:399
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:399
const fileName & instance() const
Definition: IOobject.H:378
static label nProcs(const label communicator=0)
Number of processes in parallel run.
Definition: UPstream.H:411
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:360
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
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:303
fileType type(const fileName &, const bool checkVariants=true, const bool followLink=true)
Return the file type: directory or file.
Definition: POSIX.C:488
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:1606
void set(entry *)
Assign a new entry, overwrite any existing entry.
Definition: dictionary.C:948
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 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.
virtual bool write(const bool write=true) const
Write using setting from DB.
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.
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:103
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:101
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
const labelList & reverseFaceMap() const
Reverse face map.
Definition: mapPolyMesh.H:490
Namespace for OpenFOAM.
static const NamedEnum< IOwriteType, 5 > IOwriteTypeNames
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:284