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