FreeFOAM The Cross-Platform CFD Toolkit
meshRefinementRefine.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 1991-2010 OpenCFD Ltd.
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"
28 #include <OpenFOAM/syncTools.H>
29 #include <OpenFOAM/Time.H>
31 #include <autoMesh/shellSurfaces.H>
32 #include <meshTools/faceSet.H>
38 #include <lagrangian/Cloud.H>
39 //#include <OpenFOAM/globalIndex.H>
40 
41 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42 
43 // Get faces (on the new mesh) that have in some way been affected by the
44 // mesh change. Picks up all faces but those that are between two
45 // unrefined faces. (Note that of an unchanged face the edge still might be
46 // split but that does not change any face centre or cell centre.
47 Foam::labelList Foam::meshRefinement::getChangedFaces
48 (
49  const mapPolyMesh& map,
50  const labelList& oldCellsToRefine
51 )
52 {
53  const polyMesh& mesh = map.mesh();
54 
55  labelList changedFaces;
56  {
57  // Mark any face on a cell which has been added or changed
58  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
59  // Note that refining a face changes the face centre (for a warped face)
60  // which changes the cell centre. This again changes the cellcentre-
61  // cellcentre edge across all faces of the cell.
62  // Note: this does not happen for unwarped faces but unfortunately
63  // we don't have this information.
64 
65  const labelList& faceOwner = mesh.faceOwner();
66  const labelList& faceNeighbour = mesh.faceNeighbour();
67  const cellList& cells = mesh.cells();
68  const label nInternalFaces = mesh.nInternalFaces();
69 
70  // Mark refined cells on old mesh
71  PackedBoolList oldRefineCell(map.nOldCells());
72 
73  forAll(oldCellsToRefine, i)
74  {
75  oldRefineCell.set(oldCellsToRefine[i], 1u);
76  }
77 
78  // Mark refined faces
79  PackedBoolList refinedInternalFace(nInternalFaces);
80 
81  // 1. Internal faces
82 
83  for (label faceI = 0; faceI < nInternalFaces; faceI++)
84  {
85  label oldOwn = map.cellMap()[faceOwner[faceI]];
86  label oldNei = map.cellMap()[faceNeighbour[faceI]];
87 
88  if
89  (
90  oldOwn >= 0
91  && oldRefineCell.get(oldOwn) == 0u
92  && oldNei >= 0
93  && oldRefineCell.get(oldNei) == 0u
94  )
95  {
96  // Unaffected face since both neighbours were not refined.
97  }
98  else
99  {
100  refinedInternalFace.set(faceI, 1u);
101  }
102  }
103 
104 
105  // 2. Boundary faces
106 
107  boolList refinedBoundaryFace(mesh.nFaces()-nInternalFaces, false);
108 
109  forAll(mesh.boundaryMesh(), patchI)
110  {
111  const polyPatch& pp = mesh.boundaryMesh()[patchI];
112 
113  label faceI = pp.start();
114 
115  forAll(pp, i)
116  {
117  label oldOwn = map.cellMap()[faceOwner[faceI]];
118 
119  if (oldOwn >= 0 && oldRefineCell.get(oldOwn) == 0u)
120  {
121  // owner did exist and wasn't refined.
122  }
123  else
124  {
125  refinedBoundaryFace[faceI-nInternalFaces] = true;
126  }
127  faceI++;
128  }
129  }
130 
131  // Synchronise refined face status
133  (
134  mesh,
135  refinedBoundaryFace,
136  orEqOp<bool>(),
137  false
138  );
139 
140 
141  // 3. Mark all faces affected by refinement. Refined faces are in
142  // - refinedInternalFace
143  // - refinedBoundaryFace
144  boolList changedFace(mesh.nFaces(), false);
145 
146  forAll(refinedInternalFace, faceI)
147  {
148  if (refinedInternalFace.get(faceI) == 1u)
149  {
150  const cell& ownFaces = cells[faceOwner[faceI]];
151  forAll(ownFaces, ownI)
152  {
153  changedFace[ownFaces[ownI]] = true;
154  }
155  const cell& neiFaces = cells[faceNeighbour[faceI]];
156  forAll(neiFaces, neiI)
157  {
158  changedFace[neiFaces[neiI]] = true;
159  }
160  }
161  }
162 
163  forAll(refinedBoundaryFace, i)
164  {
165  if (refinedBoundaryFace[i])
166  {
167  const cell& ownFaces = cells[faceOwner[i+nInternalFaces]];
168  forAll(ownFaces, ownI)
169  {
170  changedFace[ownFaces[ownI]] = true;
171  }
172  }
173  }
174 
176  (
177  mesh,
178  changedFace,
179  orEqOp<bool>(),
180  false
181  );
182 
183 
184  // Now we have in changedFace marked all affected faces. Pack.
185  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
186 
187  label nChanged = 0;
188 
189  forAll(changedFace, faceI)
190  {
191  if (changedFace[faceI])
192  {
193  nChanged++;
194  }
195  }
196 
197  changedFaces.setSize(nChanged);
198  nChanged = 0;
199 
200  forAll(changedFace, faceI)
201  {
202  if (changedFace[faceI])
203  {
204  changedFaces[nChanged++] = faceI;
205  }
206  }
207  }
208 
209  label nChangedFaces = changedFaces.size();
210  reduce(nChangedFaces, sumOp<label>());
211 
212  if (debug)
213  {
214  Pout<< "getChangedFaces : Detected "
215  << " local:" << changedFaces.size()
216  << " global:" << nChangedFaces
217  << " changed faces out of " << mesh.globalData().nTotalFaces()
218  << endl;
219 
220  faceSet changedFacesSet(mesh, "changedFaces", changedFaces);
221  Pout<< "getChangedFaces : Writing " << changedFaces.size()
222  << " changed faces to faceSet " << changedFacesSet.name()
223  << endl;
224  changedFacesSet.write();
225  }
226 
227  return changedFaces;
228 }
229 
230 
231 // Mark cell for refinement (if not already marked). Return false if
232 // refinelimit hit. Keeps running count (in nRefine) of cells marked for
233 // refinement
234 bool Foam::meshRefinement::markForRefine
235 (
236  const label markValue,
237  const label nAllowRefine,
238 
239  label& cellValue,
240  label& nRefine
241 )
242 {
243  if (cellValue == -1)
244  {
245  cellValue = markValue;
246  nRefine++;
247  }
248 
249  return nRefine <= nAllowRefine;
250 }
251 
252 
253 // Calculates list of cells to refine based on intersection with feature edge.
254 Foam::label Foam::meshRefinement::markFeatureRefinement
255 (
256  const point& keepPoint,
257  const PtrList<featureEdgeMesh>& featureMeshes,
258  const labelList& featureLevels,
259  const label nAllowRefine,
260 
261  labelList& refineCell,
262  label& nRefine
263 ) const
264 {
265  // We want to refine all cells containing a feature edge.
266  // - don't want to search over whole mesh
267  // - don't want to build octree for whole mesh
268  // - so use tracking to follow the feature edges
269  //
270  // 1. find non-manifold points on feature edges (i.e. start of feature edge
271  // or 'knot')
272  // 2. seed particle starting at keepPoint going to this non-manifold point
273  // 3. track particles to their non-manifold point
274  //
275  // 4. track particles across their connected feature edges, marking all
276  // visited cells with their level (through trackingData)
277  // 5. do 4 until all edges have been visited.
278 
279 
280  // Find all start cells of features. Is done by tracking from keepPoint.
281  Cloud<trackedParticle> cloud(mesh_, IDLList<trackedParticle>());
282 
283  // Create particles on whichever processor holds the keepPoint.
284  label cellI = mesh_.findCell(keepPoint);
285 
286  if (cellI != -1)
287  {
288  forAll(featureMeshes, featI)
289  {
290  const featureEdgeMesh& featureMesh = featureMeshes[featI];
291  const labelListList& pointEdges = featureMesh.pointEdges();
292 
293  forAll(pointEdges, pointI)
294  {
295  if (pointEdges[pointI].size() != 2)
296  {
297  if (debug)
298  {
299  Pout<< "Adding particle from point:" << pointI
300  << " coord:" << featureMesh.points()[pointI]
301  << " pEdges:" << pointEdges[pointI]
302  << endl;
303  }
304 
305  // Non-manifold point. Create particle.
306  cloud.addParticle
307  (
308  new trackedParticle
309  (
310  cloud,
311  keepPoint,
312  cellI,
313  featureMesh.points()[pointI], // endpos
314  featureLevels[featI], // level
315  featI, // featureMesh
316  pointI // end point
317  )
318  );
319  }
320  }
321  }
322  }
323 
324 
325  // Largest refinement level of any feature passed through
326  labelList maxFeatureLevel(mesh_.nCells(), -1);
327 
328  // Database to pass into trackedParticle::move
329  trackedParticle::trackData td(cloud, maxFeatureLevel);
330 
331  // Track all particles to their end position (= starting feature point)
332  cloud.move(td);
333 
334  // Reset level
335  maxFeatureLevel = -1;
336 
337  // Whether edge has been visited.
338  List<PackedBoolList> featureEdgeVisited(featureMeshes.size());
339 
340  forAll(featureMeshes, featI)
341  {
342  featureEdgeVisited[featI].setSize(featureMeshes[featI].edges().size());
343  featureEdgeVisited[featI] = 0u;
344  }
345 
346  while (true)
347  {
348  label nParticles = 0;
349 
350  // Make particle follow edge.
351  forAllIter(Cloud<trackedParticle>, cloud, iter)
352  {
353  trackedParticle& tp = iter();
354 
355  label featI = tp.i();
356  label pointI = tp.j();
357 
358  const featureEdgeMesh& featureMesh = featureMeshes[featI];
359  const labelList& pEdges = featureMesh.pointEdges()[pointI];
360 
361  // Particle now at pointI. Check connected edges to see which one
362  // we have to visit now.
363 
364  bool keepParticle = false;
365 
366  forAll(pEdges, i)
367  {
368  label edgeI = pEdges[i];
369 
370  if (featureEdgeVisited[featI].set(edgeI, 1u))
371  {
372  // Unvisited edge. Make the particle go to the other point
373  // on the edge.
374 
375  const edge& e = featureMesh.edges()[edgeI];
376  label otherPointI = e.otherVertex(pointI);
377 
378  tp.end() = featureMesh.points()[otherPointI];
379  tp.j() = otherPointI;
380  keepParticle = true;
381  break;
382  }
383  }
384 
385  if (!keepParticle)
386  {
387  // Particle at 'knot' where another particle already has been
388  // seeded. Delete particle.
389  cloud.deleteParticle(tp);
390  }
391  else
392  {
393  // Keep particle
394  nParticles++;
395  }
396  }
397 
398  reduce(nParticles, sumOp<label>());
399  if (nParticles == 0)
400  {
401  break;
402  }
403 
404  // Track all particles to their end position.
405  cloud.move(td);
406  }
407 
408 
409  // See which cells to refine. maxFeatureLevel will hold highest level
410  // of any feature edge that passed through.
411 
412  const labelList& cellLevel = meshCutter_.cellLevel();
413 
414  label oldNRefine = nRefine;
415 
416  forAll(maxFeatureLevel, cellI)
417  {
418  if (maxFeatureLevel[cellI] > cellLevel[cellI])
419  {
420  // Mark
421  if
422  (
423  !markForRefine
424  (
425  0, // surface (n/a)
426  nAllowRefine,
427  refineCell[cellI],
428  nRefine
429  )
430  )
431  {
432  // Reached limit
433  break;
434  }
435  }
436  }
437 
438  if
439  (
440  returnReduce(nRefine, sumOp<label>())
441  > returnReduce(nAllowRefine, sumOp<label>())
442  )
443  {
444  Info<< "Reached refinement limit." << endl;
445  }
446 
447  return returnReduce(nRefine-oldNRefine, sumOp<label>());
448 }
449 
450 
451 // Mark cells for non-surface intersection based refinement.
452 Foam::label Foam::meshRefinement::markInternalRefinement
453 (
454  const label nAllowRefine,
455 
456  labelList& refineCell,
457  label& nRefine
458 ) const
459 {
460  const labelList& cellLevel = meshCutter_.cellLevel();
461  const pointField& cellCentres = mesh_.cellCentres();
462 
463  label oldNRefine = nRefine;
464 
465  // Collect cells to test
466  pointField testCc(cellLevel.size()-nRefine);
467  labelList testLevels(cellLevel.size()-nRefine);
468  label testI = 0;
469 
470  forAll(cellLevel, cellI)
471  {
472  if (refineCell[cellI] == -1)
473  {
474  testCc[testI] = cellCentres[cellI];
475  testLevels[testI] = cellLevel[cellI];
476  testI++;
477  }
478  }
479 
480  // Do test to see whether cells is inside/outside shell with higher level
481  labelList maxLevel;
482  shells_.findHigherLevel(testCc, testLevels, maxLevel);
483 
484  // Mark for refinement. Note that we didn't store the original cellID so
485  // now just reloop in same order.
486  testI = 0;
487  forAll(cellLevel, cellI)
488  {
489  if (refineCell[cellI] == -1)
490  {
491  if (maxLevel[testI] > testLevels[testI])
492  {
493  bool reachedLimit = !markForRefine
494  (
495  maxLevel[testI], // mark with any positive value
496  nAllowRefine,
497  refineCell[cellI],
498  nRefine
499  );
500 
501  if (reachedLimit)
502  {
503  if (debug)
504  {
505  Pout<< "Stopped refining internal cells"
506  << " since reaching my cell limit of "
507  << mesh_.nCells()+7*nRefine << endl;
508  }
509  break;
510  }
511  }
512  testI++;
513  }
514  }
515 
516  if
517  (
518  returnReduce(nRefine, sumOp<label>())
519  > returnReduce(nAllowRefine, sumOp<label>())
520  )
521  {
522  Info<< "Reached refinement limit." << endl;
523  }
524 
525  return returnReduce(nRefine-oldNRefine, sumOp<label>());
526 }
527 
528 
529 // Collect faces that are intersected and whose neighbours aren't yet marked
530 // for refinement.
531 Foam::labelList Foam::meshRefinement::getRefineCandidateFaces
532 (
533  const labelList& refineCell
534 ) const
535 {
536  labelList testFaces(mesh_.nFaces());
537 
538  label nTest = 0;
539 
540  forAll(surfaceIndex_, faceI)
541  {
542  if (surfaceIndex_[faceI] != -1)
543  {
544  label own = mesh_.faceOwner()[faceI];
545 
546  if (mesh_.isInternalFace(faceI))
547  {
548  label nei = mesh_.faceNeighbour()[faceI];
549 
550  if (refineCell[own] == -1 || refineCell[nei] == -1)
551  {
552  testFaces[nTest++] = faceI;
553  }
554  }
555  else
556  {
557  if (refineCell[own] == -1)
558  {
559  testFaces[nTest++] = faceI;
560  }
561  }
562  }
563  }
564  testFaces.setSize(nTest);
565 
566  return testFaces;
567 }
568 
569 
570 // Mark cells for surface intersection based refinement.
571 Foam::label Foam::meshRefinement::markSurfaceRefinement
572 (
573  const label nAllowRefine,
574  const labelList& neiLevel,
575  const pointField& neiCc,
576 
577  labelList& refineCell,
578  label& nRefine
579 ) const
580 {
581  const labelList& cellLevel = meshCutter_.cellLevel();
582  const pointField& cellCentres = mesh_.cellCentres();
583 
584  label oldNRefine = nRefine;
585 
586  // Use cached surfaceIndex_ to detect if any intersection. If so
587  // re-intersect to determine level wanted.
588 
589  // Collect candidate faces
590  // ~~~~~~~~~~~~~~~~~~~~~~~
591 
592  labelList testFaces(getRefineCandidateFaces(refineCell));
593 
594  // Collect segments
595  // ~~~~~~~~~~~~~~~~
596 
597  pointField start(testFaces.size());
598  pointField end(testFaces.size());
599  labelList minLevel(testFaces.size());
600 
601  forAll(testFaces, i)
602  {
603  label faceI = testFaces[i];
604 
605  label own = mesh_.faceOwner()[faceI];
606 
607  if (mesh_.isInternalFace(faceI))
608  {
609  label nei = mesh_.faceNeighbour()[faceI];
610 
611  start[i] = cellCentres[own];
612  end[i] = cellCentres[nei];
613  minLevel[i] = min(cellLevel[own], cellLevel[nei]);
614  }
615  else
616  {
617  label bFaceI = faceI - mesh_.nInternalFaces();
618 
619  start[i] = cellCentres[own];
620  end[i] = neiCc[bFaceI];
621  minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
622  }
623  }
624 
625 
626  // Do test for higher intersections
627  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
628 
629  labelList surfaceHit;
630  labelList surfaceMinLevel;
631  surfaces_.findHigherIntersection
632  (
633  start,
634  end,
635  minLevel,
636 
637  surfaceHit,
638  surfaceMinLevel
639  );
640 
641 
642  // Mark cells for refinement
643  // ~~~~~~~~~~~~~~~~~~~~~~~~~
644 
645  forAll(testFaces, i)
646  {
647  label faceI = testFaces[i];
648 
649  label surfI = surfaceHit[i];
650 
651  if (surfI != -1)
652  {
653  // Found intersection with surface with higher wanted
654  // refinement. Check if the level field on that surface
655  // specifies an even higher level. Note:this is weird. Should
656  // do the check with the surfaceMinLevel whilst intersecting the
657  // surfaces?
658 
659  label own = mesh_.faceOwner()[faceI];
660 
661  if (surfaceMinLevel[i] > cellLevel[own])
662  {
663  // Owner needs refining
664  if
665  (
666  !markForRefine
667  (
668  surfI,
669  nAllowRefine,
670  refineCell[own],
671  nRefine
672  )
673  )
674  {
675  break;
676  }
677  }
678 
679  if (mesh_.isInternalFace(faceI))
680  {
681  label nei = mesh_.faceNeighbour()[faceI];
682  if (surfaceMinLevel[i] > cellLevel[nei])
683  {
684  // Neighbour needs refining
685  if
686  (
687  !markForRefine
688  (
689  surfI,
690  nAllowRefine,
691  refineCell[nei],
692  nRefine
693  )
694  )
695  {
696  break;
697  }
698  }
699  }
700  }
701  }
702 
703  if
704  (
705  returnReduce(nRefine, sumOp<label>())
706  > returnReduce(nAllowRefine, sumOp<label>())
707  )
708  {
709  Info<< "Reached refinement limit." << endl;
710  }
711 
712  return returnReduce(nRefine-oldNRefine, sumOp<label>());
713 }
714 
715 
716 // Checks if multiple intersections of a cell (by a surface with a higher
717 // max than the cell level) and if so if the normals at these intersections
718 // make a large angle.
719 // Returns false if the nRefine limit has been reached, true otherwise.
720 bool Foam::meshRefinement::checkCurvature
721 (
722  const scalar curvature,
723  const label nAllowRefine,
724 
725  const label surfaceLevel, // current intersection max level
726  const vector& surfaceNormal,// current intersection normal
727 
728  const label cellI,
729 
730  label& cellMaxLevel, // cached max surface level for this cell
731  vector& cellMaxNormal, // cached surface normal for this cell
732 
733  labelList& refineCell,
734  label& nRefine
735 ) const
736 {
737  const labelList& cellLevel = meshCutter_.cellLevel();
738 
739  // Test if surface applicable
740  if (surfaceLevel > cellLevel[cellI])
741  {
742  if (cellMaxLevel == -1)
743  {
744  // First visit of cell. Store
745  cellMaxLevel = surfaceLevel;
746  cellMaxNormal = surfaceNormal;
747  }
748  else
749  {
750  // Second or more visit. Check curvature.
751  if ((cellMaxNormal & surfaceNormal) < curvature)
752  {
753  return markForRefine
754  (
755  surfaceLevel, // mark with any non-neg number.
756  nAllowRefine,
757  refineCell[cellI],
758  nRefine
759  );
760  }
761 
762  // Set normal to that of highest surface. Not really necessary
763  // over here but we reuse cellMax info when doing coupled faces.
764  if (surfaceLevel > cellMaxLevel)
765  {
766  cellMaxLevel = surfaceLevel;
767  cellMaxNormal = surfaceNormal;
768  }
769  }
770  }
771 
772  // Did not reach refinement limit.
773  return true;
774 }
775 
776 
777 // Mark cells for surface curvature based refinement.
778 Foam::label Foam::meshRefinement::markSurfaceCurvatureRefinement
779 (
780  const scalar curvature,
781  const label nAllowRefine,
782  const labelList& neiLevel,
783  const pointField& neiCc,
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  // 1. local test: any cell on more than one surface gets refined
795  // (if its current level is < max of the surface max level)
796 
797  // 2. 'global' test: any cell on only one surface with a neighbour
798  // on a different surface gets refined (if its current level etc.)
799 
800 
801  // Collect candidate faces (i.e. intersecting any surface and
802  // owner/neighbour not yet refined.
803  labelList testFaces(getRefineCandidateFaces(refineCell));
804 
805  // Collect segments
806  pointField start(testFaces.size());
807  pointField end(testFaces.size());
808  labelList minLevel(testFaces.size());
809 
810  forAll(testFaces, i)
811  {
812  label faceI = testFaces[i];
813 
814  label own = mesh_.faceOwner()[faceI];
815 
816  if (mesh_.isInternalFace(faceI))
817  {
818  label nei = mesh_.faceNeighbour()[faceI];
819 
820  start[i] = cellCentres[own];
821  end[i] = cellCentres[nei];
822  minLevel[i] = min(cellLevel[own], cellLevel[nei]);
823  }
824  else
825  {
826  label bFaceI = faceI - mesh_.nInternalFaces();
827 
828  start[i] = cellCentres[own];
829  end[i] = neiCc[bFaceI];
830  minLevel[i] = min(cellLevel[own], neiLevel[bFaceI]);
831  }
832  }
833 
834  // Test for all intersections (with surfaces of higher max level than
835  // minLevel) and cache per cell the max surface level and the local normal
836  // on that surface.
837  labelList cellMaxLevel(mesh_.nCells(), -1);
838  vectorField cellMaxNormal(mesh_.nCells(), vector::zero);
839 
840  {
841  // Per segment the normals of the surfaces
842  List<vectorList> surfaceNormal;
843  // Per segment the list of levels of the surfaces
844  labelListList surfaceLevel;
845 
846  surfaces_.findAllHigherIntersections
847  (
848  start,
849  end,
850  minLevel, // max level of surface has to be bigger
851  // than min level of neighbouring cells
852  surfaceNormal,
853  surfaceLevel
854  );
855  // Clear out unnecessary data
856  start.clear();
857  end.clear();
858  minLevel.clear();
859 
860  // Extract per cell information on the surface with the highest max
861  forAll(testFaces, i)
862  {
863  label faceI = testFaces[i];
864  label own = mesh_.faceOwner()[faceI];
865 
866  const vectorList& fNormals = surfaceNormal[i];
867  const labelList& fLevels = surfaceLevel[i];
868 
869  forAll(fLevels, hitI)
870  {
871  checkCurvature
872  (
873  curvature,
874  nAllowRefine,
875 
876  fLevels[hitI],
877  fNormals[hitI],
878 
879  own,
880  cellMaxLevel[own],
881  cellMaxNormal[own],
882 
883  refineCell,
884  nRefine
885  );
886  }
887 
888  if (mesh_.isInternalFace(faceI))
889  {
890  label nei = mesh_.faceNeighbour()[faceI];
891 
892  forAll(fLevels, hitI)
893  {
894  checkCurvature
895  (
896  curvature,
897  nAllowRefine,
898 
899  fLevels[hitI],
900  fNormals[hitI],
901 
902  nei,
903  cellMaxLevel[nei],
904  cellMaxNormal[nei],
905 
906  refineCell,
907  nRefine
908  );
909  }
910  }
911  }
912  }
913 
914  // 2. Find out a measure of surface curvature and region edges.
915  // Send over surface region and surface normal to neighbour cell.
916 
917  labelList neiBndMaxLevel(mesh_.nFaces()-mesh_.nInternalFaces());
918  vectorField neiBndMaxNormal(mesh_.nFaces()-mesh_.nInternalFaces());
919 
920  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
921  {
922  label bFaceI = faceI-mesh_.nInternalFaces();
923  label own = mesh_.faceOwner()[faceI];
924 
925  neiBndMaxLevel[bFaceI] = cellMaxLevel[own];
926  neiBndMaxNormal[bFaceI] = cellMaxNormal[own];
927  }
928  syncTools::swapBoundaryFaceList(mesh_, neiBndMaxLevel, false);
929  syncTools::swapBoundaryFaceList(mesh_, neiBndMaxNormal, false);
930 
931  // Loop over all faces. Could only be checkFaces.. except if they're coupled
932 
933  // Internal faces
934  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
935  {
936  label own = mesh_.faceOwner()[faceI];
937  label nei = mesh_.faceNeighbour()[faceI];
938 
939  if (cellMaxLevel[own] != -1 && cellMaxLevel[nei] != -1)
940  {
941  // Have valid data on both sides. Check curvature.
942  if ((cellMaxNormal[own] & cellMaxNormal[nei]) < curvature)
943  {
944  // See which side to refine
945  if (cellLevel[own] < cellMaxLevel[own])
946  {
947  if
948  (
949  !markForRefine
950  (
951  cellMaxLevel[own],
952  nAllowRefine,
953  refineCell[own],
954  nRefine
955  )
956  )
957  {
958  if (debug)
959  {
960  Pout<< "Stopped refining since reaching my cell"
961  << " limit of " << mesh_.nCells()+7*nRefine
962  << endl;
963  }
964  break;
965  }
966  }
967 
968  if (cellLevel[nei] < cellMaxLevel[nei])
969  {
970  if
971  (
972  !markForRefine
973  (
974  cellMaxLevel[nei],
975  nAllowRefine,
976  refineCell[nei],
977  nRefine
978  )
979  )
980  {
981  if (debug)
982  {
983  Pout<< "Stopped refining since reaching my cell"
984  << " limit of " << mesh_.nCells()+7*nRefine
985  << endl;
986  }
987  break;
988  }
989  }
990  }
991  }
992  }
993  // Boundary faces
994  for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
995  {
996  label own = mesh_.faceOwner()[faceI];
997  label bFaceI = faceI - mesh_.nInternalFaces();
998 
999  if (cellLevel[own] < cellMaxLevel[own] && neiBndMaxLevel[bFaceI] != -1)
1000  {
1001  // Have valid data on both sides. Check curvature.
1002  if ((cellMaxNormal[own] & neiBndMaxNormal[bFaceI]) < curvature)
1003  {
1004  if
1005  (
1006  !markForRefine
1007  (
1008  cellMaxLevel[own],
1009  nAllowRefine,
1010  refineCell[own],
1011  nRefine
1012  )
1013  )
1014  {
1015  if (debug)
1016  {
1017  Pout<< "Stopped refining since reaching my cell"
1018  << " limit of " << mesh_.nCells()+7*nRefine
1019  << endl;
1020  }
1021  break;
1022  }
1023  }
1024  }
1025  }
1026 
1027  if
1028  (
1029  returnReduce(nRefine, sumOp<label>())
1030  > returnReduce(nAllowRefine, sumOp<label>())
1031  )
1032  {
1033  Info<< "Reached refinement limit." << endl;
1034  }
1035 
1036  return returnReduce(nRefine-oldNRefine, sumOp<label>());
1037 }
1038 
1039 
1040 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1041 
1042 // Calculate list of cells to refine. Gets for any edge (start - end)
1043 // whether it intersects the surface. Does more accurate test and checks
1044 // the wanted level on the surface intersected.
1045 // Does approximate precalculation of how many cells can be refined before
1046 // hitting overall limit maxGlobalCells.
1049  const point& keepPoint,
1050  const scalar curvature,
1051 
1052  const PtrList<featureEdgeMesh>& featureMeshes,
1053  const labelList& featureLevels,
1054 
1055  const bool featureRefinement,
1056  const bool internalRefinement,
1057  const bool surfaceRefinement,
1058  const bool curvatureRefinement,
1059  const label maxGlobalCells,
1060  const label maxLocalCells
1061 ) const
1062 {
1063  label totNCells = mesh_.globalData().nTotalCells();
1064 
1065  labelList cellsToRefine;
1066 
1067  if (totNCells >= maxGlobalCells)
1068  {
1069  Info<< "No cells marked for refinement since reached limit "
1070  << maxGlobalCells << '.' << endl;
1071  }
1072  else
1073  {
1074  // Every cell I refine adds 7 cells. Estimate the number of cells
1075  // I am allowed to refine.
1076  // Assume perfect distribution so can only refine as much the fraction
1077  // of the mesh I hold. This prediction step prevents us having to do
1078  // lots of reduces to keep count of the total number of cells selected
1079  // for refinement.
1080 
1081  //scalar fraction = scalar(mesh_.nCells())/totNCells;
1082  //label myMaxCells = label(maxGlobalCells*fraction);
1083  //label nAllowRefine = (myMaxCells - mesh_.nCells())/7;
1085  //
1086  //Pout<< "refineCandidates:" << nl
1087  // << " total cells:" << totNCells << nl
1088  // << " local cells:" << mesh_.nCells() << nl
1089  // << " local fraction:" << fraction << nl
1090  // << " local allowable cells:" << myMaxCells << nl
1091  // << " local allowable refinement:" << nAllowRefine << nl
1092  // << endl;
1093 
1094  //- Disable refinement shortcut. nAllowRefine is per processor limit.
1095  label nAllowRefine = labelMax / Pstream::nProcs();
1096 
1097  // Marked for refinement (>= 0) or not (-1). Actual value is the
1098  // index of the surface it intersects.
1099  labelList refineCell(mesh_.nCells(), -1);
1100  label nRefine = 0;
1101 
1102 
1103  // Swap neighbouring cell centres and cell level
1104  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
1105  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
1106  calcNeighbourData(neiLevel, neiCc);
1107 
1108 
1109 
1110  // Cells pierced by feature lines
1111  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1112 
1113  if (featureRefinement)
1114  {
1115  label nFeatures = markFeatureRefinement
1116  (
1117  keepPoint,
1118  featureMeshes,
1119  featureLevels,
1120  nAllowRefine,
1121 
1122  refineCell,
1123  nRefine
1124  );
1125 
1126  Info<< "Marked for refinement due to explicit features : "
1127  << nFeatures << " cells." << endl;
1128  }
1129 
1130  // Inside refinement shells
1131  // ~~~~~~~~~~~~~~~~~~~~~~~~
1132 
1133  if (internalRefinement)
1134  {
1135  label nShell = markInternalRefinement
1136  (
1137  nAllowRefine,
1138 
1139  refineCell,
1140  nRefine
1141  );
1142  Info<< "Marked for refinement due to refinement shells : "
1143  << nShell << " cells." << endl;
1144  }
1145 
1146  // Refinement based on intersection of surface
1147  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1148 
1149  if (surfaceRefinement)
1150  {
1151  label nSurf = markSurfaceRefinement
1152  (
1153  nAllowRefine,
1154  neiLevel,
1155  neiCc,
1156 
1157  refineCell,
1158  nRefine
1159  );
1160  Info<< "Marked for refinement due to surface intersection : "
1161  << nSurf << " cells." << endl;
1162  }
1163 
1164  // Refinement based on curvature of surface
1165  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1166 
1167  if
1168  (
1169  curvatureRefinement
1170  && (curvature >= -1 && curvature <= 1)
1171  && (surfaces_.minLevel() != surfaces_.maxLevel())
1172  )
1173  {
1174  label nCurv = markSurfaceCurvatureRefinement
1175  (
1176  curvature,
1177  nAllowRefine,
1178  neiLevel,
1179  neiCc,
1180 
1181  refineCell,
1182  nRefine
1183  );
1184  Info<< "Marked for refinement due to curvature/regions : "
1185  << nCurv << " cells." << endl;
1186  }
1187 
1188  // Pack cells-to-refine
1189  // ~~~~~~~~~~~~~~~~~~~~
1190 
1191  cellsToRefine.setSize(nRefine);
1192  nRefine = 0;
1193 
1194  forAll(refineCell, cellI)
1195  {
1196  if (refineCell[cellI] != -1)
1197  {
1198  cellsToRefine[nRefine++] = cellI;
1199  }
1200  }
1201  }
1202 
1203  return cellsToRefine;
1204 }
1205 
1206 
1209  const labelList& cellsToRefine
1210 )
1211 {
1212  // Mesh changing engine.
1213  polyTopoChange meshMod(mesh_);
1214 
1215  // Play refinement commands into mesh changer.
1216  meshCutter_.setRefinement(cellsToRefine, meshMod);
1217 
1218  // Create mesh (no inflation), return map from old to new mesh.
1219  autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh_, false);
1220 
1221  // Update fields
1222  mesh_.updateMesh(map);
1223 
1224  // Optionally inflate mesh
1225  if (map().hasMotionPoints())
1226  {
1227  mesh_.movePoints(map().preMotionPoints());
1228  }
1229  else
1230  {
1231  // Delete mesh volumes.
1232  mesh_.clearOut();
1233  }
1234 
1235  if (overwrite())
1236  {
1237  mesh_.setInstance(oldInstance());
1238  }
1239 
1240  // Update intersection info
1241  updateMesh(map, getChangedFaces(map, cellsToRefine));
1242 
1243  return map;
1244 }
1245 
1246 
1247 // Do refinement of consistent set of cells followed by truncation and
1248 // load balancing.
1252  const string& msg,
1253  decompositionMethod& decomposer,
1254  fvMeshDistribute& distributor,
1255  const labelList& cellsToRefine,
1256  const scalar maxLoadUnbalance
1257 )
1258 {
1259  // Do all refinement
1260  refine(cellsToRefine);
1261 
1262  if (debug)
1263  {
1264  Pout<< "Writing refined but unbalanced " << msg
1265  << " mesh to time " << timeName() << endl;
1266  write
1267  (
1268  debug,
1269  mesh_.time().path()
1270  /timeName()
1271  );
1272  Pout<< "Dumped debug data in = "
1273  << mesh_.time().cpuTimeIncrement() << " s" << endl;
1274 
1275  // test all is still synced across proc patches
1276  checkData();
1277  }
1278 
1279  Info<< "Refined mesh in = "
1280  << mesh_.time().cpuTimeIncrement() << " s" << endl;
1281  printMeshInfo(debug, "After refinement " + msg);
1282 
1283 
1284  // Load balancing
1285  // ~~~~~~~~~~~~~~
1286 
1288 
1289  if (Pstream::nProcs() > 1)
1290  {
1291  scalar nIdealCells =
1292  mesh_.globalData().nTotalCells()
1293  / Pstream::nProcs();
1294 
1295  scalar unbalance = returnReduce
1296  (
1297  mag(1.0-mesh_.nCells()/nIdealCells),
1298  maxOp<scalar>()
1299  );
1300 
1301  if (unbalance <= maxLoadUnbalance)
1302  {
1303  Info<< "Skipping balancing since max unbalance " << unbalance
1304  << " is less than allowable " << maxLoadUnbalance
1305  << endl;
1306  }
1307  else
1308  {
1309  scalarField cellWeights(mesh_.nCells(), 1);
1310 
1311  distMap = balance
1312  (
1313  false, //keepZoneFaces
1314  false, //keepBaffles
1315  cellWeights,
1316  decomposer,
1317  distributor
1318  );
1319 
1320  Info<< "Balanced mesh in = "
1321  << mesh_.time().cpuTimeIncrement() << " s" << endl;
1322 
1323  printMeshInfo(debug, "After balancing " + msg);
1324 
1325 
1326  if (debug)
1327  {
1328  Pout<< "Writing balanced " << msg
1329  << " mesh to time " << timeName() << endl;
1330  write
1331  (
1332  debug,
1333  mesh_.time().path()/timeName()
1334  );
1335  Pout<< "Dumped debug data in = "
1336  << mesh_.time().cpuTimeIncrement() << " s" << endl;
1337 
1338  // test all is still synced across proc patches
1339  checkData();
1340  }
1341  }
1342  }
1343 
1344  return distMap;
1345 }
1346 
1347 
1348 // Do load balancing followed by refinement of consistent set of cells.
1352  const string& msg,
1353  decompositionMethod& decomposer,
1354  fvMeshDistribute& distributor,
1355  const labelList& initCellsToRefine,
1356  const scalar maxLoadUnbalance
1357 )
1358 {
1359  labelList cellsToRefine(initCellsToRefine);
1360 
1361  //{
1362  // globalIndex globalCells(mesh_.nCells());
1363  //
1364  // Info<< "** Distribution before balancing/refining:" << endl;
1365  // for (label procI = 0; procI < Pstream::nProcs(); procI++)
1366  // {
1367  // Info<< " " << procI << '\t'
1368  // << globalCells.localSize(procI) << endl;
1369  // }
1370  // Info<< endl;
1371  //}
1372  //{
1373  // globalIndex globalCells(cellsToRefine.size());
1374  //
1375  // Info<< "** Cells to be refined:" << endl;
1376  // for (label procI = 0; procI < Pstream::nProcs(); procI++)
1377  // {
1378  // Info<< " " << procI << '\t'
1379  // << globalCells.localSize(procI) << endl;
1380  // }
1381  // Info<< endl;
1382  //}
1383 
1384 
1385  // Load balancing
1386  // ~~~~~~~~~~~~~~
1387 
1389 
1390  if (Pstream::nProcs() > 1)
1391  {
1392  // First check if we need to balance at all. Precalculate number of
1393  // cells after refinement and see what maximum difference is.
1394  scalar nNewCells = scalar(mesh_.nCells() + 7*cellsToRefine.size());
1395  scalar nIdealNewCells =
1396  returnReduce(nNewCells, sumOp<scalar>())
1397  / Pstream::nProcs();
1398  scalar unbalance = returnReduce
1399  (
1400  mag(1.0-nNewCells/nIdealNewCells),
1401  maxOp<scalar>()
1402  );
1403 
1404  if (unbalance <= maxLoadUnbalance)
1405  {
1406  Info<< "Skipping balancing since max unbalance " << unbalance
1407  << " is less than allowable " << maxLoadUnbalance
1408  << endl;
1409  }
1410  else
1411  {
1412  scalarField cellWeights(mesh_.nCells(), 1);
1413  forAll(cellsToRefine, i)
1414  {
1415  cellWeights[cellsToRefine[i]] += 7;
1416  }
1417 
1418  distMap = balance
1419  (
1420  false, //keepZoneFaces
1421  false, //keepBaffles
1422  cellWeights,
1423  decomposer,
1424  distributor
1425  );
1426 
1427  // Update cells to refine
1428  distMap().distributeCellIndices(cellsToRefine);
1429 
1430  Info<< "Balanced mesh in = "
1431  << mesh_.time().cpuTimeIncrement() << " s" << endl;
1432  }
1433 
1434  //{
1435  // globalIndex globalCells(mesh_.nCells());
1436  //
1437  // Info<< "** Distribution after balancing:" << endl;
1438  // for (label procI = 0; procI < Pstream::nProcs(); procI++)
1439  // {
1440  // Info<< " " << procI << '\t'
1441  // << globalCells.localSize(procI) << endl;
1442  // }
1443  // Info<< endl;
1444  //}
1445 
1446  printMeshInfo(debug, "After balancing " + msg);
1447 
1448  if (debug)
1449  {
1450  Pout<< "Writing balanced " << msg
1451  << " mesh to time " << timeName() << endl;
1452  write
1453  (
1454  debug,
1455  mesh_.time().path()/timeName()
1456  );
1457  Pout<< "Dumped debug data in = "
1458  << mesh_.time().cpuTimeIncrement() << " s" << endl;
1459 
1460  // test all is still synced across proc patches
1461  checkData();
1462  }
1463  }
1464 
1465 
1466  // Refinement
1467  // ~~~~~~~~~~
1468 
1469  refine(cellsToRefine);
1470 
1471  if (debug)
1472  {
1473  Pout<< "Writing refined " << msg
1474  << " mesh to time " << timeName() << endl;
1475  write
1476  (
1477  debug,
1478  mesh_.time().path()
1479  /timeName()
1480  );
1481  Pout<< "Dumped debug data in = "
1482  << mesh_.time().cpuTimeIncrement() << " s" << endl;
1483 
1484  // test all is still synced across proc patches
1485  checkData();
1486  }
1487 
1488  Info<< "Refined mesh in = "
1489  << mesh_.time().cpuTimeIncrement() << " s" << endl;
1490 
1491  //{
1492  // globalIndex globalCells(mesh_.nCells());
1493  //
1494  // Info<< "** After refinement distribution:" << endl;
1495  // for (label procI = 0; procI < Pstream::nProcs(); procI++)
1496  // {
1497  // Info<< " " << procI << '\t'
1498  // << globalCells.localSize(procI) << endl;
1499  // }
1500  // Info<< endl;
1501  //}
1502 
1503  printMeshInfo(debug, "After refinement " + msg);
1504 
1505  return distMap;
1506 }
1507 
1508 
1509 // ************************ vim: set sw=4 sts=4 et: ************************ //