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