meshRefinementRefine.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | www.openfoam.com
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8  Copyright (C) 2011-2017 OpenFOAM Foundation
9  Copyright (C) 2015-2021 OpenCFD Ltd.
10 -------------------------------------------------------------------------------
11 License
12  This file is part of OpenFOAM.
13 
14  OpenFOAM is free software: you can redistribute it and/or modify it
15  under the terms of the GNU General Public License as published by
16  the Free Software Foundation, either version 3 of the License, or
17  (at your option) any later version.
18 
19  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22  for more details.
23 
24  You should have received a copy of the GNU General Public License
25  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26 
27 \*---------------------------------------------------------------------------*/
28 
29 #include "meshRefinement.H"
30 #include "trackedParticle.H"
31 #include "syncTools.H"
32 #include "Time.H"
33 #include "refinementSurfaces.H"
34 #include "refinementFeatures.H"
35 #include "shellSurfaces.H"
36 #include "faceSet.H"
37 #include "decompositionMethod.H"
38 #include "fvMeshDistribute.H"
39 #include "polyTopoChange.H"
40 #include "mapDistributePolyMesh.H"
41 #include "Cloud.H"
42 //#include "globalIndex.H"
43 #include "OBJstream.H"
44 #include "cellSet.H"
45 #include "treeDataCell.H"
46 
47 #include "cellCuts.H"
48 #include "refineCell.H"
49 #include "hexCellLooper.H"
50 #include "meshCutter.H"
51 
52 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
53 
54 namespace Foam
55 {
56  //- To compare normals
57  static bool less(const vector& x, const vector& y)
58  {
59  for (direction i = 0; i < vector::nComponents; i++)
60  {
61  if (x[i] < y[i])
62  {
63  return true;
64  }
65  else if (x[i] > y[i])
66  {
67  return false;
68  }
69  }
70  // All components the same
71  return false;
72  }
73 
74 
75  //- To compare normals
76  class normalLess
77  {
78  const vectorList& values_;
79 
80  public:
81 
83  :
84  values_(values)
85  {}
86 
87  bool operator()(const label a, const label b) const
88  {
89  return less(values_[a], values_[b]);
90  }
91  };
92 
93 
94  //- Template specialization for pTraits<labelList> so we can have fields
95  template<>
97  {
98 
99  public:
100 
101  //- Component type
103  };
104 
105  //- Template specialization for pTraits<labelList> so we can have fields
106  template<>
108  {
109 
110  public:
111 
112  //- Component type
114  };
115 }
116 
117 
118 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
119 
120 // Get faces (on the new mesh) that have in some way been affected by the
121 // mesh change. Picks up all faces but those that are between two
122 // unrefined faces. (Note that of an unchanged face the edge still might be
123 // split but that does not change any face centre or cell centre.
124 Foam::labelList Foam::meshRefinement::getChangedFaces
125 (
126  const mapPolyMesh& map,
127  const labelList& oldCellsToRefine
128 )
129 {
130  const polyMesh& mesh = map.mesh();
131 
132  labelList changedFaces;
133 
134  // For reporting: number of masterFaces changed
135  label nMasterChanged = 0;
136  bitSet isMasterFace(syncTools::getMasterFaces(mesh));
137 
138  {
139  // Mark any face on a cell which has been added or changed
140  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
141  // Note that refining a face changes the face centre (for a warped face)
142  // which changes the cell centre. This again changes the cellcentre-
143  // cellcentre edge across all faces of the cell.
144  // Note: this does not happen for unwarped faces but unfortunately
145  // we don't have this information.
146 
147  const labelList& faceOwner = mesh.faceOwner();
148  const labelList& faceNeighbour = mesh.faceNeighbour();
149  const cellList& cells = mesh.cells();
150  const label nInternalFaces = mesh.nInternalFaces();
151 
152  // Mark refined cells on old mesh
153  bitSet oldRefineCell(map.nOldCells(), oldCellsToRefine);
154 
155  // Mark refined faces
156  bitSet refinedInternalFace(nInternalFaces);
157 
158  // 1. Internal faces
159 
160  for (label faceI = 0; faceI < nInternalFaces; faceI++)
161  {
162  label oldOwn = map.cellMap()[faceOwner[faceI]];
163  label oldNei = map.cellMap()[faceNeighbour[faceI]];
164 
165  if
166  (
167  oldOwn >= 0 && !oldRefineCell.test(oldOwn)
168  && oldNei >= 0 && !oldRefineCell.test(oldNei)
169  )
170  {
171  // Unaffected face since both neighbours were not refined.
172  }
173  else
174  {
175  refinedInternalFace.set(faceI);
176  }
177  }
178 
179 
180  // 2. Boundary faces
181 
182  boolList refinedBoundaryFace(mesh.nBoundaryFaces(), false);
183 
184  forAll(mesh.boundaryMesh(), patchI)
185  {
186  const polyPatch& pp = mesh.boundaryMesh()[patchI];
187 
188  label faceI = pp.start();
189 
190  forAll(pp, i)
191  {
192  label oldOwn = map.cellMap()[faceOwner[faceI]];
193 
194  if (oldOwn >= 0 && !oldRefineCell.test(oldOwn))
195  {
196  // owner did exist and wasn't refined.
197  }
198  else
199  {
200  refinedBoundaryFace[faceI-nInternalFaces] = true;
201  }
202  faceI++;
203  }
204  }
205 
206  // Synchronise refined face status
208  (
209  mesh,
210  refinedBoundaryFace,
211  orEqOp<bool>()
212  );
213 
214 
215  // 3. Mark all faces affected by refinement. Refined faces are in
216  // - refinedInternalFace
217  // - refinedBoundaryFace
218  boolList changedFace(mesh.nFaces(), false);
219 
220  forAll(refinedInternalFace, faceI)
221  {
222  if (refinedInternalFace.test(faceI))
223  {
224  const cell& ownFaces = cells[faceOwner[faceI]];
225  forAll(ownFaces, ownI)
226  {
227  changedFace[ownFaces[ownI]] = true;
228  }
229  const cell& neiFaces = cells[faceNeighbour[faceI]];
230  forAll(neiFaces, neiI)
231  {
232  changedFace[neiFaces[neiI]] = true;
233  }
234  }
235  }
236 
237  forAll(refinedBoundaryFace, i)
238  {
239  if (refinedBoundaryFace[i])
240  {
241  const cell& ownFaces = cells[faceOwner[i+nInternalFaces]];
242  forAll(ownFaces, ownI)
243  {
244  changedFace[ownFaces[ownI]] = true;
245  }
246  }
247  }
248 
250  (
251  mesh,
252  changedFace,
253  orEqOp<bool>()
254  );
255 
256 
257  // Now we have in changedFace marked all affected faces. Pack.
258  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
259 
260  changedFaces = findIndices(changedFace, true);
261 
262  // Count changed master faces.
263  nMasterChanged = 0;
264 
265  forAll(changedFace, faceI)
266  {
267  if (changedFace[faceI] && isMasterFace.test(faceI))
268  {
269  nMasterChanged++;
270  }
271  }
272 
273  }
274 
276  {
277  Pout<< "getChangedFaces : Detected "
278  << " local:" << changedFaces.size()
279  << " global:" << returnReduce(nMasterChanged, sumOp<label>())
280  << " changed faces out of " << mesh.globalData().nTotalFaces()
281  << endl;
282 
283  faceSet changedFacesSet(mesh, "changedFaces", changedFaces);
284  Pout<< "getChangedFaces : Writing " << changedFaces.size()
285  << " changed faces to faceSet " << changedFacesSet.name()
286  << endl;
287  changedFacesSet.write();
288  }
289 
290  return changedFaces;
291 }
292 
293 
294 // Mark cell for refinement (if not already marked). Return false if
295 // refinelimit hit. Keeps running count (in nRefine) of cells marked for
296 // refinement
297 bool Foam::meshRefinement::markForRefine
298 (
299  const label markValue,
300  const label nAllowRefine,
301 
302  label& cellValue,
303  label& nRefine
304 )
305 {
306  if (cellValue == -1)
307  {
308  cellValue = markValue;
309  nRefine++;
310  }
311 
312  return nRefine <= nAllowRefine;
313 }
314 
315 
316 void Foam::meshRefinement::markFeatureCellLevel
317 (
318  const pointField& keepPoints,
319  labelList& maxFeatureLevel
320 ) const
321 {
322  // We want to refine all cells containing a feature edge.
323  // - don't want to search over whole mesh
324  // - don't want to build octree for whole mesh
325  // - so use tracking to follow the feature edges
326  //
327  // 1. find non-manifold points on feature edges (i.e. start of feature edge
328  // or 'knot')
329  // 2. seed particle starting at keepPoint going to this non-manifold point
330  // 3. track particles to their non-manifold point
331  //
332  // 4. track particles across their connected feature edges, marking all
333  // visited cells with their level (through trackingData)
334  // 5. do 4 until all edges have been visited.
335 
336 
337  // Find all start cells of features. Is done by tracking from keepPoint.
338  Cloud<trackedParticle> startPointCloud
339  (
340  mesh_,
341  "startPointCloud",
342  IDLList<trackedParticle>()
343  );
344 
345 
346  // Features are identical on all processors. Number them so we know
347  // what to seed. Do this on only the processor that
348  // holds the keepPoint.
349 
350  for (const point& keepPoint : keepPoints)
351  {
352  const label celli = mesh_.cellTree().findInside(keepPoint);
353 
354  if (celli != -1)
355  {
356  // I am the processor that holds the keepPoint
357 
358  forAll(features_, feati)
359  {
360  const edgeMesh& featureMesh = features_[feati];
361  const label featureLevel = features_.levels()[feati][0];
362  const labelListList& pointEdges = featureMesh.pointEdges();
363 
364  // Find regions on edgeMesh
365  labelList edgeRegion;
366  label nRegions = featureMesh.regions(edgeRegion);
367 
368 
369  bitSet regionVisited(nRegions);
370 
371 
372  // 1. Seed all 'knots' in edgeMesh
373 
374 
375  forAll(pointEdges, pointi)
376  {
377  if (pointEdges[pointi].size() != 2)
378  {
380  {
381  Pout<< "Adding particle from point:" << pointi
382  << " coord:" << featureMesh.points()[pointi]
383  << " since number of emanating edges:"
384  << pointEdges[pointi].size()
385  << endl;
386  }
387 
388  // Non-manifold point. Create particle.
389  startPointCloud.addParticle
390  (
391  new trackedParticle
392  (
393  mesh_,
394  keepPoint,
395  celli,
396  featureMesh.points()[pointi], // endpos
397  featureLevel, // level
398  feati, // featureMesh
399  pointi, // end point
400  -1 // feature edge
401  )
402  );
403 
404  // Mark
405  if (pointEdges[pointi].size() > 0)
406  {
407  label e0 = pointEdges[pointi][0];
408  label regioni = edgeRegion[e0];
409  regionVisited.set(regioni);
410  }
411  }
412  }
413 
414 
415  // 2. Any regions that have not been visited at all? These can
416  // only be circular regions!
417  forAll(featureMesh.edges(), edgei)
418  {
419  if (regionVisited.set(edgeRegion[edgei]))
420  {
421  const edge& e = featureMesh.edges()[edgei];
422  label pointi = e.start();
424  {
425  Pout<< "Adding particle from point:" << pointi
426  << " coord:" << featureMesh.points()[pointi]
427  << " on circular region:" << edgeRegion[edgei]
428  << endl;
429  }
430 
431  // Non-manifold point. Create particle.
432  startPointCloud.addParticle
433  (
434  new trackedParticle
435  (
436  mesh_,
437  keepPoint,
438  celli,
439  featureMesh.points()[pointi], // endpos
440  featureLevel, // level
441  feati, // featureMesh
442  pointi, // end point
443  -1 // feature edge
444  )
445  );
446  }
447  }
448  }
449  }
450  }
451 
452 
453  // Largest refinement level of any feature passed through
454  maxFeatureLevel = labelList(mesh_.nCells(), -1);
455 
456  // Whether edge has been visited.
457  List<bitSet> featureEdgeVisited(features_.size());
458 
459  forAll(features_, featI)
460  {
461  featureEdgeVisited[featI].setSize(features_[featI].edges().size());
462  featureEdgeVisited[featI] = false;
463  }
464 
465  // Database to pass into trackedParticle::move
466  trackedParticle::trackingData td
467  (
468  startPointCloud,
469  maxFeatureLevel,
470  featureEdgeVisited
471  );
472 
473 
474  // Track all particles to their end position (= starting feature point)
475  // Note that the particle might have started on a different processor
476  // so this will transfer across nicely until we can start tracking proper.
477  scalar maxTrackLen = 2.0*mesh_.bounds().mag();
478 
480  {
481  Pout<< "Tracking " << startPointCloud.size()
482  << " particles over distance " << maxTrackLen
483  << " to find the starting cell" << endl;
484  }
485  startPointCloud.move(startPointCloud, td, maxTrackLen);
486 
487 
488  // Reset levels
489  maxFeatureLevel = -1;
490  forAll(features_, featI)
491  {
492  featureEdgeVisited[featI] = false;
493  }
494 
495 
496  Cloud<trackedParticle> cloud
497  (
498  mesh_,
499  "featureCloud",
500  IDLList<trackedParticle>()
501  );
502 
504  {
505  Pout<< "Constructing cloud for cell marking" << endl;
506  }
507 
508  for (const trackedParticle& startTp : startPointCloud)
509  {
510  const label featI = startTp.i();
511  const label pointI = startTp.j();
512 
513  const edgeMesh& featureMesh = features_[featI];
514  const labelList& pEdges = featureMesh.pointEdges()[pointI];
515 
516  // Now shoot particles down all pEdges.
517  forAll(pEdges, pEdgeI)
518  {
519  label edgeI = pEdges[pEdgeI];
520 
521  if (featureEdgeVisited[featI].set(edgeI))
522  {
523  // Unvisited edge. Make the particle go to the other point
524  // on the edge.
525 
526  const edge& e = featureMesh.edges()[edgeI];
527  label otherPointi = e.otherVertex(pointI);
528 
529  trackedParticle* tp(new trackedParticle(startTp));
530  tp->start() = tp->position();
531  tp->end() = featureMesh.points()[otherPointi];
532  tp->j() = otherPointi;
533  tp->k() = edgeI;
534 
536  {
537  Pout<< "Adding particle for point:" << pointI
538  << " coord:" << tp->position()
539  << " feature:" << featI
540  << " to track to:" << tp->end()
541  << endl;
542  }
543 
544  cloud.addParticle(tp);
545  }
546  }
547  }
548 
549  startPointCloud.clear();
550 
551 
552  while (true)
553  {
554  // Track all particles to their end position.
556  {
557  Pout<< "Tracking " << cloud.size()
558  << " particles over distance " << maxTrackLen
559  << " to mark cells" << endl;
560  }
561  cloud.move(cloud, td, maxTrackLen);
562 
563  // Make particle follow edge.
564  for (trackedParticle& tp : cloud)
565  {
566  const label featI = tp.i();
567  const label pointI = tp.j();
568 
569  const edgeMesh& featureMesh = features_[featI];
570  const labelList& pEdges = featureMesh.pointEdges()[pointI];
571 
572  // Particle now at pointI. Check connected edges to see which one
573  // we have to visit now.
574 
575  bool keepParticle = false;
576 
577  forAll(pEdges, i)
578  {
579  label edgeI = pEdges[i];
580 
581  if (featureEdgeVisited[featI].set(edgeI))
582  {
583  // Unvisited edge. Make the particle go to the other point
584  // on the edge.
585 
586  const edge& e = featureMesh.edges()[edgeI];
587  label otherPointi = e.otherVertex(pointI);
588 
589  tp.start() = tp.position();
590  tp.end() = featureMesh.points()[otherPointi];
591  tp.j() = otherPointi;
592  tp.k() = edgeI;
593  keepParticle = true;
594  break;
595  }
596  }
597 
598  if (!keepParticle)
599  {
600  // Particle at 'knot' where another particle already has been
601  // seeded. Delete particle.
602  cloud.deleteParticle(tp);
603  }
604  }
605 
606 
608  {
609  Pout<< "Remaining particles " << cloud.size() << endl;
610  }
611 
612  if (returnReduce(cloud.size(), sumOp<label>()) == 0)
613  {
614  break;
615  }
616  }
617 
618 
619 
620  //if (debug&meshRefinement::FEATURESEEDS)
621  //{
622  // forAll(maxFeatureLevel, cellI)
623  // {
624  // if (maxFeatureLevel[cellI] != -1)
625  // {
626  // Pout<< "Feature went through cell:" << cellI
627  // << " coord:" << mesh_.cellCentres()[cellI]
628  // << " level:" << maxFeatureLevel[cellI]
629  // << endl;
630  // }
631  // }
632  //}
633 }
634 
635 
636 // Calculates list of cells to refine based on intersection with feature edge.
637 Foam::label Foam::meshRefinement::markFeatureRefinement
638 (
639  const pointField& keepPoints,
640  const label nAllowRefine,
641 
642  labelList& refineCell,
643  label& nRefine
644 ) const
645 {
646  // Largest refinement level of any feature passed through
647  labelList maxFeatureLevel;
648  markFeatureCellLevel(keepPoints, maxFeatureLevel);
649 
650  // See which cells to refine. maxFeatureLevel will hold highest level
651  // of any feature edge that passed through.
652 
653  const labelList& cellLevel = meshCutter_.cellLevel();
654 
655  label oldNRefine = nRefine;
656 
657  forAll(maxFeatureLevel, cellI)
658  {
659  if (maxFeatureLevel[cellI] > cellLevel[cellI])
660  {
661  // Mark
662  if
663  (
664  !markForRefine
665  (
666  0, // surface (n/a)
667  nAllowRefine,
668  refineCell[cellI],
669  nRefine
670  )
671  )
672  {
673  // Reached limit
674  break;
675  }
676  }
677  }
678 
679  if
680  (
681  returnReduce(nRefine, sumOp<label>())
682  > returnReduce(nAllowRefine, sumOp<label>())
683  )
684  {
685  Info<< "Reached refinement limit." << endl;
686  }
687 
688  return returnReduce(nRefine-oldNRefine, sumOp<label>());
689 }
690 
691 
692 // Mark cells for distance-to-feature based refinement.
693 Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
694 (
695  const label nAllowRefine,
696 
697  labelList& refineCell,
698  label& nRefine
699 ) const
700 {
701  const labelList& cellLevel = meshCutter_.cellLevel();
702  const pointField& cellCentres = mesh_.cellCentres();
703 
704  // Detect if there are any distance shells
705  if (features_.maxDistance() <= 0.0)
706  {
707  return 0;
708  }
709 
710  label oldNRefine = nRefine;
711 
712  // Collect cells to test
713  pointField testCc(cellLevel.size()-nRefine);
714  labelList testLevels(cellLevel.size()-nRefine);
715  label testI = 0;
716 
717  forAll(cellLevel, cellI)
718  {
719  if (refineCell[cellI] == -1)
720  {
721  testCc[testI] = cellCentres[cellI];
722  testLevels[testI] = cellLevel[cellI];
723  testI++;
724  }
725  }
726 
727  // Do test to see whether cells are inside/outside shell with higher level
728  labelList maxLevel;
729  features_.findHigherLevel(testCc, testLevels, maxLevel);
730 
731  // Mark for refinement. Note that we didn't store the original cellID so
732  // now just reloop in same order.
733  testI = 0;
734  forAll(cellLevel, cellI)
735  {
736  if (refineCell[cellI] == -1)
737  {
738  if (maxLevel[testI] > testLevels[testI])
739  {
740  bool reachedLimit = !markForRefine
741  (
742  maxLevel[testI], // mark with any positive value
743  nAllowRefine,
744  refineCell[cellI],
745  nRefine
746  );
747 
748  if (reachedLimit)
749  {
750  if (debug)
751  {
752  Pout<< "Stopped refining internal cells"
753  << " since reaching my cell limit of "
754  << mesh_.nCells()+7*nRefine << endl;
755  }
756  break;
757  }
758  }
759  testI++;
760  }
761  }
762 
763  if
764  (
765  returnReduce(nRefine, sumOp<label>())
766  > returnReduce(nAllowRefine, sumOp<label>())
767  )
768  {
769  Info<< "Reached refinement limit." << endl;
770  }
771 
772  return returnReduce(nRefine-oldNRefine, sumOp<label>());
773 }
774 
775 
776 // Mark cells for non-surface intersection based refinement.
777 Foam::label Foam::meshRefinement::markInternalRefinement
778 (
779  const label nAllowRefine,
780 
781  labelList& refineCell,
782  label& nRefine
783 ) const
784 {
785  const labelList& cellLevel = meshCutter_.cellLevel();
786  const pointField& cellCentres = mesh_.cellCentres();
787 
788  label oldNRefine = nRefine;
789 
790  // Collect cells to test
791  pointField testCc(cellLevel.size()-nRefine);
792  labelList testLevels(cellLevel.size()-nRefine);
793  label testI = 0;
794 
795  forAll(cellLevel, cellI)
796  {
797  if (refineCell[cellI] == -1)
798  {
799  testCc[testI] = cellCentres[cellI];
800  testLevels[testI] = cellLevel[cellI];
801  testI++;
802  }
803  }
804 
805  // Do test to see whether cells are inside/outside shell with higher level
806  labelList maxLevel;
807  shells_.findHigherLevel(testCc, testLevels, maxLevel);
808 
809  // Mark for refinement. Note that we didn't store the original cellID so
810  // now just reloop in same order.
811  testI = 0;
812  forAll(cellLevel, cellI)
813  {
814  if (refineCell[cellI] == -1)
815  {
816  if (maxLevel[testI] > testLevels[testI])
817  {
818  bool reachedLimit = !markForRefine
819  (
820  maxLevel[testI], // mark with any positive value
821  nAllowRefine,
822  refineCell[cellI],
823  nRefine
824  );
825 
826  if (reachedLimit)
827  {
828  if (debug)
829  {
830  Pout<< "Stopped refining internal cells"
831  << " since reaching my cell limit of "
832  << mesh_.nCells()+7*nRefine << endl;
833  }
834  break;
835  }
836  }
837  testI++;
838  }
839  }
840 
841  if
842  (
843  returnReduce(nRefine, sumOp<label>())
844  > returnReduce(nAllowRefine, sumOp<label>())
845  )
846  {
847  Info<< "Reached refinement limit." << endl;
848  }
849 
850  return returnReduce(nRefine-oldNRefine, sumOp<label>());
851 }
852 
853 
854 Foam::label Foam::meshRefinement::unmarkInternalRefinement
855 (
856  labelList& refineCell,
857  label& nRefine
858 ) const
859 {
860  const labelList& cellLevel = meshCutter_.cellLevel();
861  const pointField& cellCentres = mesh_.cellCentres();
862 
863  label oldNRefine = nRefine;
864 
865  // Collect cells to test
866  pointField testCc(nRefine);
867  labelList testLevels(nRefine);
868  label testI = 0;
869 
870  forAll(cellLevel, cellI)
871  {
872  if (refineCell[cellI] >= 0)
873  {
874  testCc[testI] = cellCentres[cellI];
875  testLevels[testI] = cellLevel[cellI];
876  testI++;
877  }
878  }
879 
880  // Do test to see whether cells are inside/outside shell with higher level
881  labelList levelShell;
882  limitShells_.findLevel(testCc, testLevels, levelShell);
883 
884  // Mark for refinement. Note that we didn't store the original cellID so
885  // now just reloop in same order.
886  testI = 0;
887  forAll(cellLevel, cellI)
888  {
889  if (refineCell[cellI] >= 0)
890  {
891  if (levelShell[testI] != -1)
892  {
893  refineCell[cellI] = -1;
894  nRefine--;
895  }
896  testI++;
897  }
898  }
899 
900  return returnReduce(oldNRefine-nRefine, sumOp<label>());
901 }
902 
903 
904 // Collect faces that are intersected and whose neighbours aren't yet marked
905 // for refinement.
906 Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
907 (
908  const labelList& refineCell
909 ) const
910 {
911  labelList testFaces(mesh_.nFaces());
912 
913  label nTest = 0;
914 
915  const labelList& surfIndex = surfaceIndex();
916 
917  forAll(surfIndex, faceI)
918  {
919  if (surfIndex[faceI] != -1)
920  {
921  label own = mesh_.faceOwner()[faceI];
922 
923  if (mesh_.isInternalFace(faceI))
924  {
925  label nei = mesh_.faceNeighbour()[faceI];
926 
927  if (refineCell[own] == -1 || refineCell[nei] == -1)
928  {
929  testFaces[nTest++] = faceI;
930  }
931  }
932  else
933  {
934  if (refineCell[own] == -1)
935  {
936  testFaces[nTest++] = faceI;
937  }
938  }
939  }
940  }
941  testFaces.setSize(nTest);
942 
943  return testFaces;
944 }
945 
946 
947 // Mark cells for surface intersection based refinement.
948 Foam::label Foam::meshRefinement::markSurfaceRefinement
949 (
950  const label nAllowRefine,
951  const labelList& neiLevel,
952  const pointField& neiCc,
953 
954  labelList& refineCell,
955  label& nRefine
956 ) const
957 {
958  const labelList& cellLevel = meshCutter_.cellLevel();
959 
960  label oldNRefine = nRefine;
961 
962  // Use cached surfaceIndex_ to detect if any intersection. If so
963  // re-intersect to determine level wanted.
964 
965  // Collect candidate faces
966  // ~~~~~~~~~~~~~~~~~~~~~~~
967 
968  labelList testFaces(getRefineCandidateFaces(refineCell));
969 
970  // Collect segments
971  // ~~~~~~~~~~~~~~~~
972 
973  pointField start(testFaces.size());
974  pointField end(testFaces.size());
975  labelList minLevel(testFaces.size());
976 
977  calcCellCellRays
978  (
979  neiCc,
980  neiLevel,
981  testFaces,
982  start,
983  end,
984  minLevel
985  );
986 
987 
988  // Do test for higher intersections
989  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
990 
991  labelList surfaceHit;
992  labelList surfaceMinLevel;
993  surfaces_.findHigherIntersection
994  (
995  shells_,
996  start,
997  end,
998  minLevel,
999 
1000  surfaceHit,
1001  surfaceMinLevel
1002  );
1003 
1004 
1005  // Mark cells for refinement
1006  // ~~~~~~~~~~~~~~~~~~~~~~~~~
1007 
1008  forAll(testFaces, i)
1009  {
1010  label faceI = testFaces[i];
1011 
1012  label surfI = surfaceHit[i];
1013 
1014  if (surfI != -1)
1015  {
1016  // Found intersection with surface with higher wanted
1017  // refinement. Check if the level field on that surface
1018  // specifies an even higher level. Note:this is weird. Should
1019  // do the check with the surfaceMinLevel whilst intersecting the
1020  // surfaces?
1021 
1022  label own = mesh_.faceOwner()[faceI];
1023 
1024  if (surfaceMinLevel[i] > cellLevel[own])
1025  {
1026  // Owner needs refining
1027  if
1028  (
1029  !markForRefine
1030  (
1031  surfI,
1032  nAllowRefine,
1033  refineCell[own],
1034  nRefine
1035  )
1036  )
1037  {
1038  break;
1039  }
1040  }
1041 
1042  if (mesh_.isInternalFace(faceI))
1043  {
1044  label nei = mesh_.faceNeighbour()[faceI];
1045  if (surfaceMinLevel[i] > cellLevel[nei])
1046  {
1047  // Neighbour needs refining
1048  if
1049  (
1050  !markForRefine
1051  (
1052  surfI,
1053  nAllowRefine,
1054  refineCell[nei],
1055  nRefine
1056  )
1057  )
1058  {
1059  break;
1060  }
1061  }
1062  }
1063  }
1064  }
1065 
1066  if
1067  (
1068  returnReduce(nRefine, sumOp<label>())
1069  > returnReduce(nAllowRefine, sumOp<label>())
1070  )
1071  {
1072  Info<< "Reached refinement limit." << endl;
1073  }
1074 
1075  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1076 }
1077 
1078 
1079 // Count number of matches of first argument in second argument
1080 Foam::label Foam::meshRefinement::countMatches
1081 (
1082  const List<point>& normals1,
1083  const List<point>& normals2,
1084  const scalar tol
1085 ) const
1086 {
1087  label nMatches = 0;
1088 
1089  forAll(normals1, i)
1090  {
1091  const vector& n1 = normals1[i];
1092 
1093  forAll(normals2, j)
1094  {
1095  const vector& n2 = normals2[j];
1096 
1097  if (magSqr(n1-n2) < tol)
1098  {
1099  nMatches++;
1100  break;
1101  }
1102  }
1103  }
1104 
1105  return nMatches;
1106 }
1107 
1108 
1109 // Mark cells for surface curvature based refinement.
1110 Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
1111 (
1112  const scalar curvature,
1113  const label nAllowRefine,
1114  const labelList& neiLevel,
1115  const pointField& neiCc,
1116 
1117  labelList& refineCell,
1118  label& nRefine
1119 ) const
1120 {
1121  const labelList& cellLevel = meshCutter_.cellLevel();
1122  const pointField& cellCentres = mesh_.cellCentres();
1123 
1124  label oldNRefine = nRefine;
1125 
1126  // 1. local test: any cell on more than one surface gets refined
1127  // (if its current level is < max of the surface max level)
1128 
1129  // 2. 'global' test: any cell on only one surface with a neighbour
1130  // on a different surface gets refined (if its current level etc.)
1131 
1132 
1133  const bitSet isMasterFace(syncTools::getMasterFaces(mesh_));
1134 
1135 
1136  // Collect candidate faces (i.e. intersecting any surface and
1137  // owner/neighbour not yet refined.
1138  labelList testFaces(getRefineCandidateFaces(refineCell));
1139 
1140  // Collect segments
1141  pointField start(testFaces.size());
1142  pointField end(testFaces.size());
1143  labelList minLevel(testFaces.size());
1144 
1145  // Note: uses isMasterFace otherwise could be call to calcCellCellRays
1146  forAll(testFaces, i)
1147  {
1148  label faceI = testFaces[i];
1149 
1150  label own = mesh_.faceOwner()[faceI];
1151 
1152  if (mesh_.isInternalFace(faceI))
1153  {
1154  label nei = mesh_.faceNeighbour()[faceI];
1155 
1156  start[i] = cellCentres[own];
1157  end[i] = cellCentres[nei];
1158  minLevel[i] = min(cellLevel[own], cellLevel[nei]);
1159  }
1160  else
1161  {
1162  label bFaceI = faceI - mesh_.nInternalFaces();
1163 
1164  start[i] = cellCentres[own];
1165  end[i] = neiCc[bFaceI];
1166 
1167  if (!isMasterFace[faceI])
1168  {
1169  std::swap(start[i], end[i]);
1170  }
1171 
1172  minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
1173  }
1174  }
1175 
1176  // Extend segments a bit
1177  {
1178  const vectorField smallVec(ROOTSMALL*(end-start));
1179  start -= smallVec;
1180  end += smallVec;
1181  }
1182 
1183 
1184  // Test for all intersections (with surfaces of higher max level than
1185  // minLevel) and cache per cell the interesting inter
1186  labelListList cellSurfLevels(mesh_.nCells());
1187  List<vectorList> cellSurfNormals(mesh_.nCells());
1188 
1189  {
1190  // Per segment the normals of the surfaces
1191  List<vectorList> surfaceNormal;
1192  // Per segment the list of levels of the surfaces
1193  labelListList surfaceLevel;
1194 
1195  surfaces_.findAllIntersections
1196  (
1197  start,
1198  end,
1199  minLevel, // max level of surface has to be bigger
1200  // than min level of neighbouring cells
1201 
1202  labelList(surfaces_.maxLevel().size(), 0), // min level
1203  surfaces_.maxLevel(),
1204 
1205  surfaceNormal,
1206  surfaceLevel
1207  );
1208 
1209 
1210  // Sort the data according to surface location. This will guarantee
1211  // that on coupled faces both sides visit the intersections in
1212  // the same order so will decide the same
1213  labelList visitOrder;
1214  forAll(surfaceNormal, pointI)
1215  {
1216  vectorList& pNormals = surfaceNormal[pointI];
1217  labelList& pLevel = surfaceLevel[pointI];
1218 
1219  sortedOrder(pNormals, visitOrder, normalLess(pNormals));
1220 
1221  pNormals = List<point>(pNormals, visitOrder);
1222  pLevel = labelUIndList(pLevel, visitOrder);
1223  }
1224 
1225  // Clear out unnecessary data
1226  start.clear();
1227  end.clear();
1228  minLevel.clear();
1229 
1230  // Convert face-wise data to cell.
1231  forAll(testFaces, i)
1232  {
1233  label faceI = testFaces[i];
1234  label own = mesh_.faceOwner()[faceI];
1235 
1236  const vectorList& fNormals = surfaceNormal[i];
1237  const labelList& fLevels = surfaceLevel[i];
1238 
1239  forAll(fNormals, hitI)
1240  {
1241  if (fLevels[hitI] > cellLevel[own])
1242  {
1243  cellSurfLevels[own].append(fLevels[hitI]);
1244  cellSurfNormals[own].append(fNormals[hitI]);
1245  }
1246 
1247  if (mesh_.isInternalFace(faceI))
1248  {
1249  label nei = mesh_.faceNeighbour()[faceI];
1250  if (fLevels[hitI] > cellLevel[nei])
1251  {
1252  cellSurfLevels[nei].append(fLevels[hitI]);
1253  cellSurfNormals[nei].append(fNormals[hitI]);
1254  }
1255  }
1256  }
1257  }
1258  }
1259 
1260 
1261 
1262  // Bit of statistics
1263  if (debug)
1264  {
1265  label nSet = 0;
1266  label nNormals = 9;
1267  forAll(cellSurfNormals, cellI)
1268  {
1269  const vectorList& normals = cellSurfNormals[cellI];
1270  if (normals.size())
1271  {
1272  nSet++;
1273  nNormals += normals.size();
1274  }
1275  }
1276  reduce(nSet, sumOp<label>());
1277  reduce(nNormals, sumOp<label>());
1278  Info<< "markSurfaceCurvatureRefinement :"
1279  << " cells:" << mesh_.globalData().nTotalCells()
1280  << " of which with normals:" << nSet
1281  << " ; total normals stored:" << nNormals
1282  << endl;
1283  }
1284 
1285 
1286 
1287  bool reachedLimit = false;
1288 
1289 
1290  // 1. Check normals per cell
1291  // ~~~~~~~~~~~~~~~~~~~~~~~~~
1292 
1293  for
1294  (
1295  label cellI = 0;
1296  !reachedLimit && cellI < cellSurfNormals.size();
1297  cellI++
1298  )
1299  {
1300  const vectorList& normals = cellSurfNormals[cellI];
1301  const labelList& levels = cellSurfLevels[cellI];
1302 
1303  // n^2 comparison of all normals in a cell
1304  for (label i = 0; !reachedLimit && i < normals.size(); i++)
1305  {
1306  for (label j = i+1; !reachedLimit && j < normals.size(); j++)
1307  {
1308  if ((normals[i] & normals[j]) < curvature)
1309  {
1310  label maxLevel = max(levels[i], levels[j]);
1311 
1312  if (cellLevel[cellI] < maxLevel)
1313  {
1314  if
1315  (
1316  !markForRefine
1317  (
1318  maxLevel,
1319  nAllowRefine,
1320  refineCell[cellI],
1321  nRefine
1322  )
1323  )
1324  {
1325  if (debug)
1326  {
1327  Pout<< "Stopped refining since reaching my cell"
1328  << " limit of " << mesh_.nCells()+7*nRefine
1329  << endl;
1330  }
1331  reachedLimit = true;
1332  break;
1333  }
1334  }
1335  }
1336  }
1337  }
1338  }
1339 
1340 
1341 
1342  // 2. Find out a measure of surface curvature
1343  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1344  // Look at normals between neighbouring surfaces
1345  // Loop over all faces. Could only be checkFaces, except if they're coupled
1346 
1347  // Internal faces
1348  for
1349  (
1350  label faceI = 0;
1351  !reachedLimit && faceI < mesh_.nInternalFaces();
1352  faceI++
1353  )
1354  {
1355  label own = mesh_.faceOwner()[faceI];
1356  label nei = mesh_.faceNeighbour()[faceI];
1357 
1358  const vectorList& ownNormals = cellSurfNormals[own];
1359  const labelList& ownLevels = cellSurfLevels[own];
1360  const vectorList& neiNormals = cellSurfNormals[nei];
1361  const labelList& neiLevels = cellSurfLevels[nei];
1362 
1363 
1364  // Special case: owner normals set is a subset of the neighbour
1365  // normals. Do not do curvature refinement since other cell's normals
1366  // do not add any information. Avoids weird corner extensions of
1367  // refinement regions.
1368  bool ownIsSubset =
1369  countMatches(ownNormals, neiNormals)
1370  == ownNormals.size();
1371 
1372  bool neiIsSubset =
1373  countMatches(neiNormals, ownNormals)
1374  == neiNormals.size();
1375 
1376 
1377  if (!ownIsSubset && !neiIsSubset)
1378  {
1379  // n^2 comparison of between ownNormals and neiNormals
1380  for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
1381  {
1382  for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
1383  {
1384  // Have valid data on both sides. Check curvature.
1385  if ((ownNormals[i] & neiNormals[j]) < curvature)
1386  {
1387  // See which side to refine.
1388  if (cellLevel[own] < ownLevels[i])
1389  {
1390  if
1391  (
1392  !markForRefine
1393  (
1394  ownLevels[i],
1395  nAllowRefine,
1396  refineCell[own],
1397  nRefine
1398  )
1399  )
1400  {
1401  if (debug)
1402  {
1403  Pout<< "Stopped refining since reaching"
1404  << " my cell limit of "
1405  << mesh_.nCells()+7*nRefine << endl;
1406  }
1407  reachedLimit = true;
1408  break;
1409  }
1410  }
1411  if (cellLevel[nei] < neiLevels[j])
1412  {
1413  if
1414  (
1415  !markForRefine
1416  (
1417  neiLevels[j],
1418  nAllowRefine,
1419  refineCell[nei],
1420  nRefine
1421  )
1422  )
1423  {
1424  if (debug)
1425  {
1426  Pout<< "Stopped refining since reaching"
1427  << " my cell limit of "
1428  << mesh_.nCells()+7*nRefine << endl;
1429  }
1430  reachedLimit = true;
1431  break;
1432  }
1433  }
1434  }
1435  }
1436  }
1437  }
1438  }
1439 
1440 
1441  // Send over surface normal to neighbour cell.
1442  List<vectorList> neiSurfaceNormals;
1443  syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
1444 
1445  // Boundary faces
1446  for
1447  (
1448  label faceI = mesh_.nInternalFaces();
1449  !reachedLimit && faceI < mesh_.nFaces();
1450  faceI++
1451  )
1452  {
1453  label own = mesh_.faceOwner()[faceI];
1454  label bFaceI = faceI - mesh_.nInternalFaces();
1455 
1456  const vectorList& ownNormals = cellSurfNormals[own];
1457  const labelList& ownLevels = cellSurfLevels[own];
1458  const vectorList& neiNormals = neiSurfaceNormals[bFaceI];
1459 
1460  // Special case: owner normals set is a subset of the neighbour
1461  // normals. Do not do curvature refinement since other cell's normals
1462  // do not add any information. Avoids weird corner extensions of
1463  // refinement regions.
1464  bool ownIsSubset =
1465  countMatches(ownNormals, neiNormals)
1466  == ownNormals.size();
1467 
1468  bool neiIsSubset =
1469  countMatches(neiNormals, ownNormals)
1470  == neiNormals.size();
1471 
1472 
1473  if (!ownIsSubset && !neiIsSubset)
1474  {
1475  // n^2 comparison of between ownNormals and neiNormals
1476  for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
1477  {
1478  for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
1479  {
1480  // Have valid data on both sides. Check curvature.
1481  if ((ownNormals[i] & neiNormals[j]) < curvature)
1482  {
1483  if (cellLevel[own] < ownLevels[i])
1484  {
1485  if
1486  (
1487  !markForRefine
1488  (
1489  ownLevels[i],
1490  nAllowRefine,
1491  refineCell[own],
1492  nRefine
1493  )
1494  )
1495  {
1496  if (debug)
1497  {
1498  Pout<< "Stopped refining since reaching"
1499  << " my cell limit of "
1500  << mesh_.nCells()+7*nRefine
1501  << endl;
1502  }
1503  reachedLimit = true;
1504  break;
1505  }
1506  }
1507  }
1508  }
1509  }
1510  }
1511  }
1512 
1513 
1514  if
1515  (
1516  returnReduce(nRefine, sumOp<label>())
1517  > returnReduce(nAllowRefine, sumOp<label>())
1518  )
1519  {
1520  Info<< "Reached refinement limit." << endl;
1521  }
1522 
1523  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1524 }
1525 
1526 
1529  const scalar planarCos,
1530  const vector& point0,
1531  const vector& normal0,
1532 
1533  const vector& point1,
1534  const vector& normal1
1535 ) const
1536 {
1537  //- Hits differ and angles are oppositeish and
1538  // hits have a normal distance
1539  vector d = point1-point0;
1540  scalar magD = mag(d);
1541 
1542  if (magD > mergeDistance())
1543  {
1544  scalar cosAngle = (normal0 & normal1);
1545 
1546  vector avg = Zero;
1547  if (cosAngle < (-1+planarCos))
1548  {
1549  // Opposite normals
1550  avg = 0.5*(normal0-normal1);
1551  }
1552  else if (cosAngle > (1-planarCos))
1553  {
1554  avg = 0.5*(normal0+normal1);
1555  }
1556 
1557  if (avg != vector::zero)
1558  {
1559  avg /= mag(avg);
1560 
1561  // Check normal distance of intersection locations
1562  if (mag(avg&d) > mergeDistance())
1563  {
1564  return true;
1565  }
1566  }
1567  }
1568 
1569  return false;
1570 }
1571 
1572 
1573 // Mark small gaps
1576  const scalar planarCos,
1577  const label level0,
1578  const vector& point0,
1579  const vector& normal0,
1580 
1581  const label level1,
1582  const vector& point1,
1583  const vector& normal1
1584 ) const
1585 {
1586  //const scalar edge0Len = meshCutter_.level0EdgeLength();
1587 
1588  //- Hits differ and angles are oppositeish and
1589  // hits have a normal distance
1590  vector d = point1-point0;
1591  scalar magD = mag(d);
1592 
1593  if (magD > mergeDistance())
1594  {
1595  scalar cosAngle = (normal0 & normal1);
1596 
1597  vector avgNormal = Zero;
1598  if (cosAngle < (-1+planarCos))
1599  {
1600  // Opposite normals
1601  avgNormal = 0.5*(normal0-normal1);
1602  }
1603  else if (cosAngle > (1-planarCos))
1604  {
1605  avgNormal = 0.5*(normal0+normal1);
1606  }
1607 
1608  if (avgNormal != vector::zero)
1609  {
1610  avgNormal /= mag(avgNormal);
1611 
1612  //scalar normalDist = mag(d&avgNormal);
1613  //const scalar avgCellSize = edge0Len/pow(2.0, 0.5*(level0+level1));
1614  //if (normalDist > 1*avgCellSize)
1615  //{
1616  // Pout<< "** skipping large distance " << endl;
1617  // return false;
1618  //}
1619 
1620  // Check average normal with respect to intersection locations
1621  if (mag(avgNormal&d/magD) > Foam::cos(degToRad(45.0)))
1622  {
1623  return true;
1624  }
1625  }
1626  }
1627 
1628  return false;
1629 }
1630 
1631 
1632 bool Foam::meshRefinement::checkProximity
1633 (
1634  const scalar planarCos,
1635  const label nAllowRefine,
1636 
1637  const label surfaceLevel, // current intersection max level
1638  const vector& surfaceLocation, // current intersection location
1639  const vector& surfaceNormal, // current intersection normal
1640 
1641  const label cellI,
1642 
1643  label& cellMaxLevel, // cached max surface level for this cell
1644  vector& cellMaxLocation, // cached surface location for this cell
1645  vector& cellMaxNormal, // cached surface normal for this cell
1646 
1648  label& nRefine
1649 ) const
1650 {
1651  const labelList& cellLevel = meshCutter_.cellLevel();
1652 
1653  // Test if surface applicable
1654  if (surfaceLevel > cellLevel[cellI])
1655  {
1656  if (cellMaxLevel == -1)
1657  {
1658  // First visit of cell. Store
1659  cellMaxLevel = surfaceLevel;
1660  cellMaxLocation = surfaceLocation;
1661  cellMaxNormal = surfaceNormal;
1662  }
1663  else
1664  {
1665  // Second or more visit.
1666  // Check if
1667  // - different location
1668  // - opposite surface
1669 
1670  bool closeSurfaces = isNormalGap
1671  (
1672  planarCos,
1673  cellLevel[cellI], //cellMaxLevel,
1674  cellMaxLocation,
1675  cellMaxNormal,
1676  cellLevel[cellI], //surfaceLevel,
1677  surfaceLocation,
1678  surfaceNormal
1679  );
1680 
1681  // Set normal to that of highest surface. Not really necessary
1682  // over here but we reuse cellMax info when doing coupled faces.
1683  if (surfaceLevel > cellMaxLevel)
1684  {
1685  cellMaxLevel = surfaceLevel;
1686  cellMaxLocation = surfaceLocation;
1687  cellMaxNormal = surfaceNormal;
1688  }
1689 
1690 
1691  if (closeSurfaces)
1692  {
1693  //Pout<< "Found gap:" << nl
1694  // << " location:" << surfaceLocation
1695  // << "\tnormal :" << surfaceNormal << nl
1697  // << "\tnormal :" << cellMaxNormal << nl
1698  // << "\tcos :" << (surfaceNormal&cellMaxNormal) << nl
1699  // << endl;
1700 
1701  return markForRefine
1702  (
1703  surfaceLevel, // mark with any non-neg number.
1704  nAllowRefine,
1705  refineCell[cellI],
1706  nRefine
1707  );
1708  }
1709  }
1710  }
1711 
1712  // Did not reach refinement limit.
1713  return true;
1714 }
1715 
1716 
1717 Foam::label Foam::meshRefinement::markProximityRefinement
1718 (
1719  const scalar planarCos,
1720 
1721  const labelList& surfaceMinLevel,
1722  const labelList& surfaceMaxLevel,
1723 
1724  const label nAllowRefine,
1725  const labelList& neiLevel,
1726  const pointField& neiCc,
1727 
1728  labelList& refineCell,
1729  label& nRefine
1730 ) const
1731 {
1732  const labelList& cellLevel = meshCutter_.cellLevel();
1733 
1734  label oldNRefine = nRefine;
1735 
1736  // 1. local test: any cell on more than one surface gets refined
1737  // (if its current level is < max of the surface max level)
1738 
1739  // 2. 'global' test: any cell on only one surface with a neighbour
1740  // on a different surface gets refined (if its current level etc.)
1741 
1742 
1743  // Collect candidate faces (i.e. intersecting any surface and
1744  // owner/neighbour not yet refined.
1745  const labelList testFaces(getRefineCandidateFaces(refineCell));
1746 
1747  // Collect segments
1748  pointField start(testFaces.size());
1749  pointField end(testFaces.size());
1750  labelList minLevel(testFaces.size());
1751 
1752  calcCellCellRays
1753  (
1754  neiCc,
1755  neiLevel,
1756  testFaces,
1757  start,
1758  end,
1759  minLevel
1760  );
1761 
1762 
1763  // Test for all intersections (with surfaces of higher gap level than
1764  // minLevel) and cache per cell the max surface level and the local normal
1765  // on that surface.
1766  labelList cellMaxLevel(mesh_.nCells(), -1);
1767  vectorField cellMaxNormal(mesh_.nCells(), Zero);
1768  pointField cellMaxLocation(mesh_.nCells(), Zero);
1769 
1770  {
1771  // Per segment the normals of the surfaces
1772  List<vectorList> surfaceLocation;
1773  List<vectorList> surfaceNormal;
1774  // Per segment the list of levels of the surfaces
1775  labelListList surfaceLevel;
1776 
1777  surfaces_.findAllIntersections
1778  (
1779  start,
1780  end,
1781  minLevel, // gap level of surface has to be bigger
1782  // than min level of neighbouring cells
1783 
1784  surfaceMinLevel,
1785  surfaceMaxLevel,
1786 
1787  surfaceLocation,
1788  surfaceNormal,
1789  surfaceLevel
1790  );
1791  // Clear out unnecessary data
1792  start.clear();
1793  end.clear();
1794  minLevel.clear();
1795 
1797  //OBJstream str
1798  //(
1799  // mesh_.time().path()
1800  // / "findAllIntersections_"
1801  // + mesh_.time().timeName()
1802  // + ".obj"
1803  //);
1805  //OBJstream str2
1806  //(
1807  // mesh_.time().path()
1808  // / "findAllIntersections2_"
1809  // + mesh_.time().timeName()
1810  // + ".obj"
1811  //);
1812 
1813  forAll(testFaces, i)
1814  {
1815  label faceI = testFaces[i];
1816  label own = mesh_.faceOwner()[faceI];
1817 
1818  const labelList& fLevels = surfaceLevel[i];
1819  const vectorList& fPoints = surfaceLocation[i];
1820  const vectorList& fNormals = surfaceNormal[i];
1821 
1822  forAll(fLevels, hitI)
1823  {
1824  checkProximity
1825  (
1826  planarCos,
1827  nAllowRefine,
1828 
1829  fLevels[hitI],
1830  fPoints[hitI],
1831  fNormals[hitI],
1832 
1833  own,
1834  cellMaxLevel[own],
1835  cellMaxLocation[own],
1836  cellMaxNormal[own],
1837 
1838  refineCell,
1839  nRefine
1840  );
1841  }
1842 
1843  if (mesh_.isInternalFace(faceI))
1844  {
1845  label nei = mesh_.faceNeighbour()[faceI];
1846 
1847  forAll(fLevels, hitI)
1848  {
1849  checkProximity
1850  (
1851  planarCos,
1852  nAllowRefine,
1853 
1854  fLevels[hitI],
1855  fPoints[hitI],
1856  fNormals[hitI],
1857 
1858  nei,
1859  cellMaxLevel[nei],
1860  cellMaxLocation[nei],
1861  cellMaxNormal[nei],
1862 
1863  refineCell,
1864  nRefine
1865  );
1866  }
1867  }
1868  }
1869  }
1870 
1871  // 2. Find out a measure of surface curvature and region edges.
1872  // Send over surface region and surface normal to neighbour cell.
1873 
1874  labelList neiBndMaxLevel(mesh_.nBoundaryFaces());
1875  pointField neiBndMaxLocation(mesh_.nBoundaryFaces());
1876  vectorField neiBndMaxNormal(mesh_.nBoundaryFaces());
1877 
1878  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
1879  {
1880  label bFaceI = faceI-mesh_.nInternalFaces();
1881  label own = mesh_.faceOwner()[faceI];
1882 
1883  neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
1884  neiBndMaxLocation[bFaceI] = cellMaxLocation[own];
1885  neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
1886  }
1887  syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
1888  syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLocation);
1889  syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
1890 
1891  // Loop over all faces. Could only be checkFaces.. except if they're coupled
1892 
1893  // Internal faces
1894  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
1895  {
1896  label own = mesh_.faceOwner()[faceI];
1897  label nei = mesh_.faceNeighbour()[faceI];
1898 
1899  if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
1900  {
1901  // Have valid data on both sides. Check planarCos.
1902  if
1903  (
1904  isNormalGap
1905  (
1906  planarCos,
1907  cellLevel[own], //cellMaxLevel[own],
1908  cellMaxLocation[own],
1909  cellMaxNormal[own],
1910  cellLevel[nei], //cellMaxLevel[nei],
1911  cellMaxLocation[nei],
1912  cellMaxNormal[nei]
1913  )
1914  )
1915  {
1916  // See which side to refine
1917  if (cellLevel[own] < cellMaxLevel[own])
1918  {
1919  if
1920  (
1921  !markForRefine
1922  (
1923  cellMaxLevel[own],
1924  nAllowRefine,
1925  refineCell[own],
1926  nRefine
1927  )
1928  )
1929  {
1930  if (debug)
1931  {
1932  Pout<< "Stopped refining since reaching my cell"
1933  << " limit of " << mesh_.nCells()+7*nRefine
1934  << endl;
1935  }
1936  break;
1937  }
1938  }
1939 
1940  if (cellLevel[nei] < cellMaxLevel[nei])
1941  {
1942  if
1943  (
1944  !markForRefine
1945  (
1946  cellMaxLevel[nei],
1947  nAllowRefine,
1948  refineCell[nei],
1949  nRefine
1950  )
1951  )
1952  {
1953  if (debug)
1954  {
1955  Pout<< "Stopped refining since reaching my cell"
1956  << " limit of " << mesh_.nCells()+7*nRefine
1957  << endl;
1958  }
1959  break;
1960  }
1961  }
1962  }
1963  }
1964  }
1965  // Boundary faces
1966  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
1967  {
1968  label own = mesh_.faceOwner()[faceI];
1969  label bFaceI = faceI - mesh_.nInternalFaces();
1970 
1971  if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
1972  {
1973  // Have valid data on both sides. Check planarCos.
1974  if
1975  (
1976  isNormalGap
1977  (
1978  planarCos,
1979  cellLevel[own], //cellMaxLevel[own],
1980  cellMaxLocation[own],
1981  cellMaxNormal[own],
1982  neiLevel[bFaceI], //neiBndMaxLevel[bFaceI],
1983  neiBndMaxLocation[bFaceI],
1984  neiBndMaxNormal[bFaceI]
1985  )
1986  )
1987  {
1988  if
1989  (
1990  !markForRefine
1991  (
1992  cellMaxLevel[own],
1993  nAllowRefine,
1994  refineCell[own],
1995  nRefine
1996  )
1997  )
1998  {
1999  if (debug)
2000  {
2001  Pout<< "Stopped refining since reaching my cell"
2002  << " limit of " << mesh_.nCells()+7*nRefine
2003  << endl;
2004  }
2005  break;
2006  }
2007  }
2008  }
2009  }
2010 
2011  if
2012  (
2013  returnReduce(nRefine, sumOp<label>())
2014  > returnReduce(nAllowRefine, sumOp<label>())
2015  )
2016  {
2017  Info<< "Reached refinement limit." << endl;
2018  }
2019 
2020  return returnReduce(nRefine-oldNRefine, sumOp<label>());
2021 }
2022 
2023 
2024 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2025 
2026 // Calculate list of cells to refine. Gets for any edge (start - end)
2027 // whether it intersects the surface. Does more accurate test and checks
2028 // the wanted level on the surface intersected.
2029 // Does approximate precalculation of how many cells can be refined before
2030 // hitting overall limit maxGlobalCells.
2033  const pointField& keepPoints,
2034  const scalar curvature,
2035  const scalar planarAngle,
2036 
2037  const bool featureRefinement,
2038  const bool featureDistanceRefinement,
2039  const bool internalRefinement,
2040  const bool surfaceRefinement,
2041  const bool curvatureRefinement,
2042  const bool smallFeatureRefinement,
2043  const bool gapRefinement,
2044  const bool bigGapRefinement,
2045  const bool spreadGapSize,
2046  const label maxGlobalCells,
2047  const label maxLocalCells
2048 ) const
2049 {
2050  label totNCells = mesh_.globalData().nTotalCells();
2051 
2052  labelList cellsToRefine;
2053 
2054  if (totNCells >= maxGlobalCells)
2055  {
2056  Info<< "No cells marked for refinement since reached limit "
2057  << maxGlobalCells << '.' << endl;
2058  }
2059  else
2060  {
2061  // Every cell I refine adds 7 cells. Estimate the number of cells
2062  // I am allowed to refine.
2063  // Assume perfect distribution so can only refine as much the fraction
2064  // of the mesh I hold. This prediction step prevents us having to do
2065  // lots of reduces to keep count of the total number of cells selected
2066  // for refinement.
2067 
2068  //scalar fraction = scalar(mesh_.nCells())/totNCells;
2069  //label myMaxCells = label(maxGlobalCells*fraction);
2070  //label nAllowRefine = (myMaxCells - mesh_.nCells())/7;
2072  //
2073  //Pout<< "refineCandidates:" << nl
2074  // << " total cells:" << totNCells << nl
2075  // << " local cells:" << mesh_.nCells() << nl
2076  // << " local fraction:" << fraction << nl
2077  // << " local allowable cells:" << myMaxCells << nl
2078  // << " local allowable refinement:" << nAllowRefine << nl
2079  // << endl;
2080 
2081  //- Disable refinement shortcut. nAllowRefine is per processor limit.
2082  label nAllowRefine = labelMax / Pstream::nProcs();
2083 
2084  // Marked for refinement (>= 0) or not (-1). Actual value is the
2085  // index of the surface it intersects / shell it is inside.
2086  labelList refineCell(mesh_.nCells(), -1);
2087  label nRefine = 0;
2088 
2089 
2090  // Swap neighbouring cell centres and cell level
2091  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2092  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2093  calcNeighbourData(neiLevel, neiCc);
2094 
2095 
2096  const scalar planarCos = Foam::cos(degToRad(planarAngle));
2097 
2098 
2099  // Cells pierced by feature lines
2100  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2101 
2102  if (featureRefinement)
2103  {
2104  label nFeatures = markFeatureRefinement
2105  (
2106  keepPoints,
2107  nAllowRefine,
2108 
2109  refineCell,
2110  nRefine
2111  );
2112 
2113  Info<< "Marked for refinement due to explicit features "
2114  << ": " << nFeatures << " cells." << endl;
2115  }
2116 
2117  // Inside distance-to-feature shells
2118  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2119 
2120  if (featureDistanceRefinement)
2121  {
2122  label nShell = markInternalDistanceToFeatureRefinement
2123  (
2124  nAllowRefine,
2125 
2126  refineCell,
2127  nRefine
2128  );
2129  Info<< "Marked for refinement due to distance to explicit features "
2130  ": " << nShell << " cells." << endl;
2131  }
2132 
2133  // Inside refinement shells
2134  // ~~~~~~~~~~~~~~~~~~~~~~~~
2135 
2136  if (internalRefinement)
2137  {
2138  label nShell = markInternalRefinement
2139  (
2140  nAllowRefine,
2141 
2142  refineCell,
2143  nRefine
2144  );
2145  Info<< "Marked for refinement due to refinement shells "
2146  << ": " << nShell << " cells." << endl;
2147  }
2148 
2149  // Refinement based on intersection of surface
2150  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2151 
2152  if (surfaceRefinement)
2153  {
2154  label nSurf = markSurfaceRefinement
2155  (
2156  nAllowRefine,
2157  neiLevel,
2158  neiCc,
2159 
2160  refineCell,
2161  nRefine
2162  );
2163  Info<< "Marked for refinement due to surface intersection "
2164  << ": " << nSurf << " cells." << endl;
2165 
2166  // Refine intersected-cells only inside gaps. See
2167  // markInternalGapRefinement to refine all cells inside gaps.
2168  if
2169  (
2170  planarCos >= -1
2171  && planarCos <= 1
2172  && max(shells_.maxGapLevel()) > 0
2173  )
2174  {
2175  label nGapSurf = markSurfaceGapRefinement
2176  (
2177  planarCos,
2178  nAllowRefine,
2179  neiLevel,
2180  neiCc,
2181 
2182  refineCell,
2183  nRefine
2184  );
2185  Info<< "Marked for refinement due to surface intersection"
2186  << " (at gaps)"
2187  << ": " << nGapSurf << " cells." << endl;
2188  }
2189  }
2190 
2191  // Refinement based on curvature of surface
2192  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2193 
2194  if
2195  (
2196  curvatureRefinement
2197  && (curvature >= -1 && curvature <= 1)
2198  && (surfaces_.minLevel() != surfaces_.maxLevel())
2199  )
2200  {
2201  label nCurv = markSurfaceCurvatureRefinement
2202  (
2203  curvature,
2204  nAllowRefine,
2205  neiLevel,
2206  neiCc,
2207 
2208  refineCell,
2209  nRefine
2210  );
2211  Info<< "Marked for refinement due to curvature/regions "
2212  << ": " << nCurv << " cells." << endl;
2213  }
2214 
2215 
2216  // Refinement based on features smaller than cell size
2217  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2218 
2219  if
2220  (
2221  smallFeatureRefinement
2222  && (planarCos >= -1 && planarCos <= 1)
2223  && max(shells_.maxGapLevel()) > 0
2224  )
2225  {
2226  label nGap = markSmallFeatureRefinement
2227  (
2228  planarCos,
2229  nAllowRefine,
2230  neiLevel,
2231  neiCc,
2232 
2233  refineCell,
2234  nRefine
2235  );
2236  Info<< "Marked for refinement due to close opposite surfaces "
2237  << ": " << nGap << " cells." << endl;
2238  }
2239 
2240 
2241  // Refinement based on gap (only neighbouring cells)
2242  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2243 
2244  const labelList& surfaceGapLevel = surfaces_.gapLevel();
2245 
2246  if
2247  (
2248  gapRefinement
2249  && (planarCos >= -1 && planarCos <= 1)
2250  && (max(surfaceGapLevel) > -1)
2251  )
2252  {
2253  Info<< "Specified gap level : " << max(surfaceGapLevel)
2254  << ", planar angle " << planarAngle << endl;
2255 
2256  label nGap = markProximityRefinement
2257  (
2258  planarCos,
2259 
2260  labelList(surfaceGapLevel.size(), -1), // surface min level
2261  surfaceGapLevel, // surface max level
2262 
2263  nAllowRefine,
2264  neiLevel,
2265  neiCc,
2266 
2267  refineCell,
2268  nRefine
2269  );
2270  Info<< "Marked for refinement due to close opposite surfaces "
2271  << ": " << nGap << " cells." << endl;
2272  }
2273 
2274 
2275  // Refinement based on gaps larger than cell size
2276  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2277 
2278  if
2279  (
2280  bigGapRefinement
2281  && (planarCos >= -1 && planarCos <= 1)
2282  && max(shells_.maxGapLevel()) > 0
2283  )
2284  {
2285  // Refine based on gap information provided by shell and nearest
2286  // surface
2287  labelList numGapCells;
2288  scalarField gapSize;
2289  label nGap = markInternalGapRefinement
2290  (
2291  planarCos,
2292  spreadGapSize,
2293  nAllowRefine,
2294 
2295  refineCell,
2296  nRefine,
2297  numGapCells,
2298  gapSize
2299  );
2300  Info<< "Marked for refinement due to opposite surfaces"
2301  << " "
2302  << ": " << nGap << " cells." << endl;
2303  }
2304 
2305 
2306  // Limit refinement
2307  // ~~~~~~~~~~~~~~~~
2308 
2309  {
2310  label nUnmarked = unmarkInternalRefinement(refineCell, nRefine);
2311  if (nUnmarked > 0)
2312  {
2313  Info<< "Unmarked for refinement due to limit shells"
2314  << " : " << nUnmarked << " cells." << endl;
2315  }
2316  }
2317 
2318 
2319 
2320  // Pack cells-to-refine
2321  // ~~~~~~~~~~~~~~~~~~~~
2322 
2323  cellsToRefine.setSize(nRefine);
2324  nRefine = 0;
2325 
2326  forAll(refineCell, cellI)
2327  {
2328  if (refineCell[cellI] != -1)
2329  {
2330  cellsToRefine[nRefine++] = cellI;
2331  }
2332  }
2333  }
2334 
2335  return cellsToRefine;
2336 }
2337 
2338 
2341  const labelList& cellsToRefine
2342 )
2343 {
2344  // Mesh changing engine.
2345  polyTopoChange meshMod(mesh_);
2346 
2347  // Play refinement commands into mesh changer.
2348  meshCutter_.setRefinement(cellsToRefine, meshMod);
2349 
2350  // Create mesh (no inflation), return map from old to new mesh.
2351  autoPtr<mapPolyMesh> mapPtr = meshMod.changeMesh(mesh_, false);
2352  mapPolyMesh& map = *mapPtr;
2353 
2354  // Update fields
2355  mesh_.updateMesh(map);
2356 
2357  // Optionally inflate mesh
2358  if (map.hasMotionPoints())
2359  {
2360  mesh_.movePoints(map.preMotionPoints());
2361  }
2362  else
2363  {
2364  // Delete mesh volumes.
2365  mesh_.clearOut();
2366  }
2367 
2368  // Reset the instance for if in overwrite mode
2369  mesh_.setInstance(timeName());
2370 
2371  // Update intersection info
2372  updateMesh(map, getChangedFaces(map, cellsToRefine));
2373 
2374  return mapPtr;
2375 }
2376 
2377 
2378 // Do refinement of consistent set of cells followed by truncation and
2379 // load balancing.
2383  const string& msg,
2384  decompositionMethod& decomposer,
2385  fvMeshDistribute& distributor,
2386  const labelList& cellsToRefine,
2387  const scalar maxLoadUnbalance
2388 )
2389 {
2390  // Do all refinement
2391  refine(cellsToRefine);
2392 
2394  {
2395  Pout<< "Writing refined but unbalanced " << msg
2396  << " mesh to time " << timeName() << endl;
2397  write
2398  (
2399  debugType(debug),
2400  writeType(writeLevel() | WRITEMESH),
2401  mesh_.time().path()/timeName()
2402  );
2403  Pout<< "Dumped debug data in = "
2404  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2405 
2406  // test all is still synced across proc patches
2407  checkData();
2408  }
2409 
2410  Info<< "Refined mesh in = "
2411  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2412  printMeshInfo(debug, "After refinement " + msg);
2413 
2414 
2415  // Load balancing
2416  // ~~~~~~~~~~~~~~
2417 
2419 
2420  if (Pstream::nProcs() > 1)
2421  {
2422  scalar nIdealCells =
2423  mesh_.globalData().nTotalCells()
2424  / Pstream::nProcs();
2425 
2426  scalar unbalance = returnReduce
2427  (
2428  mag(1.0-mesh_.nCells()/nIdealCells),
2429  maxOp<scalar>()
2430  );
2431 
2432  if (unbalance <= maxLoadUnbalance)
2433  {
2434  Info<< "Skipping balancing since max unbalance " << unbalance
2435  << " is less than allowable " << maxLoadUnbalance
2436  << endl;
2437  }
2438  else
2439  {
2440  scalarField cellWeights(mesh_.nCells(), 1);
2441 
2442  distMap = balance
2443  (
2444  false, //keepZoneFaces
2445  false, //keepBaffles
2446  cellWeights,
2447  decomposer,
2448  distributor
2449  );
2450 
2451  Info<< "Balanced mesh in = "
2452  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2453 
2454  printMeshInfo(debug, "After balancing " + msg);
2455 
2456 
2458  {
2459  Pout<< "Writing balanced " << msg
2460  << " mesh to time " << timeName() << endl;
2461  write
2462  (
2463  debugType(debug),
2464  writeType(writeLevel() | WRITEMESH),
2465  mesh_.time().path()/timeName()
2466  );
2467  Pout<< "Dumped debug data in = "
2468  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2469 
2470  // test all is still synced across proc patches
2471  checkData();
2472  }
2473  }
2474  }
2475 
2476  return distMap;
2477 }
2478 
2479 
2480 // Do load balancing followed by refinement of consistent set of cells.
2484  const string& msg,
2485  decompositionMethod& decomposer,
2486  fvMeshDistribute& distributor,
2487  const labelList& initCellsToRefine,
2488  const scalar maxLoadUnbalance
2489 )
2490 {
2491  labelList cellsToRefine(initCellsToRefine);
2492 
2493  //{
2494  // globalIndex globalCells(mesh_.nCells());
2495  //
2496  // Info<< "** Distribution before balancing/refining:" << endl;
2497  // for (const int procI : Pstream::allProcs())
2498  // {
2499  // Info<< " " << procI << '\t'
2500  // << globalCells.localSize(procI) << endl;
2501  // }
2502  // Info<< endl;
2503  //}
2504  //{
2505  // globalIndex globalCells(cellsToRefine.size());
2506  //
2507  // Info<< "** Cells to be refined:" << endl;
2508  // for (const int procI : Pstream::allProcs())
2509  // {
2510  // Info<< " " << procI << '\t'
2511  // << globalCells.localSize(procI) << endl;
2512  // }
2513  // Info<< endl;
2514  //}
2515 
2516 
2517  // Load balancing
2518  // ~~~~~~~~~~~~~~
2519 
2521 
2522  if (Pstream::nProcs() > 1)
2523  {
2524  // First check if we need to balance at all. Precalculate number of
2525  // cells after refinement and see what maximum difference is.
2526  scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
2527  scalar nIdealNewCells =
2528  returnReduce(nNewCells, sumOp<scalar>())
2529  / Pstream::nProcs();
2530  scalar unbalance = returnReduce
2531  (
2532  mag(1.0-nNewCells/nIdealNewCells),
2533  maxOp<scalar>()
2534  );
2535 
2536  if (unbalance <= maxLoadUnbalance)
2537  {
2538  Info<< "Skipping balancing since max unbalance " << unbalance
2539  << " is less than allowable " << maxLoadUnbalance
2540  << endl;
2541  }
2542  else
2543  {
2544  scalarField cellWeights(mesh_.nCells(), 1);
2545  forAll(cellsToRefine, i)
2546  {
2547  cellWeights[cellsToRefine[i]] += 7;
2548  }
2549 
2550  distMap = balance
2551  (
2552  false, //keepZoneFaces
2553  false, //keepBaffles
2554  cellWeights,
2555  decomposer,
2556  distributor
2557  );
2558 
2559  // Update cells to refine
2560  distMap().distributeCellIndices(cellsToRefine);
2561 
2562  Info<< "Balanced mesh in = "
2563  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2564  }
2565 
2566  //{
2567  // globalIndex globalCells(mesh_.nCells());
2568  //
2569  // Info<< "** Distribution after balancing:" << endl;
2570  // for (const int procI : Pstream::allProcs())
2571  // {
2572  // Info<< " " << procI << '\t'
2573  // << globalCells.localSize(procI) << endl;
2574  // }
2575  // Info<< endl;
2576  //}
2577 
2578  printMeshInfo(debug, "After balancing " + msg);
2579 
2581  {
2582  Pout<< "Writing balanced " << msg
2583  << " mesh to time " << timeName() << endl;
2584  write
2585  (
2586  debugType(debug),
2587  writeType(writeLevel() | WRITEMESH),
2588  mesh_.time().path()/timeName()
2589  );
2590  Pout<< "Dumped debug data in = "
2591  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2592 
2593  // test all is still synced across proc patches
2594  checkData();
2595  }
2596  }
2597 
2598 
2599  // Refinement
2600  // ~~~~~~~~~~
2601 
2602  refine(cellsToRefine);
2603 
2605  {
2606  Pout<< "Writing refined " << msg
2607  << " mesh to time " << timeName() << endl;
2608  write
2609  (
2610  debugType(debug),
2611  writeType(writeLevel() | WRITEMESH),
2612  mesh_.time().path()/timeName()
2613  );
2614  Pout<< "Dumped debug data in = "
2615  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2616 
2617  // test all is still synced across proc patches
2618  checkData();
2619  }
2620 
2621  Info<< "Refined mesh in = "
2622  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2623 
2624  //{
2625  // globalIndex globalCells(mesh_.nCells());
2626  //
2627  // Info<< "** After refinement distribution:" << endl;
2628  // for (const int procI : Pstream::allProcs())
2629  // {
2630  // Info<< " " << procI << '\t'
2631  // << globalCells.localSize(procI) << endl;
2632  // }
2633  // Info<< endl;
2634  //}
2635 
2636  printMeshInfo(debug, "After refinement " + msg);
2637 
2638  return distMap;
2639 }
2640 
2641 
2644  const label maxGlobalCells,
2645  const label maxLocalCells,
2646  const labelList& currentLevel,
2647  const direction dir
2648 ) const
2649 {
2650  const labelList& cellLevel = meshCutter_.cellLevel();
2651  const pointField& cellCentres = mesh_.cellCentres();
2652 
2653  label totNCells = mesh_.globalData().nTotalCells();
2654 
2655  labelList cellsToRefine;
2656 
2657  if (totNCells >= maxGlobalCells)
2658  {
2659  Info<< "No cells marked for refinement since reached limit "
2660  << maxGlobalCells << '.' << endl;
2661  }
2662  else
2663  {
2664  // Disable refinement shortcut. nAllowRefine is per processor limit.
2665  label nAllowRefine = labelMax / Pstream::nProcs();
2666 
2667  // Marked for refinement (>= 0) or not (-1). Actual value is the
2668  // index of the surface it intersects / shell it is inside
2669  labelList refineCell(mesh_.nCells(), -1);
2670  label nRefine = 0;
2671 
2672  // Find cells inside the shells with directional levels
2673  labelList insideShell;
2674  shells_.findDirectionalLevel
2675  (
2676  cellCentres,
2677  cellLevel,
2678  currentLevel, // current directional level
2679  dir,
2680  insideShell
2681  );
2682 
2683  // Mark for refinement
2684  forAll(insideShell, celli)
2685  {
2686  if (insideShell[celli] >= 0)
2687  {
2688  bool reachedLimit = !markForRefine
2689  (
2690  insideShell[celli], // mark with any positive value
2691  nAllowRefine,
2692  refineCell[celli],
2693  nRefine
2694  );
2695 
2696  if (reachedLimit)
2697  {
2698  if (debug)
2699  {
2700  Pout<< "Stopped refining cells"
2701  << " since reaching my cell limit of "
2702  << mesh_.nCells()+nAllowRefine << endl;
2703  }
2704  break;
2705  }
2706  }
2707  }
2708 
2709  // Limit refinement
2710  // ~~~~~~~~~~~~~~~~
2711 
2712  {
2713  label nUnmarked = unmarkInternalRefinement(refineCell, nRefine);
2714  if (nUnmarked > 0)
2715  {
2716  Info<< "Unmarked for refinement due to limit shells"
2717  << " : " << nUnmarked << " cells." << endl;
2718  }
2719  }
2720 
2721 
2722 
2723  // Pack cells-to-refine
2724  // ~~~~~~~~~~~~~~~~~~~~
2725 
2726  cellsToRefine.setSize(nRefine);
2727  nRefine = 0;
2728 
2729  forAll(refineCell, cellI)
2730  {
2731  if (refineCell[cellI] != -1)
2732  {
2733  cellsToRefine[nRefine++] = cellI;
2734  }
2735  }
2736  }
2737 
2738  return cellsToRefine;
2739 }
2740 
2741 
2744  const string& msg,
2745  const direction cmpt,
2746  const labelList& cellsToRefine
2747 )
2748 {
2749  // Set splitting direction
2750  vector refDir(Zero);
2751  refDir[cmpt] = 1;
2752  List<refineCell> refCells(cellsToRefine.size());
2753  forAll(cellsToRefine, i)
2754  {
2755  refCells[i] = refineCell(cellsToRefine[i], refDir);
2756  }
2757 
2758  // How to walk circumference of cells
2759  hexCellLooper cellWalker(mesh_);
2760 
2761  // Analyse cuts
2762  cellCuts cuts(mesh_, cellWalker, refCells);
2763 
2764  // Cell cutter
2765  Foam::meshCutter meshRefiner(mesh_);
2766 
2767  polyTopoChange meshMod(mesh_);
2768 
2769  // Insert mesh refinement into polyTopoChange.
2770  meshRefiner.setRefinement(cuts, meshMod);
2771 
2772  autoPtr<mapPolyMesh> morphMap = meshMod.changeMesh(mesh_, false);
2773 
2774  // Update fields
2775  mesh_.updateMesh(*morphMap);
2776 
2777  // Move mesh (since morphing does not do this)
2778  if (morphMap().hasMotionPoints())
2779  {
2780  mesh_.movePoints(morphMap().preMotionPoints());
2781  }
2782  else
2783  {
2784  // Delete mesh volumes.
2785  mesh_.clearOut();
2786  }
2787 
2788  // Reset the instance for if in overwrite mode
2789  mesh_.setInstance(timeName());
2790 
2791  // Update stored refinement pattern
2792  meshRefiner.updateMesh(*morphMap);
2793 
2794  // Update intersection info
2795  updateMesh(*morphMap, getChangedFaces(*morphMap, cellsToRefine));
2796 
2797  return morphMap;
2798 }
2799 
2800 
2801 // ************************************************************************* //
Foam::expressions::patchExpr::debug
int debug
Static debugging option.
shellSurfaces.H
Foam::labelList
List< label > labelList
A List of labels.
Definition: List.H:67
Foam::pointField
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:44
Foam::pTraits< vectorList >::cmptType
vectorList cmptType
Component type.
Definition: meshRefinementRefine.C:113
Foam::maxOp
Definition: ops.H:223
meshCutter.H
Foam::meshCutter
Cuts (splits) cells.
Definition: meshCutter.H:138
Foam::BitOps::set
void set(List< bool > &bools, const labelRange &range)
Set the specified range 'on' in a boolList.
Definition: BitOps.C:37
Foam::meshRefinement::isNormalGap
bool isNormalGap(const scalar planarCos, const label level0, const vector &point0, const vector &normal0, const label level1, const vector &point1, const vector &normal1) const
Is local topology a small gap normal to the test vector.
Definition: meshRefinementRefine.C:1575
Foam::labelMax
constexpr label labelMax
Definition: label.H:61
Foam::returnReduce
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Definition: PstreamReduceOps.H:94
Foam::Zero
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
Foam::normalLess::normalLess
normalLess(const vectorList &values)
Definition: meshRefinementRefine.C:82
Foam::HashTableOps::values
List< T > values(const HashTable< T, Key, Hash > &tbl, const bool doSort=false)
List of values from HashTable, optionally sorted.
Definition: HashOps.H:149
Foam::less
static bool less(const vector &x, const vector &y)
To compare normals.
Definition: meshRefinementRefine.C:57
Foam::checkData
label checkData(const fvMesh &mesh, const instantList &timeDirs, wordList &objectNames)
Check if fields are good to use (available at all times)
Foam::syncTools::syncBoundaryFaceList
static void syncBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop, const TransformOp &top, const bool parRun=Pstream::parRun())
Synchronize values on boundary faces only.
Definition: syncToolsTemplates.C:998
polyTopoChange.H
Foam::primitiveMesh::cells
const cellList & cells() const
Definition: primitiveMeshCells.C:138
Foam::surfaceLocation
Contains information about location on a triSurface.
Definition: surfaceLocation.H:71
Foam::polyTopoChange
Direct mesh changes based on v1.3 polyTopoChange syntax.
Definition: polyTopoChange.H:99
Cloud.H
Foam::syncTools::swapBoundaryFaceList
static void swapBoundaryFaceList(const polyMesh &mesh, UList< T > &faceValues)
Swap coupled boundary face values. Uses eqOp.
Definition: syncTools.H:445
Foam::boolList
List< bool > boolList
A List of bools.
Definition: List.H:65
Foam::polyMesh::boundaryMesh
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:444
Foam::pTraits< labelList >::cmptType
labelList cmptType
Component type.
Definition: meshRefinementRefine.C:102
Foam::endl
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:369
decompositionMethod.H
Foam::mapPolyMesh::preMotionPoints
const pointField & preMotionPoints() const
Pre-motion point positions.
Definition: mapPolyMesh.H:613
Foam::syncTools::getMasterFaces
static bitSet getMasterFaces(const polyMesh &mesh)
Definition: syncTools.C:126
Foam::meshRefinement::balanceAndRefine
autoPtr< mapDistributePolyMesh > balanceAndRefine(const string &msg, decompositionMethod &decomposer, fvMeshDistribute &distributor, const labelList &cellsToRefine, const scalar maxLoadUnbalance)
Balance before refining some cells.
Definition: meshRefinementRefine.C:2483
refinementFeatures.H
Foam::Pout
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
Foam::meshRefinement::refineAndBalance
autoPtr< mapDistributePolyMesh > refineAndBalance(const string &msg, decompositionMethod &decomposer, fvMeshDistribute &distributor, const labelList &cellsToRefine, const scalar maxLoadUnbalance)
Refine some cells and rebalance.
Definition: meshRefinementRefine.C:2382
Foam::meshCutter::updateMesh
void updateMesh(const mapPolyMesh &)
Force recalculation of locally stored data on topological change.
Definition: meshCutter.C:995
Foam::meshRefinement::mesh
const fvMesh & mesh() const
Reference to mesh.
Definition: meshRefinement.H:944
syncTools.H
Foam::polyBoundaryMesh::start
label start() const
The start label of the boundary faces in the polyMesh face list.
Definition: polyBoundaryMesh.C:611
Foam::min
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Foam::sumOp
Definition: ops.H:213
forAll
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:296
Foam::vectorField
Field< vector > vectorField
Specialisation of Field<T> for vector.
Definition: primitiveFieldsFwd.H:54
Foam::magSqr
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
Foam::meshRefinement::refineCandidates
labelList refineCandidates(const pointField &keepPoints, const scalar curvature, const scalar planarAngle, const bool featureRefinement, const bool featureDistanceRefinement, const bool internalRefinement, const bool surfaceRefinement, const bool curvatureRefinement, const bool smallFeatureRefinement, const bool gapRefinement, const bool bigGapRefinement, const bool spreadGapSize, const label maxGlobalCells, const label maxLocalCells) const
Calculate list of cells to refine.
Definition: meshRefinementRefine.C:2032
mapDistributePolyMesh.H
Foam::reduce
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
Definition: PstreamReduceOps.H:51
Foam::syncTools::syncFaceList
static void syncFaceList(const polyMesh &mesh, UList< T > &faceValues, const CombineOp &cop)
Synchronize values on all mesh faces.
Definition: syncTools.H:396
Foam::meshRefinement::FEATURESEEDS
Definition: meshRefinement.H:96
Foam::constant::physicoChemical::b
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Foam::Field< vector >
Foam::meshRefinement::directionalRefine
autoPtr< mapPolyMesh > directionalRefine(const string &msg, const direction cmpt, const labelList &cellsToRefine)
Directionally refine in direction cmpt.
Definition: meshRefinementRefine.C:2743
Foam::Info
messageStream Info
Information stream (stdout output on master, null elsewhere)
refinementSurfaces.H
Foam::OSstream::name
virtual const fileName & name() const
Get the name of the stream.
Definition: OSstream.H:107
Foam::List::setSize
void setSize(const label n)
Alias for resize()
Definition: List.H:222
Foam::polyMesh::faceOwner
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1107
faceSet.H
Foam::cellList
List< cell > cellList
A List of cells.
Definition: cellListFwd.H:47
Foam::meshRefinement::directionalRefineCandidates
labelList directionalRefineCandidates(const label maxGlobalCells, const label maxLocalCells, const labelList &currentLevel, const direction dir) const
Calculate list of cells to directionally refine.
Definition: meshRefinementRefine.C:2643
Foam::primitiveMesh::nBoundaryFaces
label nBoundaryFaces() const noexcept
Number of boundary faces (== nFaces - nInternalFaces)
Definition: primitiveMeshI.H:84
refineCell.H
Foam::max
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
Foam::polyTopoChange::changeMesh
autoPtr< mapPolyMesh > changeMesh(polyMesh &mesh, const labelUList &patchMap, const bool inflate, const bool syncParallel=true, const bool orderCells=false, const bool orderPoints=false)
Inplace changes mesh without change of patches.
Definition: polyTopoChange.C:2980
Foam::syncTools::swapBoundaryCellList
static void swapBoundaryCellList(const polyMesh &mesh, const UList< T > &cellData, List< T > &neighbourCellData)
Swap to obtain neighbour cell values for all boundary faces.
Definition: syncToolsTemplates.C:1392
timeName
word timeName
Definition: getTimeIndex.H:3
cellCuts.H
Foam::normalLess::operator()
bool operator()(const label a, const label b) const
Definition: meshRefinementRefine.C:87
meshRefinement.H
Foam::decompositionMethod
Abstract base class for domain decomposition.
Definition: decompositionMethod.H:51
stdFoam::end
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:121
Foam::Ostream::write
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
Foam
Namespace for OpenFOAM.
Definition: atmBoundaryLayer.C:33
treeDataCell.H
trackedParticle.H
Foam::vector
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:51
Foam::degToRad
constexpr scalar degToRad(const scalar deg) noexcept
Conversion from degrees to radians.
Definition: unitConversion.H:48
Foam::meshRefinement::isGap
bool isGap(const scalar, const vector &, const vector &, const vector &, const vector &) const
Is local topology a small gap?
Definition: meshRefinementRefine.C:1528
Time.H
Foam::autoPtr
Pointer management similar to std::unique_ptr, with some additional methods and type checking.
Definition: HashPtrTable.H:53
Foam::labelListList
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:56
Foam::vectorList
List< vector > vectorList
A List of vectors.
Definition: vectorList.H:64
Foam::refineCell
Container with cells to refine. Refinement given as single direction.
Definition: refineCell.H:56
fvMeshDistribute.H
Foam::Vector< scalar >
Foam::meshRefinement::writeType
writeType
Enumeration for what to write. Used as a bit-pattern.
Definition: meshRefinement.H:112
Foam::normalLess
To compare normals.
Definition: meshRefinementRefine.C:76
Foam::meshCutter::setRefinement
void setRefinement(const cellCuts &cuts, polyTopoChange &meshMod)
Do actual cutting with cut description. Inserts mesh changes.
Definition: meshCutter.C:522
Foam::primitiveMesh::nInternalFaces
label nInternalFaces() const noexcept
Number of internal faces.
Definition: primitiveMeshI.H:78
Foam::List< vector >
Foam::pTraits
A traits class, which is primarily used for primitives.
Definition: pTraits.H:56
Foam::mag
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
Foam::constant::electromagnetic::e
const dimensionedScalar e
Elementary charge.
Definition: createFields.H:11
Foam::List::clear
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:116
Foam::direction
uint8_t direction
Definition: direction.H:52
x
x
Definition: LISASMDCalcMethod2.H:52
Foam::vtk::write
void write(vtk::formatter &fmt, const Type &val, const label n=1)
Component-wise write of a value (N times)
Definition: foamVtkOutputTemplates.C:36
Foam::mapPolyMesh
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:161
Foam::mapPolyMesh::hasMotionPoints
bool hasMotionPoints() const
Has valid preMotionPoints?
Definition: mapPolyMesh.H:619
Foam::VectorSpace< Vector< scalar >, scalar, 3 >::zero
static const Vector< scalar > zero
Definition: VectorSpace.H:115
Foam::findIndices
labelList findIndices(const ListType &input, typename ListType::const_reference val, label start=0)
Linear search to find all occurrences of given element.
Foam::primitiveMesh::nFaces
label nFaces() const noexcept
Number of mesh faces.
Definition: primitiveMeshI.H:90
Foam::meshRefinement::refine
autoPtr< mapPolyMesh > refine(const labelList &cellsToRefine)
Refine some cells.
Definition: meshRefinementRefine.C:2340
cells
const cellShapeList & cells
Definition: gmvOutputHeader.H:3
Foam::sortedOrder
labelList sortedOrder(const UList< T > &input)
Return the (stable) sort order for the list.
cellSet.H
Foam::polyMesh::globalData
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1295
Foam::point
vector point
Point is a vector.
Definition: point.H:43
Foam::hexCellLooper
Implementation of cellLooper.
Definition: hexCellLooper.H:62
Foam::fvMeshDistribute
Sends/receives parts of mesh+fvfields to neighbouring processors. Used in load balancing.
Definition: fvMeshDistribute.H:70
hexCellLooper.H
OBJstream.H
Foam::polyMesh::faceNeighbour
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1113
Foam::meshRefinement::MESH
Definition: meshRefinement.H:94
Foam::UPstream::nProcs
static label nProcs(const label communicator=worldComm)
Number of processes in parallel run, and 1 for serial run.
Definition: UPstream.H:445
Foam::VectorSpace< Vector< scalar >, scalar, 3 >::nComponents
static constexpr direction nComponents
Number of components in this vector space.
Definition: VectorSpace.H:101
Foam::globalMeshData::nTotalFaces
label nTotalFaces() const noexcept
Return total number of faces in decomposed mesh. Not.
Definition: globalMeshData.H:365
y
scalar y
Definition: LISASMDCalcMethod1.H:14
Foam::cellCuts
Description of cuts across cells.
Definition: cellCuts.H:110
Foam::cos
dimensionedScalar cos(const dimensionedScalar &ds)
Definition: dimensionedScalar.C:265
Foam::labelUIndList
UIndirectList< label > labelUIndList
UIndirectList of labels.
Definition: UIndirectList.H:58
Foam::meshRefinement::debugType
debugType
Enumeration for what to debug. Used as a bit-pattern.
Definition: meshRefinement.H:92