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-2018 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "meshRefinement.H"
27 #include "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  const label celli = mesh_.cellTree().findInside(keepPoint);
354 
355  if (celli != -1)
356  {
357  // I am the processor that holds the keepPoint
358 
359  forAll(features_, featI)
360  {
361  const edgeMesh& featureMesh = features_[featI];
362  const label featureLevel = features_.levels()[featI][0];
363  const labelListList& pointEdges = featureMesh.pointEdges();
364 
365  // Find regions on edgeMesh
366  labelList edgeRegion;
367  label nRegions = featureMesh.regions(edgeRegion);
368 
369 
370  PackedBoolList regionVisited(nRegions);
371 
372 
373  // 1. Seed all 'knots' in edgeMesh
374 
375 
376  forAll(pointEdges, pointi)
377  {
378  if (pointEdges[pointi].size() != 2)
379  {
381  {
382  Pout<< "Adding particle from point:" << pointi
383  << " coord:" << featureMesh.points()[pointi]
384  << " since number of emanating edges:"
385  << pointEdges[pointi].size()
386  << endl;
387  }
388 
389  // Non-manifold point. Create particle.
390  startPointCloud.addParticle
391  (
392  new trackedParticle
393  (
394  mesh_,
395  keepPoint,
396  celli,
397  featureMesh.points()[pointi], // endpos
398  featureLevel, // level
399  featI, // featureMesh
400  pointi, // end point
401  -1 // feature edge
402  )
403  );
404 
405  // Mark
406  if (pointEdges[pointi].size() > 0)
407  {
408  label e0 = pointEdges[pointi][0];
409  label regionI = edgeRegion[e0];
410  regionVisited[regionI] = 1u;
411  }
412  }
413  }
414 
415 
416  // 2. Any regions that have not been visited at all? These can
417  // only be circular regions!
418  forAll(featureMesh.edges(), edgeI)
419  {
420  if (regionVisited.set(edgeRegion[edgeI], 1u))
421  {
422  const edge& e = featureMesh.edges()[edgeI];
423  label pointi = e.start();
425  {
426  Pout<< "Adding particle from point:" << pointi
427  << " coord:" << featureMesh.points()[pointi]
428  << " on circular region:" << edgeRegion[edgeI]
429  << endl;
430  }
431 
432  // Non-manifold point. Create particle.
433  startPointCloud.addParticle
434  (
435  new trackedParticle
436  (
437  mesh_,
438  keepPoint,
439  celli,
440  featureMesh.points()[pointi], // endpos
441  featureLevel, // level
442  featI, // featureMesh
443  pointi, // end point
444  -1 // feature edge
445  )
446  );
447  }
448  }
449  }
450  }
451  }
452 
453 
454  // Largest refinement level of any feature passed through
455  maxFeatureLevel = labelList(mesh_.nCells(), -1);
456 
457  // Whether edge has been visited.
458  List<PackedBoolList> featureEdgeVisited(features_.size());
459 
460  forAll(features_, featI)
461  {
462  featureEdgeVisited[featI].setSize(features_[featI].edges().size());
463  featureEdgeVisited[featI] = 0u;
464  }
465 
466  // Database to pass into trackedParticle::move
468  (
469  startPointCloud,
470  maxFeatureLevel,
471  featureEdgeVisited
472  );
473 
474 
475  // Track all particles to their end position (= starting feature point)
476  // Note that the particle might have started on a different processor
477  // so this will transfer across nicely until we can start tracking proper.
478  scalar maxTrackLen = 2.0*mesh_.bounds().mag();
479 
481  {
482  Pout<< "Tracking " << startPointCloud.size()
483  << " particles over distance " << maxTrackLen
484  << " to find the starting cell" << endl;
485  }
486  startPointCloud.move(startPointCloud, td, maxTrackLen);
487 
488 
489  // Reset levels
490  maxFeatureLevel = -1;
491  forAll(features_, featI)
492  {
493  featureEdgeVisited[featI] = 0u;
494  }
495 
496 
498  (
499  mesh_,
500  "featureCloud",
502  );
503 
504  if (debug&meshRefinement::FEATURESEEDS)
505  {
506  Pout<< "Constructing cloud for cell marking" << endl;
507  }
508 
509  forAllIter(Cloud<trackedParticle>, startPointCloud, iter)
510  {
511  const trackedParticle& startTp = iter();
512 
513  label featI = startTp.i();
514  label pointi = startTp.j();
515 
516  const edgeMesh& featureMesh = features_[featI];
517  const labelList& pEdges = featureMesh.pointEdges()[pointi];
518 
519  // Now shoot particles down all pEdges.
520  forAll(pEdges, pEdgeI)
521  {
522  label edgeI = pEdges[pEdgeI];
523 
524  if (featureEdgeVisited[featI].set(edgeI, 1u))
525  {
526  // Unvisited edge. Make the particle go to the other point
527  // on the edge.
528 
529  const edge& e = featureMesh.edges()[edgeI];
530  label otherPointi = e.otherVertex(pointi);
531 
532  trackedParticle* tp(new trackedParticle(startTp));
533  tp->start() = tp->position();
534  tp->end() = featureMesh.points()[otherPointi];
535  tp->j() = otherPointi;
536  tp->k() = edgeI;
537 
538  if (debug&meshRefinement::FEATURESEEDS)
539  {
540  Pout<< "Adding particle for point:" << pointi
541  << " coord:" << tp->position()
542  << " feature:" << featI
543  << " to track to:" << tp->end()
544  << endl;
545  }
546 
547  cloud.addParticle(tp);
548  }
549  }
550  }
551 
552  startPointCloud.clear();
553 
554 
555  while (true)
556  {
557  // Track all particles to their end position.
558  if (debug&meshRefinement::FEATURESEEDS)
559  {
560  Pout<< "Tracking " << cloud.size()
561  << " particles over distance " << maxTrackLen
562  << " to mark cells" << endl;
563  }
564  cloud.move(cloud, td, maxTrackLen);
565 
566  // Make particle follow edge.
567  forAllIter(Cloud<trackedParticle>, cloud, iter)
568  {
569  trackedParticle& tp = iter();
570 
571  label featI = tp.i();
572  label pointi = tp.j();
573 
574  const edgeMesh& featureMesh = features_[featI];
575  const labelList& pEdges = featureMesh.pointEdges()[pointi];
576 
577  // Particle now at pointi. Check connected edges to see which one
578  // we have to visit now.
579 
580  bool keepParticle = false;
581 
582  forAll(pEdges, i)
583  {
584  label edgeI = pEdges[i];
585 
586  if (featureEdgeVisited[featI].set(edgeI, 1u))
587  {
588  // Unvisited edge. Make the particle go to the other point
589  // on the edge.
590 
591  const edge& e = featureMesh.edges()[edgeI];
592  label otherPointi = e.otherVertex(pointi);
593 
594  tp.start() = tp.position();
595  tp.end() = featureMesh.points()[otherPointi];
596  tp.j() = otherPointi;
597  tp.k() = edgeI;
598  keepParticle = true;
599  break;
600  }
601  }
602 
603  if (!keepParticle)
604  {
605  // Particle at 'knot' where another particle already has been
606  // seeded. Delete particle.
607  cloud.deleteParticle(tp);
608  }
609  }
610 
611 
612  if (debug&meshRefinement::FEATURESEEDS)
613  {
614  Pout<< "Remaining particles " << cloud.size() << endl;
615  }
616 
617  if (returnReduce(cloud.size(), sumOp<label>()) == 0)
618  {
619  break;
620  }
621  }
622 
623 
624 
625  // if (debug&meshRefinement::FEATURESEEDS)
626  //{
627  // forAll(maxFeatureLevel, celli)
628  // {
629  // if (maxFeatureLevel[celli] != -1)
630  // {
631  // Pout<< "Feature went through cell:" << celli
632  // << " coord:" << mesh_.cellCentres()[celli]
633  // << " leve:" << maxFeatureLevel[celli]
634  // << endl;
635  // }
636  // }
637  //}
638 }
639 
640 
641 // Calculates list of cells to refine based on intersection with feature edge.
642 Foam::label Foam::meshRefinement::markFeatureRefinement
643 (
644  const pointField& keepPoints,
645  const label nAllowRefine,
646 
647  labelList& refineCell,
648  label& nRefine
649 ) const
650 {
651  // Largest refinement level of any feature passed through
652  labelList maxFeatureLevel;
653  markFeatureCellLevel(keepPoints, maxFeatureLevel);
654 
655  // See which cells to refine. maxFeatureLevel will hold highest level
656  // of any feature edge that passed through.
657 
658  const labelList& cellLevel = meshCutter_.cellLevel();
659 
660  label oldNRefine = nRefine;
661 
662  forAll(maxFeatureLevel, celli)
663  {
664  if (maxFeatureLevel[celli] > cellLevel[celli])
665  {
666  // Mark
667  if
668  (
669  !markForRefine
670  (
671  0, // surface (n/a)
672  nAllowRefine,
673  refineCell[celli],
674  nRefine
675  )
676  )
677  {
678  // Reached limit
679  break;
680  }
681  }
682  }
683 
684  if
685  (
686  returnReduce(nRefine, sumOp<label>())
687  > returnReduce(nAllowRefine, sumOp<label>())
688  )
689  {
690  Info<< "Reached refinement limit." << endl;
691  }
692 
693  return returnReduce(nRefine-oldNRefine, sumOp<label>());
694 }
695 
696 
697 // Mark cells for distance-to-feature based refinement.
698 Foam::label Foam::meshRefinement::markInternalDistanceToFeatureRefinement
699 (
700  const label nAllowRefine,
701 
702  labelList& refineCell,
703  label& nRefine
704 ) const
705 {
706  const labelList& cellLevel = meshCutter_.cellLevel();
707  const pointField& cellCentres = mesh_.cellCentres();
708 
709  // Detect if there are any distance shells
710  if (features_.maxDistance() <= 0.0)
711  {
712  return 0;
713  }
714 
715  label oldNRefine = nRefine;
716 
717  // Collect cells to test
718  pointField testCc(cellLevel.size()-nRefine);
719  labelList testLevels(cellLevel.size()-nRefine);
720  label testI = 0;
721 
722  forAll(cellLevel, celli)
723  {
724  if (refineCell[celli] == -1)
725  {
726  testCc[testI] = cellCentres[celli];
727  testLevels[testI] = cellLevel[celli];
728  testI++;
729  }
730  }
731 
732  // Do test to see whether cells is inside/outside shell with higher level
733  labelList maxLevel;
734  features_.findHigherLevel(testCc, testLevels, maxLevel);
735 
736  // Mark for refinement. Note that we didn't store the original cellID so
737  // now just reloop in same order.
738  testI = 0;
739  forAll(cellLevel, celli)
740  {
741  if (refineCell[celli] == -1)
742  {
743  if (maxLevel[testI] > testLevels[testI])
744  {
745  bool reachedLimit = !markForRefine
746  (
747  maxLevel[testI], // mark with any positive value
748  nAllowRefine,
749  refineCell[celli],
750  nRefine
751  );
752 
753  if (reachedLimit)
754  {
755  if (debug)
756  {
757  Pout<< "Stopped refining internal cells"
758  << " since reaching my cell limit of "
759  << mesh_.nCells()+7*nRefine << endl;
760  }
761  break;
762  }
763  }
764  testI++;
765  }
766  }
767 
768  if
769  (
770  returnReduce(nRefine, sumOp<label>())
771  > returnReduce(nAllowRefine, sumOp<label>())
772  )
773  {
774  Info<< "Reached refinement limit." << endl;
775  }
776 
777  return returnReduce(nRefine-oldNRefine, sumOp<label>());
778 }
779 
780 
781 // Mark cells for non-surface intersection based refinement.
782 Foam::label Foam::meshRefinement::markInternalRefinement
783 (
784  const label nAllowRefine,
785 
786  labelList& refineCell,
787  label& nRefine
788 ) const
789 {
790  const labelList& cellLevel = meshCutter_.cellLevel();
791  const pointField& cellCentres = mesh_.cellCentres();
792 
793  label oldNRefine = nRefine;
794 
795  // Collect cells to test
796  pointField testCc(cellLevel.size()-nRefine);
797  labelList testLevels(cellLevel.size()-nRefine);
798  label testI = 0;
799 
800  forAll(cellLevel, celli)
801  {
802  if (refineCell[celli] == -1)
803  {
804  testCc[testI] = cellCentres[celli];
805  testLevels[testI] = cellLevel[celli];
806  testI++;
807  }
808  }
809 
810  // Do test to see whether cells is inside/outside shell with higher level
811  labelList maxLevel;
812  shells_.findHigherLevel(testCc, testLevels, maxLevel);
813 
814  // Mark for refinement. Note that we didn't store the original cellID so
815  // now just reloop in same order.
816  testI = 0;
817  forAll(cellLevel, celli)
818  {
819  if (refineCell[celli] == -1)
820  {
821  if (maxLevel[testI] > testLevels[testI])
822  {
823  bool reachedLimit = !markForRefine
824  (
825  maxLevel[testI], // mark with any positive value
826  nAllowRefine,
827  refineCell[celli],
828  nRefine
829  );
830 
831  if (reachedLimit)
832  {
833  if (debug)
834  {
835  Pout<< "Stopped refining internal cells"
836  << " since reaching my cell limit of "
837  << mesh_.nCells()+7*nRefine << endl;
838  }
839  break;
840  }
841  }
842  testI++;
843  }
844  }
845 
846  if
847  (
848  returnReduce(nRefine, sumOp<label>())
849  > returnReduce(nAllowRefine, sumOp<label>())
850  )
851  {
852  Info<< "Reached refinement limit." << endl;
853  }
854 
855  return returnReduce(nRefine-oldNRefine, sumOp<label>());
856 }
857 
858 
859 // Collect faces that are intersected and whose neighbours aren't yet marked
860 // for refinement.
861 Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
862 (
863  const labelList& refineCell
864 ) const
865 {
866  labelList testFaces(mesh_.nFaces());
867 
868  label nTest = 0;
869 
870  forAll(surfaceIndex_, facei)
871  {
872  if (surfaceIndex_[facei] != -1)
873  {
874  label own = mesh_.faceOwner()[facei];
875 
876  if (mesh_.isInternalFace(facei))
877  {
878  label nei = mesh_.faceNeighbour()[facei];
879 
880  if (refineCell[own] == -1 || refineCell[nei] == -1)
881  {
882  testFaces[nTest++] = facei;
883  }
884  }
885  else
886  {
887  if (refineCell[own] == -1)
888  {
889  testFaces[nTest++] = facei;
890  }
891  }
892  }
893  }
894  testFaces.setSize(nTest);
895 
896  return testFaces;
897 }
898 
899 
900 // Mark cells for surface intersection based refinement.
901 Foam::label Foam::meshRefinement::markSurfaceRefinement
902 (
903  const label nAllowRefine,
904  const labelList& neiLevel,
905  const pointField& neiCc,
906 
907  labelList& refineCell,
908  label& nRefine
909 ) const
910 {
911  const labelList& cellLevel = meshCutter_.cellLevel();
912  const pointField& cellCentres = mesh_.cellCentres();
913 
914  label oldNRefine = nRefine;
915 
916  // Use cached surfaceIndex_ to detect if any intersection. If so
917  // re-intersect to determine level wanted.
918 
919  // Collect candidate faces
920  // ~~~~~~~~~~~~~~~~~~~~~~~
921 
922  labelList testFaces(getRefineCandidateFaces(refineCell));
923 
924  // Collect segments
925  // ~~~~~~~~~~~~~~~~
926 
927  pointField start(testFaces.size());
928  pointField end(testFaces.size());
929  labelList minLevel(testFaces.size());
930 
931  forAll(testFaces, i)
932  {
933  label facei = testFaces[i];
934 
935  label own = mesh_.faceOwner()[facei];
936 
937  if (mesh_.isInternalFace(facei))
938  {
939  label nei = mesh_.faceNeighbour()[facei];
940 
941  start[i] = cellCentres[own];
942  end[i] = cellCentres[nei];
943  minLevel[i] = min(cellLevel[own], cellLevel[nei]);
944  }
945  else
946  {
947  label bFacei = facei - mesh_.nInternalFaces();
948 
949  start[i] = cellCentres[own];
950  end[i] = neiCc[bFacei];
951  minLevel[i] = min(cellLevel[own], neiLevel[bFacei]);
952  }
953  }
954 
955  // Extend segments a bit
956  {
957  const vectorField smallVec(rootSmall*(end-start));
958  start -= smallVec;
959  end += smallVec;
960  }
961 
962 
963  // Do test for higher intersections
964  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
965 
966  labelList surfaceHit;
967  labelList surfaceMinLevel;
968  surfaces_.findHigherIntersection
969  (
970  start,
971  end,
972  minLevel,
973 
974  surfaceHit,
975  surfaceMinLevel
976  );
977 
978 
979  // Mark cells for refinement
980  // ~~~~~~~~~~~~~~~~~~~~~~~~~
981 
982  forAll(testFaces, i)
983  {
984  label facei = testFaces[i];
985 
986  label surfI = surfaceHit[i];
987 
988  if (surfI != -1)
989  {
990  // Found intersection with surface with higher wanted
991  // refinement. Check if the level field on that surface
992  // specifies an even higher level. Note:this is weird. Should
993  // do the check with the surfaceMinLevel whilst intersecting the
994  // surfaces?
995 
996  label own = mesh_.faceOwner()[facei];
997 
998  if (surfaceMinLevel[i] > cellLevel[own])
999  {
1000  // Owner needs refining
1001  if
1002  (
1003  !markForRefine
1004  (
1005  surfI,
1006  nAllowRefine,
1007  refineCell[own],
1008  nRefine
1009  )
1010  )
1011  {
1012  break;
1013  }
1014  }
1015 
1016  if (mesh_.isInternalFace(facei))
1017  {
1018  label nei = mesh_.faceNeighbour()[facei];
1019  if (surfaceMinLevel[i] > cellLevel[nei])
1020  {
1021  // Neighbour needs refining
1022  if
1023  (
1024  !markForRefine
1025  (
1026  surfI,
1027  nAllowRefine,
1028  refineCell[nei],
1029  nRefine
1030  )
1031  )
1032  {
1033  break;
1034  }
1035  }
1036  }
1037  }
1038  }
1039 
1040  if
1041  (
1042  returnReduce(nRefine, sumOp<label>())
1043  > returnReduce(nAllowRefine, sumOp<label>())
1044  )
1045  {
1046  Info<< "Reached refinement limit." << endl;
1047  }
1048 
1049  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1050 }
1051 
1052 
1053 // Count number of matches of first argument in second argument
1054 Foam::label Foam::meshRefinement::countMatches
1055 (
1056  const List<point>& normals1,
1057  const List<point>& normals2,
1058  const scalar tol
1059 ) const
1060 {
1061  label nMatches = 0;
1062 
1063  forAll(normals1, i)
1064  {
1065  const vector& n1 = normals1[i];
1066 
1067  forAll(normals2, j)
1068  {
1069  const vector& n2 = normals2[j];
1070 
1071  if (magSqr(n1-n2) < tol)
1072  {
1073  nMatches++;
1074  break;
1075  }
1076  }
1077  }
1078 
1079  return nMatches;
1080 }
1081 
1082 
1083 // Mark cells for surface curvature based refinement.
1084 Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
1085 (
1086  const scalar curvature,
1087  const label nAllowRefine,
1088  const labelList& neiLevel,
1089  const pointField& neiCc,
1090 
1091  labelList& refineCell,
1092  label& nRefine
1093 ) const
1094 {
1095  const labelList& cellLevel = meshCutter_.cellLevel();
1096  const pointField& cellCentres = mesh_.cellCentres();
1097 
1098  label oldNRefine = nRefine;
1099 
1100  // 1. local test: any cell on more than one surface gets refined
1101  // (if its current level is < max of the surface max level)
1102 
1103  // 2. 'global' test: any cell on only one surface with a neighbour
1104  // on a different surface gets refined (if its current level etc.)
1105 
1106 
1107  const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
1108 
1109 
1110  // Collect candidate faces (i.e. intersecting any surface and
1111  // owner/neighbour not yet refined.
1112  labelList testFaces(getRefineCandidateFaces(refineCell));
1113 
1114  // Collect segments
1115  pointField start(testFaces.size());
1116  pointField end(testFaces.size());
1117  labelList minLevel(testFaces.size());
1118 
1119  forAll(testFaces, i)
1120  {
1121  label facei = testFaces[i];
1122 
1123  label own = mesh_.faceOwner()[facei];
1124 
1125  if (mesh_.isInternalFace(facei))
1126  {
1127  label nei = mesh_.faceNeighbour()[facei];
1128 
1129  start[i] = cellCentres[own];
1130  end[i] = cellCentres[nei];
1131  minLevel[i] = min(cellLevel[own], cellLevel[nei]);
1132  }
1133  else
1134  {
1135  label bFacei = facei - mesh_.nInternalFaces();
1136 
1137  start[i] = cellCentres[own];
1138  end[i] = neiCc[bFacei];
1139 
1140  if (!isMasterFace[facei])
1141  {
1142  Swap(start[i], end[i]);
1143  }
1144 
1145  minLevel[i] = min(cellLevel[own], neiLevel[bFacei]);
1146  }
1147  }
1148 
1149  // Extend segments a bit
1150  {
1151  const vectorField smallVec(rootSmall*(end-start));
1152  start -= smallVec;
1153  end += smallVec;
1154  }
1155 
1156 
1157  // Test for all intersections (with surfaces of higher max level than
1158  // minLevel) and cache per cell the interesting inter
1159  labelListList cellSurfLevels(mesh_.nCells());
1160  List<vectorList> cellSurfNormals(mesh_.nCells());
1161 
1162  {
1163  // Per segment the normals of the surfaces
1164  List<vectorList> surfaceNormal;
1165  // Per segment the list of levels of the surfaces
1166  labelListList surfaceLevel;
1167 
1168  surfaces_.findAllHigherIntersections
1169  (
1170  start,
1171  end,
1172  minLevel, // max level of surface has to be bigger
1173  // than min level of neighbouring cells
1174 
1175  surfaces_.maxLevel(),
1176 
1177  surfaceNormal,
1178  surfaceLevel
1179  );
1180 
1181 
1182  // Sort the data according to surface location. This will guarantee
1183  // that on coupled faces both sides visit the intersections in
1184  // the same order so will decide the same
1185  labelList visitOrder;
1186  forAll(surfaceNormal, pointi)
1187  {
1188  vectorList& pNormals = surfaceNormal[pointi];
1189  labelList& pLevel = surfaceLevel[pointi];
1190 
1191  sortedOrder(pNormals, visitOrder, normalLess(pNormals));
1192 
1193  pNormals = List<point>(pNormals, visitOrder);
1194  pLevel = UIndirectList<label>(pLevel, visitOrder);
1195  }
1196 
1197  // Clear out unnecessary data
1198  start.clear();
1199  end.clear();
1200  minLevel.clear();
1201 
1202  // Convert face-wise data to cell.
1203  forAll(testFaces, i)
1204  {
1205  label facei = testFaces[i];
1206  label own = mesh_.faceOwner()[facei];
1207 
1208  const vectorList& fNormals = surfaceNormal[i];
1209  const labelList& fLevels = surfaceLevel[i];
1210 
1211  forAll(fNormals, hitI)
1212  {
1213  if (fLevels[hitI] > cellLevel[own])
1214  {
1215  cellSurfLevels[own].append(fLevels[hitI]);
1216  cellSurfNormals[own].append(fNormals[hitI]);
1217  }
1218 
1219  if (mesh_.isInternalFace(facei))
1220  {
1221  label nei = mesh_.faceNeighbour()[facei];
1222  if (fLevels[hitI] > cellLevel[nei])
1223  {
1224  cellSurfLevels[nei].append(fLevels[hitI]);
1225  cellSurfNormals[nei].append(fNormals[hitI]);
1226  }
1227  }
1228  }
1229  }
1230  }
1231 
1232 
1233 
1234  // Bit of statistics
1235  if (debug)
1236  {
1237  label nSet = 0;
1238  label nNormals = 9;
1239  forAll(cellSurfNormals, celli)
1240  {
1241  const vectorList& normals = cellSurfNormals[celli];
1242  if (normals.size())
1243  {
1244  nSet++;
1245  nNormals += normals.size();
1246  }
1247  }
1248  reduce(nSet, sumOp<label>());
1249  reduce(nNormals, sumOp<label>());
1250  Info<< "markSurfaceCurvatureRefinement :"
1251  << " cells:" << mesh_.globalData().nTotalCells()
1252  << " of which with normals:" << nSet
1253  << " ; total normals stored:" << nNormals
1254  << endl;
1255  }
1256 
1257 
1258 
1259  bool reachedLimit = false;
1260 
1261 
1262  // 1. Check normals per cell
1263  // ~~~~~~~~~~~~~~~~~~~~~~~~~
1264 
1265  for
1266  (
1267  label celli = 0;
1268  !reachedLimit && celli < cellSurfNormals.size();
1269  celli++
1270  )
1271  {
1272  const vectorList& normals = cellSurfNormals[celli];
1273  const labelList& levels = cellSurfLevels[celli];
1274 
1275  // n^2 comparison of all normals in a cell
1276  for (label i = 0; !reachedLimit && i < normals.size(); i++)
1277  {
1278  for (label j = i+1; !reachedLimit && j < normals.size(); j++)
1279  {
1280  if ((normals[i] & normals[j]) < curvature)
1281  {
1282  label maxLevel = max(levels[i], levels[j]);
1283 
1284  if (cellLevel[celli] < maxLevel)
1285  {
1286  if
1287  (
1288  !markForRefine
1289  (
1290  maxLevel,
1291  nAllowRefine,
1292  refineCell[celli],
1293  nRefine
1294  )
1295  )
1296  {
1297  if (debug)
1298  {
1299  Pout<< "Stopped refining since reaching my cell"
1300  << " limit of " << mesh_.nCells()+7*nRefine
1301  << endl;
1302  }
1303  reachedLimit = true;
1304  break;
1305  }
1306  }
1307  }
1308  }
1309  }
1310  }
1311 
1312 
1313 
1314  // 2. Find out a measure of surface curvature
1315  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1316  // Look at normals between neighbouring surfaces
1317  // Loop over all faces. Could only be checkFaces, except if they're coupled
1318 
1319  // Internal faces
1320  for
1321  (
1322  label facei = 0;
1323  !reachedLimit && facei < mesh_.nInternalFaces();
1324  facei++
1325  )
1326  {
1327  label own = mesh_.faceOwner()[facei];
1328  label nei = mesh_.faceNeighbour()[facei];
1329 
1330  const vectorList& ownNormals = cellSurfNormals[own];
1331  const labelList& ownLevels = cellSurfLevels[own];
1332  const vectorList& neiNormals = cellSurfNormals[nei];
1333  const labelList& neiLevels = cellSurfLevels[nei];
1334 
1335 
1336  // Special case: owner normals set is a subset of the neighbour
1337  // normals. Do not do curvature refinement since other cell's normals
1338  // do not add any information. Avoids weird corner extensions of
1339  // refinement regions.
1340  bool ownIsSubset =
1341  countMatches(ownNormals, neiNormals)
1342  == ownNormals.size();
1343 
1344  bool neiIsSubset =
1345  countMatches(neiNormals, ownNormals)
1346  == neiNormals.size();
1347 
1348 
1349  if (!ownIsSubset && !neiIsSubset)
1350  {
1351  // n^2 comparison of between ownNormals and neiNormals
1352  for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
1353  {
1354  for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
1355  {
1356  // Have valid data on both sides. Check curvature.
1357  if ((ownNormals[i] & neiNormals[j]) < curvature)
1358  {
1359  // See which side to refine.
1360  if (cellLevel[own] < ownLevels[i])
1361  {
1362  if
1363  (
1364  !markForRefine
1365  (
1366  ownLevels[i],
1367  nAllowRefine,
1368  refineCell[own],
1369  nRefine
1370  )
1371  )
1372  {
1373  if (debug)
1374  {
1375  Pout<< "Stopped refining since reaching"
1376  << " my cell limit of "
1377  << mesh_.nCells()+7*nRefine << endl;
1378  }
1379  reachedLimit = true;
1380  break;
1381  }
1382  }
1383  if (cellLevel[nei] < neiLevels[j])
1384  {
1385  if
1386  (
1387  !markForRefine
1388  (
1389  neiLevels[j],
1390  nAllowRefine,
1391  refineCell[nei],
1392  nRefine
1393  )
1394  )
1395  {
1396  if (debug)
1397  {
1398  Pout<< "Stopped refining since reaching"
1399  << " my cell limit of "
1400  << mesh_.nCells()+7*nRefine << endl;
1401  }
1402  reachedLimit = true;
1403  break;
1404  }
1405  }
1406  }
1407  }
1408  }
1409  }
1410  }
1411 
1412 
1413  // Send over surface normal to neighbour cell.
1414  List<vectorList> neiSurfaceNormals;
1415  syncTools::swapBoundaryCellList(mesh_, cellSurfNormals, neiSurfaceNormals);
1416 
1417  // Boundary faces
1418  for
1419  (
1420  label facei = mesh_.nInternalFaces();
1421  !reachedLimit && facei < mesh_.nFaces();
1422  facei++
1423  )
1424  {
1425  label own = mesh_.faceOwner()[facei];
1426  label bFacei = facei - mesh_.nInternalFaces();
1427 
1428  const vectorList& ownNormals = cellSurfNormals[own];
1429  const labelList& ownLevels = cellSurfLevels[own];
1430  const vectorList& neiNormals = neiSurfaceNormals[bFacei];
1431 
1432  // Special case: owner normals set is a subset of the neighbour
1433  // normals. Do not do curvature refinement since other cell's normals
1434  // do not add any information. Avoids weird corner extensions of
1435  // refinement regions.
1436  bool ownIsSubset =
1437  countMatches(ownNormals, neiNormals)
1438  == ownNormals.size();
1439 
1440  bool neiIsSubset =
1441  countMatches(neiNormals, ownNormals)
1442  == neiNormals.size();
1443 
1444 
1445  if (!ownIsSubset && !neiIsSubset)
1446  {
1447  // n^2 comparison of between ownNormals and neiNormals
1448  for (label i = 0; !reachedLimit && i < ownNormals.size(); i++)
1449  {
1450  for (label j = 0; !reachedLimit && j < neiNormals.size(); j++)
1451  {
1452  // Have valid data on both sides. Check curvature.
1453  if ((ownNormals[i] & neiNormals[j]) < curvature)
1454  {
1455  if (cellLevel[own] < ownLevels[i])
1456  {
1457  if
1458  (
1459  !markForRefine
1460  (
1461  ownLevels[i],
1462  nAllowRefine,
1463  refineCell[own],
1464  nRefine
1465  )
1466  )
1467  {
1468  if (debug)
1469  {
1470  Pout<< "Stopped refining since reaching"
1471  << " my cell limit of "
1472  << mesh_.nCells()+7*nRefine
1473  << endl;
1474  }
1475  reachedLimit = true;
1476  break;
1477  }
1478  }
1479  }
1480  }
1481  }
1482  }
1483  }
1484 
1485 
1486  if
1487  (
1488  returnReduce(nRefine, sumOp<label>())
1489  > returnReduce(nAllowRefine, sumOp<label>())
1490  )
1491  {
1492  Info<< "Reached refinement limit." << endl;
1493  }
1494 
1495  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1496 }
1497 
1498 
1501  const scalar planarCos,
1502  const vector& point0,
1503  const vector& normal0,
1504 
1505  const vector& point1,
1506  const vector& normal1
1507 ) const
1508 {
1509  //- Hits differ and angles are oppositeish and
1510  // hits have a normal distance
1511  vector d = point1-point0;
1512  scalar magD = mag(d);
1513 
1514  if (magD > mergeDistance())
1515  {
1516  scalar cosAngle = (normal0 & normal1);
1517 
1518  vector avg = Zero;
1519  if (cosAngle < (-1+planarCos))
1520  {
1521  // Opposite normals
1522  avg = 0.5*(normal0-normal1);
1523  }
1524  else if (cosAngle > (1-planarCos))
1525  {
1526  avg = 0.5*(normal0+normal1);
1527  }
1528 
1529  if (avg != vector::zero)
1530  {
1531  avg /= mag(avg);
1532 
1533  // Check normal distance of intersection locations
1534  if (mag(avg&d) > mergeDistance())
1535  {
1536  return true;
1537  }
1538  else
1539  {
1540  return false;
1541  }
1542  }
1543  else
1544  {
1545  return false;
1546  }
1547  }
1548  else
1549  {
1550  return false;
1551  }
1552 }
1553 
1554 
1555 // Mark small gaps
1558  const scalar planarCos,
1559  const vector& point0,
1560  const vector& normal0,
1561 
1562  const vector& point1,
1563  const vector& normal1
1564 ) const
1565 {
1566  //- Hits differ and angles are oppositeish and
1567  // hits have a normal distance
1568  vector d = point1-point0;
1569  scalar magD = mag(d);
1570 
1571  if (magD > mergeDistance())
1572  {
1573  scalar cosAngle = (normal0 & normal1);
1574 
1575  vector avg = Zero;
1576  if (cosAngle < (-1+planarCos))
1577  {
1578  // Opposite normals
1579  avg = 0.5*(normal0-normal1);
1580  }
1581  else if (cosAngle > (1-planarCos))
1582  {
1583  avg = 0.5*(normal0+normal1);
1584  }
1585 
1586  if (avg != vector::zero)
1587  {
1588  avg /= mag(avg);
1589  d /= magD;
1590 
1591  // Check average normal with respect to intersection locations
1592  if (mag(avg&d) > Foam::cos(degToRad(45)))
1593  {
1594  return true;
1595  }
1596  else
1597  {
1598  return false;
1599  }
1600  }
1601  else
1602  {
1603  return false;
1604  }
1605  }
1606  else
1607  {
1608  return false;
1609  }
1610 }
1611 
1612 
1613 bool Foam::meshRefinement::checkProximity
1614 (
1615  const scalar planarCos,
1616  const label nAllowRefine,
1617 
1618  const label surfaceLevel, // current intersection max level
1619  const vector& surfaceLocation, // current intersection location
1620  const vector& surfaceNormal, // current intersection normal
1621 
1622  const label celli,
1623 
1624  label& cellMaxLevel, // cached max surface level for this cell
1625  vector& cellMaxLocation, // cached surface normal for this cell
1626  vector& cellMaxNormal, // cached surface normal for this cell
1627 
1628  labelList& refineCell,
1629  label& nRefine
1630 ) const
1631 {
1632  const labelList& cellLevel = meshCutter_.cellLevel();
1633 
1634  // Test if surface applicable
1635  if (surfaceLevel > cellLevel[celli])
1636  {
1637  if (cellMaxLevel == -1)
1638  {
1639  // First visit of cell. Store
1640  cellMaxLevel = surfaceLevel;
1641  cellMaxLocation = surfaceLocation;
1642  cellMaxNormal = surfaceNormal;
1643  }
1644  else
1645  {
1646  // Second or more visit.
1647  // Check if
1648  // - different location
1649  // - opposite surface
1650 
1651  bool closeSurfaces = isNormalGap
1652  (
1653  planarCos,
1654  cellMaxLocation,
1655  cellMaxNormal,
1656  surfaceLocation,
1657  surfaceNormal
1658  );
1659 
1660  // Set normal to that of highest surface. Not really necessary
1661  // over here but we reuse cellMax info when doing coupled faces.
1662  if (surfaceLevel > cellMaxLevel)
1663  {
1664  cellMaxLevel = surfaceLevel;
1665  cellMaxLocation = surfaceLocation;
1666  cellMaxNormal = surfaceNormal;
1667  }
1668 
1669 
1670  if (closeSurfaces)
1671  {
1672  // Pout<< "Found gap:" << nl
1673  // << " location:" << surfaceLocation
1674  // << "\tnormal :" << surfaceNormal << nl
1676  // << "\tnormal :" << cellMaxNormal << nl
1677  // << "\tcos :" << (surfaceNormal&cellMaxNormal) << nl
1678  // << endl;
1679 
1680  return markForRefine
1681  (
1682  surfaceLevel, // mark with any non-neg number.
1683  nAllowRefine,
1684  refineCell[celli],
1685  nRefine
1686  );
1687  }
1688  }
1689  }
1690 
1691  // Did not reach refinement limit.
1692  return true;
1693 }
1694 
1695 
1696 Foam::label Foam::meshRefinement::markProximityRefinement
1697 (
1698  const scalar planarCos,
1699  const label nAllowRefine,
1700  const labelList& neiLevel,
1701  const pointField& neiCc,
1702 
1703  labelList& refineCell,
1704  label& nRefine
1705 ) const
1706 {
1707  const labelList& cellLevel = meshCutter_.cellLevel();
1708  const pointField& cellCentres = mesh_.cellCentres();
1709 
1710  label oldNRefine = nRefine;
1711 
1712  // 1. local test: any cell on more than one surface gets refined
1713  // (if its current level is < max of the surface max level)
1714 
1715  // 2. 'global' test: any cell on only one surface with a neighbour
1716  // on a different surface gets refined (if its current level etc.)
1717 
1718 
1719  // Collect candidate faces (i.e. intersecting any surface and
1720  // owner/neighbour not yet refined.
1721  labelList testFaces(getRefineCandidateFaces(refineCell));
1722 
1723  // Collect segments
1724  pointField start(testFaces.size());
1725  pointField end(testFaces.size());
1726  labelList minLevel(testFaces.size());
1727 
1728  forAll(testFaces, i)
1729  {
1730  label facei = testFaces[i];
1731 
1732  label own = mesh_.faceOwner()[facei];
1733 
1734  if (mesh_.isInternalFace(facei))
1735  {
1736  label nei = mesh_.faceNeighbour()[facei];
1737 
1738  start[i] = cellCentres[own];
1739  end[i] = cellCentres[nei];
1740  minLevel[i] = min(cellLevel[own], cellLevel[nei]);
1741  }
1742  else
1743  {
1744  label bFacei = facei - mesh_.nInternalFaces();
1745 
1746  start[i] = cellCentres[own];
1747  end[i] = neiCc[bFacei];
1748  minLevel[i] = min(cellLevel[own], neiLevel[bFacei]);
1749  }
1750  }
1751 
1752  // Extend segments a bit
1753  {
1754  const vectorField smallVec(rootSmall*(end-start));
1755  start -= smallVec;
1756  end += smallVec;
1757  }
1758 
1759 
1760  // Test for all intersections (with surfaces of higher gap level than
1761  // minLevel) and cache per cell the max surface level and the local normal
1762  // on that surface.
1763  labelList cellMaxLevel(mesh_.nCells(), -1);
1764  vectorField cellMaxNormal(mesh_.nCells(), Zero);
1765  pointField cellMaxLocation(mesh_.nCells(), Zero);
1766 
1767  {
1768  // Per segment the normals of the surfaces
1769  List<vectorList> surfaceLocation;
1770  List<vectorList> surfaceNormal;
1771  // Per segment the list of levels of the surfaces
1772  labelListList surfaceLevel;
1773 
1774  surfaces_.findAllHigherIntersections
1775  (
1776  start,
1777  end,
1778  minLevel, // gap level of surface has to be bigger
1779  // than min level of neighbouring cells
1780 
1781  surfaces_.gapLevel(),
1782 
1783  surfaceLocation,
1784  surfaceNormal,
1785  surfaceLevel
1786  );
1787  // Clear out unnecessary data
1788  start.clear();
1789  end.clear();
1790  minLevel.clear();
1791 
1793  // OBJstream str
1794  //(
1795  // mesh_.time().path()
1796  // / "findAllHigherIntersections_"
1797  // + mesh_.time().timeName()
1798  // + ".obj"
1799  //);
1801  // OBJstream str2
1802  //(
1803  // mesh_.time().path()
1804  // / "findAllHigherIntersections2_"
1805  // + mesh_.time().timeName()
1806  // + ".obj"
1807  //);
1808 
1809  forAll(testFaces, i)
1810  {
1811  label facei = testFaces[i];
1812  label own = mesh_.faceOwner()[facei];
1813 
1814  const labelList& fLevels = surfaceLevel[i];
1815  const vectorList& fPoints = surfaceLocation[i];
1816  const vectorList& fNormals = surfaceNormal[i];
1817 
1818  forAll(fLevels, hitI)
1819  {
1820  checkProximity
1821  (
1822  planarCos,
1823  nAllowRefine,
1824 
1825  fLevels[hitI],
1826  fPoints[hitI],
1827  fNormals[hitI],
1828 
1829  own,
1830  cellMaxLevel[own],
1831  cellMaxLocation[own],
1832  cellMaxNormal[own],
1833 
1834  refineCell,
1835  nRefine
1836  );
1837  }
1838 
1839  if (mesh_.isInternalFace(facei))
1840  {
1841  label nei = mesh_.faceNeighbour()[facei];
1842 
1843  forAll(fLevels, hitI)
1844  {
1845  checkProximity
1846  (
1847  planarCos,
1848  nAllowRefine,
1849 
1850  fLevels[hitI],
1851  fPoints[hitI],
1852  fNormals[hitI],
1853 
1854  nei,
1855  cellMaxLevel[nei],
1856  cellMaxLocation[nei],
1857  cellMaxNormal[nei],
1858 
1859  refineCell,
1860  nRefine
1861  );
1862  }
1863  }
1864  }
1865  }
1866 
1867  // 2. Find out a measure of surface curvature and region edges.
1868  // Send over surface region and surface normal to neighbour cell.
1869 
1870  labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1871  pointField neiBndMaxLocation(mesh_.nFaces()-mesh_.nInternalFaces());
1872  vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
1873 
1874  for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
1875  {
1876  label bFacei = facei-mesh_.nInternalFaces();
1877  label own = mesh_.faceOwner()[facei];
1878 
1879  neiBndMaxLevel[bFacei] = cellMaxLevel[own];
1880  neiBndMaxLocation[bFacei] = cellMaxLocation[own];
1881  neiBndMaxNormal[bFacei] = cellMaxNormal[own];
1882  }
1883  syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel);
1884  syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLocation);
1885  syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal);
1886 
1887  // Loop over all faces. Could only be checkFaces.. except if they're coupled
1888 
1889  // Internal faces
1890  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
1891  {
1892  label own = mesh_.faceOwner()[facei];
1893  label nei = mesh_.faceNeighbour()[facei];
1894 
1895  if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
1896  {
1897  // Have valid data on both sides. Check planarCos.
1898  if
1899  (
1900  isNormalGap
1901  (
1902  planarCos,
1903  cellMaxLocation[own],
1904  cellMaxNormal[own],
1905  cellMaxLocation[nei],
1906  cellMaxNormal[nei]
1907  )
1908  )
1909  {
1910  // See which side to refine
1911  if (cellLevel[own] < cellMaxLevel[own])
1912  {
1913  if
1914  (
1915  !markForRefine
1916  (
1917  cellMaxLevel[own],
1918  nAllowRefine,
1919  refineCell[own],
1920  nRefine
1921  )
1922  )
1923  {
1924  if (debug)
1925  {
1926  Pout<< "Stopped refining since reaching my cell"
1927  << " limit of " << mesh_.nCells()+7*nRefine
1928  << endl;
1929  }
1930  break;
1931  }
1932  }
1933 
1934  if (cellLevel[nei] < cellMaxLevel[nei])
1935  {
1936  if
1937  (
1938  !markForRefine
1939  (
1940  cellMaxLevel[nei],
1941  nAllowRefine,
1942  refineCell[nei],
1943  nRefine
1944  )
1945  )
1946  {
1947  if (debug)
1948  {
1949  Pout<< "Stopped refining since reaching my cell"
1950  << " limit of " << mesh_.nCells()+7*nRefine
1951  << endl;
1952  }
1953  break;
1954  }
1955  }
1956  }
1957  }
1958  }
1959  // Boundary faces
1960  for (label facei = mesh_.nInternalFaces(); facei < mesh_.nFaces(); facei++)
1961  {
1962  label own = mesh_.faceOwner()[facei];
1963  label bFacei = facei - mesh_.nInternalFaces();
1964 
1965  if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFacei] != -1)
1966  {
1967  // Have valid data on both sides. Check planarCos.
1968  if
1969  (
1970  isNormalGap
1971  (
1972  planarCos,
1973  cellMaxLocation[own],
1974  cellMaxNormal[own],
1975  neiBndMaxLocation[bFacei],
1976  neiBndMaxNormal[bFacei]
1977  )
1978  )
1979  {
1980  if
1981  (
1982  !markForRefine
1983  (
1984  cellMaxLevel[own],
1985  nAllowRefine,
1986  refineCell[own],
1987  nRefine
1988  )
1989  )
1990  {
1991  if (debug)
1992  {
1993  Pout<< "Stopped refining since reaching my cell"
1994  << " limit of " << mesh_.nCells()+7*nRefine
1995  << endl;
1996  }
1997  break;
1998  }
1999  }
2000  }
2001  }
2002 
2003  if
2004  (
2005  returnReduce(nRefine, sumOp<label>())
2006  > returnReduce(nAllowRefine, sumOp<label>())
2007  )
2008  {
2009  Info<< "Reached refinement limit." << endl;
2010  }
2011 
2012  return returnReduce(nRefine-oldNRefine, sumOp<label>());
2013 }
2014 
2015 
2016 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
2017 
2018 // Calculate list of cells to refine. Gets for any edge (start - end)
2019 // whether it intersects the surface. Does more accurate test and checks
2020 // the wanted level on the surface intersected.
2021 // Does approximate precalculation of how many cells can be refined before
2022 // hitting overall limit maxGlobalCells.
2025  const pointField& keepPoints,
2026  const scalar curvature,
2027  const scalar planarAngle,
2028 
2029  const bool featureRefinement,
2030  const bool featureDistanceRefinement,
2031  const bool internalRefinement,
2032  const bool surfaceRefinement,
2033  const bool curvatureRefinement,
2034  const bool gapRefinement,
2035  const label maxGlobalCells,
2036  const label maxLocalCells
2037 ) const
2038 {
2039  label totNCells = mesh_.globalData().nTotalCells();
2040 
2041  labelList cellsToRefine;
2042 
2043  if (totNCells >= maxGlobalCells)
2044  {
2045  Info<< "No cells marked for refinement since reached limit "
2046  << maxGlobalCells << '.' << endl;
2047  }
2048  else
2049  {
2050  // Every cell I refine adds 7 cells. Estimate the number of cells
2051  // I am allowed to refine.
2052  // Assume perfect distribution so can only refine as much the fraction
2053  // of the mesh I hold. This prediction step prevents us having to do
2054  // lots of reduces to keep count of the total number of cells selected
2055  // for refinement.
2056 
2057  // scalar fraction = scalar(mesh_.nCells())/totNCells;
2058  // label myMaxCells = label(maxGlobalCells*fraction);
2059  // label nAllowRefine = (myMaxCells - mesh_.nCells())/7;
2061  //
2062  // Pout<< "refineCandidates:" << nl
2063  // << " total cells:" << totNCells << nl
2064  // << " local cells:" << mesh_.nCells() << nl
2065  // << " local fraction:" << fraction << nl
2066  // << " local allowable cells:" << myMaxCells << nl
2067  // << " local allowable refinement:" << nAllowRefine << nl
2068  // << endl;
2069 
2070  //- Disable refinement shortcut. nAllowRefine is per processor limit.
2071  label nAllowRefine = labelMax / Pstream::nProcs();
2072 
2073  // Marked for refinement (>= 0) or not (-1). Actual value is the
2074  // index of the surface it intersects.
2075  labelList refineCell(mesh_.nCells(), -1);
2076  label nRefine = 0;
2077 
2078 
2079  // Swap neighbouring cell centres and cell level
2080  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
2081  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
2082  calcNeighbourData(neiLevel, neiCc);
2083 
2084 
2085 
2086  // Cells pierced by feature lines
2087  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2088 
2089  if (featureRefinement)
2090  {
2091  label nFeatures = markFeatureRefinement
2092  (
2093  keepPoints,
2094  nAllowRefine,
2095 
2096  refineCell,
2097  nRefine
2098  );
2099 
2100  Info<< "Marked for refinement due to explicit features "
2101  << ": " << nFeatures << " cells." << endl;
2102  }
2103 
2104  // Inside distance-to-feature shells
2105  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2106 
2107  if (featureDistanceRefinement)
2108  {
2109  label nShell = markInternalDistanceToFeatureRefinement
2110  (
2111  nAllowRefine,
2112 
2113  refineCell,
2114  nRefine
2115  );
2116  Info<< "Marked for refinement due to distance to explicit features "
2117  ": " << nShell << " cells." << endl;
2118  }
2119 
2120  // Inside refinement shells
2121  // ~~~~~~~~~~~~~~~~~~~~~~~~
2122 
2123  if (internalRefinement)
2124  {
2125  label nShell = markInternalRefinement
2126  (
2127  nAllowRefine,
2128 
2129  refineCell,
2130  nRefine
2131  );
2132  Info<< "Marked for refinement due to refinement shells "
2133  << ": " << nShell << " cells." << endl;
2134  }
2135 
2136  // Refinement based on intersection of surface
2137  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2138 
2139  if (surfaceRefinement)
2140  {
2141  label nSurf = markSurfaceRefinement
2142  (
2143  nAllowRefine,
2144  neiLevel,
2145  neiCc,
2146 
2147  refineCell,
2148  nRefine
2149  );
2150  Info<< "Marked for refinement due to surface intersection "
2151  << ": " << nSurf << " cells." << endl;
2152  }
2153 
2154  // Refinement based on curvature of surface
2155  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2156 
2157  if
2158  (
2159  curvatureRefinement
2160  && (curvature >= -1 && curvature <= 1)
2161  && (surfaces_.minLevel() != surfaces_.maxLevel())
2162  )
2163  {
2164  label nCurv = markSurfaceCurvatureRefinement
2165  (
2166  curvature,
2167  nAllowRefine,
2168  neiLevel,
2169  neiCc,
2170 
2171  refineCell,
2172  nRefine
2173  );
2174  Info<< "Marked for refinement due to curvature/regions "
2175  << ": " << nCurv << " cells." << endl;
2176  }
2177 
2178 
2179  const scalar planarCos = Foam::cos(degToRad(planarAngle));
2180 
2181  if
2182  (
2183  gapRefinement
2184  && (planarCos >= -1 && planarCos <= 1)
2185  && (max(surfaces_.gapLevel()) > -1)
2186  )
2187  {
2188  Info<< "Specified gap level : " << max(surfaces_.gapLevel())
2189  << ", planar angle " << planarAngle << endl;
2190 
2191  label nGap = markProximityRefinement
2192  (
2193  planarCos,
2194  nAllowRefine,
2195  neiLevel,
2196  neiCc,
2197 
2198  refineCell,
2199  nRefine
2200  );
2201  Info<< "Marked for refinement due to close opposite surfaces "
2202  << ": " << nGap << " cells." << endl;
2203  }
2204 
2205 
2206 
2207  // Pack cells-to-refine
2208  // ~~~~~~~~~~~~~~~~~~~~
2209 
2210  cellsToRefine.setSize(nRefine);
2211  nRefine = 0;
2212 
2213  forAll(refineCell, celli)
2214  {
2215  if (refineCell[celli] != -1)
2216  {
2217  cellsToRefine[nRefine++] = celli;
2218  }
2219  }
2220  }
2221 
2222  return cellsToRefine;
2223 }
2224 
2225 
2228  const labelList& cellsToRefine
2229 )
2230 {
2231  // Mesh changing engine.
2232  polyTopoChange meshMod(mesh_);
2233 
2234  // Play refinement commands into mesh changer.
2235  meshCutter_.setRefinement(cellsToRefine, meshMod);
2236 
2237  // Create mesh (no inflation), return map from old to new mesh.
2238  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false);
2239 
2240  // Update fields
2241  mesh_.updateMesh(map);
2242 
2243  // Optionally inflate mesh
2244  if (map().hasMotionPoints())
2245  {
2246  mesh_.movePoints(map().preMotionPoints());
2247  }
2248  else
2249  {
2250  // Delete mesh volumes.
2251  mesh_.clearOut();
2252  }
2253 
2254  // Reset the instance for if in overwrite mode
2255  mesh_.setInstance(timeName());
2256 
2257  // Update intersection info
2258  updateMesh(map, getChangedFaces(map, cellsToRefine));
2259 
2260  return map;
2261 }
2262 
2263 
2264 // Do refinement of consistent set of cells followed by truncation and
2265 // load balancing.
2269  const string& msg,
2270  decompositionMethod& decomposer,
2271  fvMeshDistribute& distributor,
2272  const labelList& cellsToRefine,
2273  const scalar maxLoadUnbalance
2274 )
2275 {
2276  // Do all refinement
2277  refine(cellsToRefine);
2278 
2279  if (debug&meshRefinement::MESH)
2280  {
2281  Pout<< "Writing refined but unbalanced " << msg
2282  << " mesh to time " << timeName() << endl;
2283  write
2284  (
2285  debugType(debug),
2286  writeType(writeLevel() | WRITEMESH),
2287  mesh_.time().path()/timeName()
2288  );
2289  Pout<< "Dumped debug data in = "
2290  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2291 
2292  // test all is still synced across proc patches
2293  checkData();
2294  }
2295 
2296  Info<< "Refined mesh in = "
2297  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2298  printMeshInfo(debug, "After refinement " + msg);
2299 
2300 
2301  // Load balancing
2302  // ~~~~~~~~~~~~~~
2303 
2305 
2306  if (Pstream::nProcs() > 1)
2307  {
2308  scalar nIdealCells =
2309  mesh_.globalData().nTotalCells()
2310  / Pstream::nProcs();
2311 
2312  scalar unbalance = returnReduce
2313  (
2314  mag(1.0-mesh_.nCells()/nIdealCells),
2315  maxOp<scalar>()
2316  );
2317 
2318  if (unbalance <= maxLoadUnbalance)
2319  {
2320  Info<< "Skipping balancing since max unbalance " << unbalance
2321  << " is less than allowable " << maxLoadUnbalance
2322  << endl;
2323  }
2324  else
2325  {
2326  scalarField cellWeights(mesh_.nCells(), 1);
2327 
2328  distMap = balance
2329  (
2330  false, // keepZoneFaces
2331  false, // keepBaffles
2332  cellWeights,
2333  decomposer,
2334  distributor
2335  );
2336 
2337  Info<< "Balanced mesh in = "
2338  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2339 
2340  printMeshInfo(debug, "After balancing " + msg);
2341 
2342 
2343  if (debug&meshRefinement::MESH)
2344  {
2345  Pout<< "Writing balanced " << msg
2346  << " mesh to time " << timeName() << endl;
2347  write
2348  (
2349  debugType(debug),
2350  writeType(writeLevel() | WRITEMESH),
2351  mesh_.time().path()/timeName()
2352  );
2353  Pout<< "Dumped debug data in = "
2354  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2355 
2356  // test all is still synced across proc patches
2357  checkData();
2358  }
2359  }
2360  }
2361 
2362  return distMap;
2363 }
2364 
2365 
2366 // Do load balancing followed by refinement of consistent set of cells.
2370  const string& msg,
2371  decompositionMethod& decomposer,
2372  fvMeshDistribute& distributor,
2373  const labelList& initCellsToRefine,
2374  const scalar maxLoadUnbalance
2375 )
2376 {
2377  labelList cellsToRefine(initCellsToRefine);
2378 
2379  //{
2380  // globalIndex globalCells(mesh_.nCells());
2381  //
2382  // Info<< "** Distribution before balancing/refining:" << endl;
2383  // for (label proci = 0; proci < Pstream::nProcs(); proci++)
2384  // {
2385  // Info<< " " << proci << '\t'
2386  // << globalCells.localSize(proci) << endl;
2387  // }
2388  // Info<< endl;
2389  //}
2390  //{
2391  // globalIndex globalCells(cellsToRefine.size());
2392  //
2393  // Info<< "** Cells to be refined:" << endl;
2394  // for (label proci = 0; proci < Pstream::nProcs(); proci++)
2395  // {
2396  // Info<< " " << proci << '\t'
2397  // << globalCells.localSize(proci) << endl;
2398  // }
2399  // Info<< endl;
2400  //}
2401 
2402 
2403  // Load balancing
2404  // ~~~~~~~~~~~~~~
2405 
2407 
2408  if (Pstream::nProcs() > 1)
2409  {
2410  // First check if we need to balance at all. Precalculate number of
2411  // cells after refinement and see what maximum difference is.
2412  scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
2413  scalar nIdealNewCells =
2414  returnReduce(nNewCells, sumOp<scalar>())
2415  / Pstream::nProcs();
2416  scalar unbalance = returnReduce
2417  (
2418  mag(1.0-nNewCells/nIdealNewCells),
2419  maxOp<scalar>()
2420  );
2421 
2422  if (unbalance <= maxLoadUnbalance)
2423  {
2424  Info<< "Skipping balancing since max unbalance " << unbalance
2425  << " is less than allowable " << maxLoadUnbalance
2426  << endl;
2427  }
2428  else
2429  {
2430  scalarField cellWeights(mesh_.nCells(), 1);
2431  forAll(cellsToRefine, i)
2432  {
2433  cellWeights[cellsToRefine[i]] += 7;
2434  }
2435 
2436  distMap = balance
2437  (
2438  false, // keepZoneFaces
2439  false, // keepBaffles
2440  cellWeights,
2441  decomposer,
2442  distributor
2443  );
2444 
2445  // Update cells to refine
2446  distMap().distributeCellIndices(cellsToRefine);
2447 
2448  Info<< "Balanced mesh in = "
2449  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2450  }
2451 
2452  //{
2453  // globalIndex globalCells(mesh_.nCells());
2454  //
2455  // Info<< "** Distribution after balancing:" << endl;
2456  // for (label proci = 0; proci < Pstream::nProcs(); proci++)
2457  // {
2458  // Info<< " " << proci << '\t'
2459  // << globalCells.localSize(proci) << endl;
2460  // }
2461  // Info<< endl;
2462  //}
2463 
2464  printMeshInfo(debug, "After balancing " + msg);
2465 
2466  if (debug&meshRefinement::MESH)
2467  {
2468  Pout<< "Writing balanced " << msg
2469  << " mesh to time " << timeName() << endl;
2470  write
2471  (
2472  debugType(debug),
2473  writeType(writeLevel() | WRITEMESH),
2474  mesh_.time().path()/timeName()
2475  );
2476  Pout<< "Dumped debug data in = "
2477  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2478 
2479  // test all is still synced across proc patches
2480  checkData();
2481  }
2482  }
2483 
2484 
2485  // Refinement
2486  // ~~~~~~~~~~
2487 
2488  refine(cellsToRefine);
2489 
2490  if (debug&meshRefinement::MESH)
2491  {
2492  Pout<< "Writing refined " << msg
2493  << " mesh to time " << timeName() << endl;
2494  write
2495  (
2496  debugType(debug),
2497  writeType(writeLevel() | WRITEMESH),
2498  mesh_.time().path()/timeName()
2499  );
2500  Pout<< "Dumped debug data in = "
2501  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2502 
2503  // test all is still synced across proc patches
2504  checkData();
2505  }
2506 
2507  Info<< "Refined mesh in = "
2508  << mesh_.time().cpuTimeIncrement() << " s" << endl;
2509 
2510  //{
2511  // globalIndex globalCells(mesh_.nCells());
2512  //
2513  // Info<< "** After refinement distribution:" << endl;
2514  // for (label proci = 0; proci < Pstream::nProcs(); proci++)
2515  // {
2516  // Info<< " " << proci << '\t'
2517  // << globalCells.localSize(proci) << endl;
2518  // }
2519  // Info<< endl;
2520  //}
2521 
2522  printMeshInfo(debug, "After refinement " + msg);
2523 
2524  return distMap;
2525 }
2526 
2527 
2528 // ************************************************************************* //
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())
Synchronize 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:295
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:163
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
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:151
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.
void write(Ostream &, const label, const dictionary &)
Write with dictionary lookup.
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.
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
void clear()
Clear the list, i.e. set size to zero.
Definition: ListI.H:124
A cloud is a collection of lagrangian particles.
Definition: cloud.H:51
const edgeList & edges() const
Return edges.
Definition: edgeMeshI.H:68
point & end()
Point to track to.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1169
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:192
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)
Synchronize 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:303
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:98
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