meshRefinementProblemCells.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2016 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 "fvMesh.H"
28 #include "syncTools.H"
29 #include "Time.H"
30 #include "refinementSurfaces.H"
31 #include "pointSet.H"
32 #include "faceSet.H"
33 #include "indirectPrimitivePatch.H"
34 #include "cellSet.H"
35 #include "searchableSurfaces.H"
36 #include "polyMeshGeometry.H"
37 #include "IOmanip.H"
38 #include "unitConversion.H"
39 #include "snappySnapDriver.H"
40 
41 #include "snapParameters.H"
42 #include "motionSmoother.H"
43 #include "topoDistanceData.H"
44 #include "FaceCellWave.H"
45 
46 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
47 
48 void Foam::meshRefinement::markBoundaryFace
49 (
50  const label facei,
51  boolList& isBoundaryFace,
52  boolList& isBoundaryEdge,
53  boolList& isBoundaryPoint
54 ) const
55 {
56  isBoundaryFace[facei] = true;
57 
58  const labelList& fEdges = mesh_.faceEdges(facei);
59 
60  forAll(fEdges, fp)
61  {
62  isBoundaryEdge[fEdges[fp]] = true;
63  }
64 
65  const face& f = mesh_.faces()[facei];
66 
67  forAll(f, fp)
68  {
69  isBoundaryPoint[f[fp]] = true;
70  }
71 }
72 
73 
74 void Foam::meshRefinement::findNearest
75 (
76  const labelList& meshFaces,
77  List<pointIndexHit>& nearestInfo,
78  labelList& nearestSurface,
79  labelList& nearestRegion,
80  vectorField& nearestNormal
81 ) const
82 {
83  pointField fc(meshFaces.size());
84  forAll(meshFaces, i)
85  {
86  fc[i] = mesh_.faceCentres()[meshFaces[i]];
87  }
88 
89  const labelList allSurfaces(identity(surfaces_.surfaces().size()));
90 
91  surfaces_.findNearest
92  (
93  allSurfaces,
94  fc,
95  scalarField(fc.size(), sqr(GREAT)), // sqr of attraction
96  nearestSurface,
97  nearestInfo
98  );
99 
100  // Do normal testing per surface.
101  nearestNormal.setSize(nearestInfo.size());
102  nearestRegion.setSize(nearestInfo.size());
103 
104  forAll(allSurfaces, surfI)
105  {
106  DynamicList<pointIndexHit> localHits;
107 
108  forAll(nearestSurface, i)
109  {
110  if (nearestSurface[i] == surfI)
111  {
112  localHits.append(nearestInfo[i]);
113  }
114  }
115 
116  label geomI = surfaces_.surfaces()[surfI];
117 
118  pointField localNormals;
119  surfaces_.geometry()[geomI].getNormal(localHits, localNormals);
120 
121  labelList localRegion;
122  surfaces_.geometry()[geomI].getRegion(localHits, localRegion);
123 
124  label localI = 0;
125  forAll(nearestSurface, i)
126  {
127  if (nearestSurface[i] == surfI)
128  {
129  nearestNormal[i] = localNormals[localI];
130  nearestRegion[i] = localRegion[localI];
131  localI++;
132  }
133  }
134  }
135 }
136 
137 
138 Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
139 (
140  const scalarField& perpendicularAngle,
141  const labelList& globalToPatch
142 ) const
143 {
144  // Construct addressing engine from all patches added for meshing.
145  autoPtr<indirectPrimitivePatch> ppPtr
146  (
148  (
149  mesh_,
150  meshedPatches()
151  )
152  );
153  const indirectPrimitivePatch& pp = ppPtr();
154 
155 
156  // 1. Collect faces to test
157  // ~~~~~~~~~~~~~~~~~~~~~~~~
158 
159  DynamicList<label> candidateFaces(pp.size()/20);
160 
161  const labelListList& edgeFaces = pp.edgeFaces();
162 
163  const labelList& cellLevel = meshCutter_.cellLevel();
164 
165  forAll(edgeFaces, edgeI)
166  {
167  const labelList& eFaces = edgeFaces[edgeI];
168 
169  if (eFaces.size() == 2)
170  {
171  label face0 = pp.addressing()[eFaces[0]];
172  label face1 = pp.addressing()[eFaces[1]];
173 
174  label cell0 = mesh_.faceOwner()[face0];
175  label cell1 = mesh_.faceOwner()[face1];
176 
177  if (cellLevel[cell0] > cellLevel[cell1])
178  {
179  // cell0 smaller.
180  const vector& n0 = pp.faceNormals()[eFaces[0]];
181  const vector& n1 = pp.faceNormals()[eFaces[1]];
182 
183  if (mag(n0 & n1) < 0.1)
184  {
185  candidateFaces.append(face0);
186  }
187  }
188  else if (cellLevel[cell1] > cellLevel[cell0])
189  {
190  // cell1 smaller.
191  const vector& n0 = pp.faceNormals()[eFaces[0]];
192  const vector& n1 = pp.faceNormals()[eFaces[1]];
193 
194  if (mag(n0 & n1) < 0.1)
195  {
196  candidateFaces.append(face1);
197  }
198  }
199  }
200  }
201  candidateFaces.shrink();
202 
203  Info<< "Testing " << returnReduce(candidateFaces.size(), sumOp<label>())
204  << " faces on edge-connected cells of differing level."
205  << endl;
206 
207  if (debug&meshRefinement::MESH)
208  {
209  faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
210  fSet.instance() = timeName();
211  Pout<< "Writing " << fSet.size()
212  << " with problematic topology to faceSet "
213  << fSet.objectPath() << endl;
214  fSet.write();
215  }
216 
217 
218  // 2. Find nearest surface on candidate faces
219  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
220 
221  List<pointIndexHit> nearestInfo;
222  labelList nearestSurface;
223  labelList nearestRegion;
224  vectorField nearestNormal;
225  findNearest
226  (
227  candidateFaces,
228  nearestInfo,
229  nearestSurface,
230  nearestRegion,
231  nearestNormal
232  );
233 
234 
235  // 3. Test angle to surface
236  // ~~~~~~~~~~~~~~~~~~~~~~~~
237 
238  Map<label> candidateCells(candidateFaces.size());
239 
240  faceSet perpFaces(mesh_, "perpendicularFaces", pp.size()/100);
241 
242  forAll(candidateFaces, i)
243  {
244  label facei = candidateFaces[i];
245 
246  vector n = mesh_.faceAreas()[facei];
247  n /= mag(n);
248 
249  label region = surfaces_.globalRegion
250  (
251  nearestSurface[i],
252  nearestRegion[i]
253  );
254 
255  scalar angle = degToRad(perpendicularAngle[region]);
256 
257  if (angle >= 0)
258  {
259  if (mag(n & nearestNormal[i]) < Foam::sin(angle))
260  {
261  perpFaces.insert(facei);
262  candidateCells.insert
263  (
264  mesh_.faceOwner()[facei],
265  globalToPatch[region]
266  );
267  }
268  }
269  }
270 
271  if (debug&meshRefinement::MESH)
272  {
273  perpFaces.instance() = timeName();
274  Pout<< "Writing " << perpFaces.size()
275  << " faces that are perpendicular to the surface to set "
276  << perpFaces.objectPath() << endl;
277  perpFaces.write();
278  }
279  return candidateCells;
280 }
281 
282 
283 // Check if moving face to new points causes a 'collapsed' face.
284 // Uses new point position only for the face, not the neighbouring
285 // cell centres
286 bool Foam::meshRefinement::isCollapsedFace
287 (
288  const pointField& points,
289  const pointField& neiCc,
290  const scalar minFaceArea,
291  const scalar maxNonOrtho,
292  const label facei
293 ) const
294 {
295  // Severe nonorthogonality threshold
296  const scalar severeNonorthogonalityThreshold =
297  ::cos(degToRad(maxNonOrtho));
298 
299  vector s = mesh_.faces()[facei].normal(points);
300  scalar magS = mag(s);
301 
302  // Check face area
303  if (magS < minFaceArea)
304  {
305  return true;
306  }
307 
308  // Check orthogonality
309  const point& ownCc = mesh_.cellCentres()[mesh_.faceOwner()[facei]];
310 
311  if (mesh_.isInternalFace(facei))
312  {
313  label nei = mesh_.faceNeighbour()[facei];
314  vector d = mesh_.cellCentres()[nei] - ownCc;
315 
316  scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
317 
318  if (dDotS < severeNonorthogonalityThreshold)
319  {
320  return true;
321  }
322  else
323  {
324  return false;
325  }
326  }
327  else
328  {
329  label patchi = mesh_.boundaryMesh().whichPatch(facei);
330 
331  if (mesh_.boundaryMesh()[patchi].coupled())
332  {
333  vector d = neiCc[facei-mesh_.nInternalFaces()] - ownCc;
334 
335  scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
336 
337  if (dDotS < severeNonorthogonalityThreshold)
338  {
339  return true;
340  }
341  else
342  {
343  return false;
344  }
345  }
346  else
347  {
348  // Collapsing normal boundary face does not cause problems with
349  // non-orthogonality
350  return false;
351  }
352  }
353 }
354 
355 
356 // Check if moving cell to new points causes it to collapse.
357 bool Foam::meshRefinement::isCollapsedCell
358 (
359  const pointField& points,
360  const scalar volFraction,
361  const label celli
362 ) const
363 {
364  scalar vol = mesh_.cells()[celli].mag(points, mesh_.faces());
365 
366  if (vol/mesh_.cellVolumes()[celli] < volFraction)
367  {
368  return true;
369  }
370  else
371  {
372  return false;
373  }
374 }
375 
376 
377 Foam::labelList Foam::meshRefinement::nearestPatch
378 (
379  const labelList& adaptPatchIDs
380 ) const
381 {
382  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
383 
384  labelList nearestAdaptPatch;
385 
386  if (adaptPatchIDs.size())
387  {
388  nearestAdaptPatch.setSize(mesh_.nFaces(), adaptPatchIDs[0]);
389 
390 
391  // Count number of faces in adaptPatchIDs
392  label nFaces = 0;
393  forAll(adaptPatchIDs, i)
394  {
395  const polyPatch& pp = patches[adaptPatchIDs[i]];
396  nFaces += pp.size();
397  }
398 
399  // Field on cells and faces.
400  List<topoDistanceData> cellData(mesh_.nCells());
401  List<topoDistanceData> faceData(mesh_.nFaces());
402 
403  // Start of changes
404  labelList patchFaces(nFaces);
405  List<topoDistanceData> patchData(nFaces);
406  nFaces = 0;
407  forAll(adaptPatchIDs, i)
408  {
409  label patchi = adaptPatchIDs[i];
410  const polyPatch& pp = patches[patchi];
411 
412  forAll(pp, i)
413  {
414  patchFaces[nFaces] = pp.start()+i;
415  patchData[nFaces] = topoDistanceData(patchi, 0);
416  nFaces++;
417  }
418  }
419 
420  // Propagate information inwards
421  FaceCellWave<topoDistanceData> deltaCalc
422  (
423  mesh_,
424  patchFaces,
425  patchData,
426  faceData,
427  cellData,
428  mesh_.globalData().nTotalCells()+1
429  );
430 
431  // And extract
432 
433  bool haveWarned = false;
434  forAll(faceData, facei)
435  {
436  if (!faceData[facei].valid(deltaCalc.data()))
437  {
438  if (!haveWarned)
439  {
441  << "Did not visit some faces, e.g. face " << facei
442  << " at " << mesh_.faceCentres()[facei] << endl
443  << "Assigning these cells to patch "
444  << adaptPatchIDs[0]
445  << endl;
446  haveWarned = true;
447  }
448  }
449  else
450  {
451  nearestAdaptPatch[facei] = faceData[facei].data();
452  }
453  }
454  }
455  else
456  {
457  // Use patch 0
458  nearestAdaptPatch.setSize(mesh_.nFaces(), 0);
459  }
460 
461  return nearestAdaptPatch;
462 }
463 
464 
465 // Returns list with for every internal face -1 or the patch they should
466 // be baffled into. Gets run after createBaffles so all the unzoned surface
467 // intersections have already been turned into baffles. (Note: zoned surfaces
468 // are not baffled at this stage)
469 // Used to remove cells by baffling all their faces and have the
470 // splitMeshRegions chuck away non used regions.
471 Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
472 (
473  const dictionary& motionDict,
474  const bool removeEdgeConnectedCells,
475  const scalarField& perpendicularAngle,
476  const labelList& globalToPatch
477 ) const
478 {
479  const labelList& cellLevel = meshCutter_.cellLevel();
480  const labelList& pointLevel = meshCutter_.pointLevel();
481  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
482 
483 
484  // Mark all points and edges on baffle patches (so not on any inlets,
485  // outlets etc.)
486  boolList isBoundaryPoint(mesh_.nPoints(), false);
487  boolList isBoundaryEdge(mesh_.nEdges(), false);
488  boolList isBoundaryFace(mesh_.nFaces(), false);
489 
490  // Fill boundary data. All elements on meshed patches get marked.
491  // Get the labels of added patches.
492  labelList adaptPatchIDs(meshedPatches());
493 
494  forAll(adaptPatchIDs, i)
495  {
496  const polyPatch& pp = patches[adaptPatchIDs[i]];
497 
498  label facei = pp.start();
499 
500  forAll(pp, j)
501  {
502  markBoundaryFace
503  (
504  facei,
505  isBoundaryFace,
506  isBoundaryEdge,
507  isBoundaryPoint
508  );
509 
510  facei++;
511  }
512  }
513 
514 
515  // Per face the nearest adaptPatch
516  const labelList nearestAdaptPatch(nearestPatch(adaptPatchIDs));
517 
518 
519  // Per internal face (boundary faces not used) the patch that the
520  // baffle should get (or -1)
521  labelList facePatch(mesh_.nFaces(), -1);
522 
523  // Swap neighbouring cell centres and cell level
524  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
525  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
526  calcNeighbourData(neiLevel, neiCc);
527 
528 
529  // Count of faces marked for baffling
530  label nBaffleFaces = 0;
531  PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
532 
533  // Count of faces not baffled since would not cause a collapse
534  label nPrevented = 0;
535 
536  if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
537  {
538  Info<< "markFacesOnProblemCells :"
539  << " Checking for edge-connected cells of highly differing sizes."
540  << endl;
541 
542  // Pick up the cells that need to be removed and (a guess for)
543  // the patch they should be patched with.
544  Map<label> problemCells
545  (
546  findEdgeConnectedProblemCells
547  (
548  perpendicularAngle,
549  globalToPatch
550  )
551  );
552 
553  // Baffle all faces of cells that need to be removed
554  forAllConstIter(Map<label>, problemCells, iter)
555  {
556  const cell& cFaces = mesh_.cells()[iter.key()];
557 
558  forAll(cFaces, i)
559  {
560  label facei = cFaces[i];
561 
562  if (facePatch[facei] == -1 && mesh_.isInternalFace(facei))
563  {
564  facePatch[facei] = nearestAdaptPatch[facei];
565  nBaffleFaces++;
566 
567  // Mark face as a 'boundary'
568  markBoundaryFace
569  (
570  facei,
571  isBoundaryFace,
572  isBoundaryEdge,
573  isBoundaryPoint
574  );
575  }
576  }
577  }
578  Info<< "markFacesOnProblemCells : Marked "
579  << returnReduce(nBaffleFaces, sumOp<label>())
580  << " additional internal faces to be converted into baffles"
581  << " due to "
582  << returnReduce(problemCells.size(), sumOp<label>())
583  << " cells edge-connected to lower level cells." << endl;
584 
585  if (debug&meshRefinement::MESH)
586  {
587  cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
588  problemCellSet.instance() = timeName();
589  Pout<< "Writing " << problemCellSet.size()
590  << " cells that are edge connected to coarser cell to set "
591  << problemCellSet.objectPath() << endl;
592  problemCellSet.write();
593  }
594  }
595 
597  (
598  mesh_,
599  isBoundaryPoint,
600  orEqOp<bool>(),
601  false // null value
602  );
603 
605  (
606  mesh_,
607  isBoundaryEdge,
608  orEqOp<bool>(),
609  false // null value
610  );
611 
613  (
614  mesh_,
615  isBoundaryFace,
616  orEqOp<bool>()
617  );
618 
619 
620  // See if checking for collapse
621  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
622 
623  // Collapse checking parameters
624  const scalar volFraction =
625  motionDict.lookupOrDefault<scalar>("minVolCollapseRatio", -1);
626 
627  const bool checkCollapse = (volFraction > 0);
628  scalar minArea = -1;
629  scalar maxNonOrtho = -1;
630 
631 
632  // Find nearest (non-baffle) surface
633  pointField newPoints;
634 
635  if (checkCollapse)
636  {
637  minArea = readScalar(motionDict.lookup("minArea"));
638  maxNonOrtho = readScalar(motionDict.lookup("maxNonOrtho"));
639 
640  Info<< "markFacesOnProblemCells :"
641  << " Deleting all-anchor surface cells only if"
642  << " snapping them violates mesh quality constraints:" << nl
643  << " snapped/original cell volume < " << volFraction << nl
644  << " face area < " << minArea << nl
645  << " non-orthogonality > " << maxNonOrtho << nl
646  << endl;
647 
648  // Construct addressing engine.
649  autoPtr<indirectPrimitivePatch> ppPtr
650  (
652  (
653  mesh_,
654  adaptPatchIDs
655  )
656  );
657  const indirectPrimitivePatch& pp = ppPtr();
658  const pointField& localPoints = pp.localPoints();
659  const labelList& meshPoints = pp.meshPoints();
660 
661  List<pointIndexHit> hitInfo;
662  labelList hitSurface;
663  surfaces_.findNearest
664  (
666  localPoints,
667  scalarField(localPoints.size(), sqr(GREAT)), // sqr of attraction
668  hitSurface,
669  hitInfo
670  );
671 
672  // Start off from current points
673  newPoints = mesh_.points();
674 
675  forAll(hitInfo, i)
676  {
677  if (hitInfo[i].hit())
678  {
679  newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
680  }
681  }
682 
683  if (debug&meshRefinement::MESH)
684  {
685  const_cast<Time&>(mesh_.time())++;
686  pointField oldPoints(mesh_.points());
687  mesh_.movePoints(newPoints);
688  Pout<< "Writing newPoints mesh to time " << timeName()
689  << endl;
690  write
691  (
692  debugType(debug),
694  mesh_.time().path()/"newPoints"
695  );
696  mesh_.movePoints(oldPoints);
697  }
698  }
699 
700 
701 
702  // For each cell count the number of anchor points that are on
703  // the boundary:
704  // 8 : check the number of (baffle) boundary faces. If 3 or more block
705  // off the cell since the cell would get squeezed down to a diamond
706  // (probably; if the 3 or more faces are unrefined (only use the
707  // anchor points))
708  // 7 : store. Used to check later on whether there are points with
709  // 3 or more of these cells. (note that on a flat surface a boundary
710  // point will only have 4 cells connected to it)
711 
712 
713  // Does cell have exactly 7 of its 8 anchor points on the boundary?
714  PackedBoolList hasSevenBoundaryAnchorPoints(mesh_.nCells());
715  // If so what is the remaining non-boundary anchor point?
716  labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
717 
718  // On-the-fly addressing storage.
719  DynamicList<label> dynFEdges;
720  DynamicList<label> dynCPoints;
721 
722  forAll(cellLevel, celli)
723  {
724  const labelList& cPoints = mesh_.cellPoints(celli, dynCPoints);
725 
726  // Get number of anchor points (pointLevel <= cellLevel)
727 
728  label nBoundaryAnchors = 0;
729  label nNonAnchorBoundary = 0;
730  label nonBoundaryAnchor = -1;
731 
732  forAll(cPoints, i)
733  {
734  label pointi = cPoints[i];
735 
736  if (pointLevel[pointi] <= cellLevel[celli])
737  {
738  // Anchor point
739  if (isBoundaryPoint[pointi])
740  {
741  nBoundaryAnchors++;
742  }
743  else
744  {
745  // Anchor point which is not on the surface
746  nonBoundaryAnchor = pointi;
747  }
748  }
749  else if (isBoundaryPoint[pointi])
750  {
751  nNonAnchorBoundary++;
752  }
753  }
754 
755  if (nBoundaryAnchors == 8)
756  {
757  const cell& cFaces = mesh_.cells()[celli];
758 
759  // Count boundary faces.
760  label nBfaces = 0;
761 
762  forAll(cFaces, cFacei)
763  {
764  if (isBoundaryFace[cFaces[cFacei]])
765  {
766  nBfaces++;
767  }
768  }
769 
770  // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
771  // We assume that this situation is where there is a single
772  // cell sticking out which would get flattened.
773 
774  // Eugene: delete cell no matter what.
775  //if (nBfaces > 1)
776  {
777  if
778  (
779  checkCollapse
780  && !isCollapsedCell(newPoints, volFraction, celli)
781  )
782  {
783  nPrevented++;
784  //Pout<< "Preventing baffling/removal of 8 anchor point"
785  // << " cell "
786  // << celli << " at " << mesh_.cellCentres()[celli]
787  // << " since new volume "
788  // << mesh_.cells()[celli].mag(newPoints, mesh_.faces())
789  // << " old volume " << mesh_.cellVolumes()[celli]
790  // << endl;
791  }
792  else
793  {
794  // Block all faces of cell
795  forAll(cFaces, cf)
796  {
797  label facei = cFaces[cf];
798 
799  if
800  (
801  facePatch[facei] == -1
802  && mesh_.isInternalFace(facei)
803  )
804  {
805  facePatch[facei] = nearestAdaptPatch[facei];
806  nBaffleFaces++;
807 
808  // Mark face as a 'boundary'
809  markBoundaryFace
810  (
811  facei,
812  isBoundaryFace,
813  isBoundaryEdge,
814  isBoundaryPoint
815  );
816  }
817  }
818  }
819  }
820  }
821  else if (nBoundaryAnchors == 7)
822  {
823  // Mark the cell. Store the (single!) non-boundary anchor point.
824  hasSevenBoundaryAnchorPoints.set(celli, 1u);
825  nonBoundaryAnchors.insert(nonBoundaryAnchor);
826  }
827  }
828 
829 
830  // Loop over all points. If a point is connected to 4 or more cells
831  // with 7 anchor points on the boundary set those cell's non-boundary faces
832  // to baffles
833 
834  DynamicList<label> dynPCells;
835 
836  forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
837  {
838  label pointi = iter.key();
839 
840  const labelList& pCells = mesh_.pointCells(pointi, dynPCells);
841 
842  // Count number of 'hasSevenBoundaryAnchorPoints' cells.
843  label n = 0;
844 
845  forAll(pCells, i)
846  {
847  if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
848  {
849  n++;
850  }
851  }
852 
853  if (n > 3)
854  {
855  // Point in danger of being what? Remove all 7-cells.
856  forAll(pCells, i)
857  {
858  label celli = pCells[i];
859 
860  if (hasSevenBoundaryAnchorPoints.get(celli) == 1u)
861  {
862  if
863  (
864  checkCollapse
865  && !isCollapsedCell(newPoints, volFraction, celli)
866  )
867  {
868  nPrevented++;
869  //Pout<< "Preventing baffling of 7 anchor cell "
870  // << celli
871  // << " at " << mesh_.cellCentres()[celli]
872  // << " since new volume "
873  // << mesh_.cells()[celli].mag
874  // (newPoints, mesh_.faces())
875  // << " old volume " << mesh_.cellVolumes()[celli]
876  // << endl;
877  }
878  else
879  {
880  const cell& cFaces = mesh_.cells()[celli];
881 
882  forAll(cFaces, cf)
883  {
884  label facei = cFaces[cf];
885 
886  if
887  (
888  facePatch[facei] == -1
889  && mesh_.isInternalFace(facei)
890  )
891  {
892  facePatch[facei] = nearestAdaptPatch[facei];
893  nBaffleFaces++;
894 
895  // Mark face as a 'boundary'
896  markBoundaryFace
897  (
898  facei,
899  isBoundaryFace,
900  isBoundaryEdge,
901  isBoundaryPoint
902  );
903  }
904  }
905  }
906  }
907  }
908  }
909  }
910 
911 
912  // Sync all. (note that pointdata and facedata not used anymore but sync
913  // anyway)
914 
916  (
917  mesh_,
918  isBoundaryPoint,
919  orEqOp<bool>(),
920  false // null value
921  );
922 
924  (
925  mesh_,
926  isBoundaryEdge,
927  orEqOp<bool>(),
928  false // null value
929  );
930 
932  (
933  mesh_,
934  isBoundaryFace,
935  orEqOp<bool>()
936  );
937 
938 
939  // Find faces with all edges on the boundary and make them baffles
940  for (label facei = 0; facei < mesh_.nInternalFaces(); facei++)
941  {
942  if (facePatch[facei] == -1)
943  {
944  const labelList& fEdges = mesh_.faceEdges(facei, dynFEdges);
945  label nFaceBoundaryEdges = 0;
946 
947  forAll(fEdges, fe)
948  {
949  if (isBoundaryEdge[fEdges[fe]])
950  {
951  nFaceBoundaryEdges++;
952  }
953  }
954 
955  if (nFaceBoundaryEdges == fEdges.size())
956  {
957  if
958  (
959  checkCollapse
960  && !isCollapsedFace
961  (
962  newPoints,
963  neiCc,
964  minArea,
965  maxNonOrtho,
966  facei
967  )
968  )
969  {
970  nPrevented++;
971  //Pout<< "Preventing baffling (to avoid collapse) of face "
972  // << facei
973  // << " with all boundary edges "
974  // << " at " << mesh_.faceCentres()[facei]
975  // << endl;
976  }
977  else
978  {
979  facePatch[facei] = nearestAdaptPatch[facei];
980  nBaffleFaces++;
981 
982  // Do NOT update boundary data since this would grow blocked
983  // faces across gaps.
984  }
985  }
986  }
987  }
988 
989  forAll(patches, patchi)
990  {
991  const polyPatch& pp = patches[patchi];
992 
993  if (pp.coupled())
994  {
995  label facei = pp.start();
996 
997  forAll(pp, i)
998  {
999  if (facePatch[facei] == -1)
1000  {
1001  const labelList& fEdges = mesh_.faceEdges(facei, dynFEdges);
1002  label nFaceBoundaryEdges = 0;
1003 
1004  forAll(fEdges, fe)
1005  {
1006  if (isBoundaryEdge[fEdges[fe]])
1007  {
1008  nFaceBoundaryEdges++;
1009  }
1010  }
1011 
1012  if (nFaceBoundaryEdges == fEdges.size())
1013  {
1014  if
1015  (
1016  checkCollapse
1017  && !isCollapsedFace
1018  (
1019  newPoints,
1020  neiCc,
1021  minArea,
1022  maxNonOrtho,
1023  facei
1024  )
1025  )
1026  {
1027  nPrevented++;
1028  //Pout<< "Preventing baffling of coupled face "
1029  // << facei
1030  // << " with all boundary edges "
1031  // << " at " << mesh_.faceCentres()[facei]
1032  // << endl;
1033  }
1034  else
1035  {
1036  facePatch[facei] = nearestAdaptPatch[facei];
1037  if (isMasterFace[facei])
1038  {
1039  nBaffleFaces++;
1040  }
1041 
1042  // Do NOT update boundary data since this would grow
1043  // blocked faces across gaps.
1044  }
1045  }
1046  }
1047 
1048  facei++;
1049  }
1050  }
1051  }
1052 
1053 
1054  // Because of isCollapsedFace one side can decide not to baffle whereas
1055  // the other side does so sync. Baffling is prefered over not baffling.
1056  if (checkCollapse) // Or always?
1057  {
1059  (
1060  mesh_,
1061  facePatch,
1062  maxEqOp<label>()
1063  );
1064  }
1065 
1066  Info<< "markFacesOnProblemCells : marked "
1067  << returnReduce(nBaffleFaces, sumOp<label>())
1068  << " additional internal faces to be converted into baffles."
1069  << endl;
1070 
1071  if (checkCollapse)
1072  {
1073  Info<< "markFacesOnProblemCells : prevented "
1074  << returnReduce(nPrevented, sumOp<label>())
1075  << " internal faces from getting converted into baffles."
1076  << endl;
1077  }
1078 
1079  return facePatch;
1080 }
1081 
1082 
1083 // Mark faces to be baffled to prevent snapping problems. Does
1084 // test to find nearest surface and checks which faces would get squashed.
1085 Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
1086 (
1087  const snapParameters& snapParams,
1088  const dictionary& motionDict
1089 ) const
1090 {
1091  pointField oldPoints(mesh_.points());
1092 
1093  // Repeat (most of) snappySnapDriver::doSnap
1094  {
1095  labelList adaptPatchIDs(meshedPatches());
1096 
1097  // Construct addressing engine.
1098  autoPtr<indirectPrimitivePatch> ppPtr
1099  (
1101  (
1102  mesh_,
1103  adaptPatchIDs
1104  )
1105  );
1106  indirectPrimitivePatch& pp = ppPtr();
1107 
1108  // Distance to attract to nearest feature on surface
1109  const scalarField snapDist
1110  (
1111  snappySnapDriver::calcSnapDistance(mesh_, snapParams, pp)
1112  );
1113 
1114 
1115  // Construct iterative mesh mover.
1116  Info<< "Constructing mesh displacer ..." << endl;
1117  Info<< "Using mesh parameters " << motionDict << nl << endl;
1118 
1119  const pointMesh& pMesh = pointMesh::New(mesh_);
1120 
1121  motionSmoother meshMover
1122  (
1123  mesh_,
1124  pp,
1125  adaptPatchIDs,
1126  meshRefinement::makeDisplacementField(pMesh, adaptPatchIDs)(),
1127  motionDict
1128  );
1129 
1130 
1131  // Check initial mesh
1132  Info<< "Checking initial mesh ..." << endl;
1133  labelHashSet wrongFaces(mesh_.nFaces()/100);
1134  motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
1135  const label nInitErrors = returnReduce
1136  (
1137  wrongFaces.size(),
1138  sumOp<label>()
1139  );
1140 
1141  Info<< "Detected " << nInitErrors << " illegal faces"
1142  << " (concave, zero area or negative cell pyramid volume)"
1143  << endl;
1144 
1145 
1146  Info<< "Checked initial mesh in = "
1147  << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl;
1148 
1149  // Pre-smooth patch vertices (so before determining nearest)
1151  (
1152  *this,
1153  snapParams,
1154  nInitErrors,
1155  List<labelPair>(0), //baffles
1156  meshMover
1157  );
1158 
1159  pointField nearestPoint;
1160  vectorField nearestNormal;
1161  const vectorField disp
1162  (
1164  (
1165  *this,
1166  snapDist, // attraction
1167  pp,
1168  nearestPoint,
1169  nearestNormal
1170  )
1171  );
1172 
1173  const labelList& meshPoints = pp.meshPoints();
1174 
1175  pointField newPoints(mesh_.points());
1176  forAll(meshPoints, i)
1177  {
1178  newPoints[meshPoints[i]] += disp[i];
1179  }
1180 
1182  (
1183  mesh_,
1184  newPoints,
1185  minMagSqrEqOp<point>(), // combine op
1186  vector(GREAT, GREAT, GREAT) // null value (note: cannot use VGREAT)
1187  );
1188 
1189  mesh_.movePoints(newPoints);
1190  }
1191 
1192 
1193  // Per face the nearest adaptPatch
1194  const labelList nearestAdaptPatch(nearestPatch(meshedPatches()));
1195 
1196  // Per face (internal or coupled!) the patch that the
1197  // baffle should get (or -1).
1198  labelList facePatch(mesh_.nFaces(), -1);
1199  // Count of baffled faces
1200  label nBaffleFaces = 0;
1201 
1202  {
1203  faceSet wrongFaces(mesh_, "wrongFaces", 100);
1204  {
1205  //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
1206 
1207  // Just check the errors from squashing
1208  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1209 
1210  const labelList allFaces(identity(mesh_.nFaces()));
1211  label nWrongFaces = 0;
1212 
1213  //const scalar minV(readScalar(motionDict.lookup("minVol", true)));
1214  //if (minV > -GREAT)
1215  //{
1216  // polyMeshGeometry::checkFacePyramids
1217  // (
1218  // false,
1219  // minV,
1220  // mesh_,
1221  // mesh_.cellCentres(),
1222  // mesh_.points(),
1223  // allFaces,
1224  // List<labelPair>(0),
1225  // &wrongFaces
1226  // );
1227  //
1228  // label nNewWrongFaces = returnReduce
1229  // (
1230  // wrongFaces.size(),
1231  // sumOp<label>()
1232  // );
1233  //
1234  // Info<< " faces with pyramid volume < "
1235  // << setw(5) << minV
1236  // << " m^3 : "
1237  // << nNewWrongFaces-nWrongFaces << endl;
1238  //
1239  // nWrongFaces = nNewWrongFaces;
1240  //}
1241 
1242  scalar minArea(readScalar(motionDict.lookup("minArea")));
1243  if (minArea > -SMALL)
1244  {
1246  (
1247  false,
1248  minArea,
1249  mesh_,
1250  mesh_.faceAreas(),
1251  allFaces,
1252  &wrongFaces
1253  );
1254 
1255  label nNewWrongFaces = returnReduce
1256  (
1257  wrongFaces.size(),
1258  sumOp<label>()
1259  );
1260 
1261  Info<< " faces with area < "
1262  << setw(5) << minArea
1263  << " m^2 : "
1264  << nNewWrongFaces-nWrongFaces << endl;
1265 
1266  nWrongFaces = nNewWrongFaces;
1267  }
1268 
1269  scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
1270  if (minDet > -1)
1271  {
1273  (
1274  false,
1275  minDet,
1276  mesh_,
1277  mesh_.faceAreas(),
1278  allFaces,
1279  polyMeshGeometry::affectedCells(mesh_, allFaces),
1280  &wrongFaces
1281  );
1282 
1283  label nNewWrongFaces = returnReduce
1284  (
1285  wrongFaces.size(),
1286  sumOp<label>()
1287  );
1288 
1289  Info<< " faces on cells with determinant < "
1290  << setw(5) << minDet << " : "
1291  << nNewWrongFaces-nWrongFaces << endl;
1292 
1293  nWrongFaces = nNewWrongFaces;
1294  }
1295  }
1296 
1297 
1298  forAllConstIter(faceSet, wrongFaces, iter)
1299  {
1300  label patchi = mesh_.boundaryMesh().whichPatch(iter.key());
1301 
1302  if (patchi == -1 || mesh_.boundaryMesh()[patchi].coupled())
1303  {
1304  facePatch[iter.key()] = nearestAdaptPatch[iter.key()];
1305  nBaffleFaces++;
1306 
1307  //Pout<< " " << iter.key()
1308  // //<< " on patch " << mesh_.boundaryMesh()[patchi].name()
1309  // << " is destined for patch " << facePatch[iter.key()]
1310  // << endl;
1311  }
1312  }
1313  }
1314 
1315 
1316  // Restore points.
1317  mesh_.movePoints(oldPoints);
1318 
1319 
1320  Info<< "markFacesOnProblemCellsGeometric : marked "
1321  << returnReduce(nBaffleFaces, sumOp<label>())
1322  << " additional internal and coupled faces"
1323  << " to be converted into baffles." << endl;
1324 
1326  (
1327  mesh_,
1328  facePatch,
1329  maxEqOp<label>()
1330  );
1331 
1332  return facePatch;
1333 }
1334 
1335 
1336 // ************************************************************************* //
static bool checkMesh(const bool report, const polyMesh &mesh, const dictionary &dict, labelHashSet &wrongFaces)
Check mesh with mesh settings in dict. Collects incorrect faces.
List< labelList > labelListList
A List of labelList.
Definition: labelList.H:57
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
const vectorField & faceAreas() const
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Unit conversion functions.
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:76
const labelListList & cellPoints() const
const searchableSurfaces & geometry() const
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
static writeType writeLevel()
Get/set write level.
labelList identity(const label len)
Create identity map (map[i] == i) of given length.
Definition: ListOps.C:104
const labelIOList & cellLevel() const
Definition: hexRef8.H:389
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
const vectorField & faceCentres() const
bool isInternalFace(const label faceIndex) const
Return true if given face label is internal to the mesh.
PrimitivePatch< face, IndirectList, const pointField & > indirectPrimitivePatch
Foam::indirectPrimitivePatch.
const cellList & cells() const
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:421
static PackedBoolList getMasterFaces(const polyMesh &)
Get per face whether it is uncoupled or a master of a.
Definition: syncTools.C:153
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
label nCells() const
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:979
static const pointMesh & New(const polyMesh &mesh)
List< bool > boolList
Bool container classes.
Definition: boolList.H:50
HashSet< label, Hash< label > > labelHashSet
A HashSet with label keys.
Definition: HashSet.H:210
bool write() const
Write mesh and all data.
gmvFile<< "tracers "<< particles.size()<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().x()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().y()<< " ";}gmvFile<< nl;forAllConstIter(Cloud< passiveParticle >, particles, iter){gmvFile<< iter().position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
vectorField pointField
pointField is a vectorField.
Definition: pointFieldFwd.H:42
const labelListList & pointCells() const
dimensionedScalar cos(const dimensionedScalar &ds)
static void syncEdgeList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh edges.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
void findNearest(const labelList &surfacesToTest, const pointField &samples, const scalarField &nearestDistSqr, labelList &surfaces, List< pointIndexHit > &) const
Find nearest point on surfaces.
const vectorField & cellCentres() const
List< label > labelList
A List of labels.
Definition: labelList.H:56
static void preSmoothPatch(const meshRefinement &meshRefiner, const snapParameters &snapParams, const label nInitErrors, const List< labelPair > &baffles, motionSmoother &)
Smooth the mesh (patch and internal) to increase visibility.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
static labelList getUnnamedSurfaces(const PtrList< surfaceZonesInfo > &surfList)
Get indices of unnamed surfaces (surfaces without faceZoneName)
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
static scalarField calcSnapDistance(const fvMesh &mesh, const snapParameters &snapParams, const indirectPrimitivePatch &)
Calculate edge length per patch point.
const globalMeshData & globalData() const
Return parallel info.
Definition: polyMesh.C:1140
prefixOSstream Pout(cout,"Pout")
Definition: IOstreams.H:53
dimensionedScalar sin(const dimensionedScalar &ds)
Istream and Ostream manipulators taking arguments.
static const char nl
Definition: Ostream.H:262
double cpuTimeIncrement() const
Return CPU time (in seconds) since last call to cpuTimeIncrement()
Definition: cpuTime.C:74
static bool checkFaceArea(const bool report, const scalar minArea, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, labelHashSet *setPtr)
Small faces.
label nEdges() const
const scalarField & cellVolumes() const
static autoPtr< indirectPrimitivePatch > makePatch(const polyMesh &, const labelList &)
Create patch from set of patches.
const labelIOList & pointLevel() const
Definition: hexRef8.H:394
virtual tmp< scalarField > movePoints(const pointField &)
Move points, returns volumes swept by faces in motion.
Definition: fvMesh.C:721
label globalRegion(const label surfI, const label regionI) const
From surface and region on surface to global region.
label nFaces() const
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:295
const PtrList< surfaceZonesInfo > & surfZones() const
label patchi
virtual const labelList & faceNeighbour() const
Return face neighbour.
Definition: polyMesh.C:1023
vector point
Point is a vector.
Definition: point.H:41
#define WarningInFunction
Report a warning using Foam::Warning.
label whichPatch(const label faceIndex) const
Return patch index for a given face label.
static void syncPointList(const polyMesh &, List< T > &, const CombineOp &cop, const T &nullValue, const TransformOp &top)
Synchronize values on all mesh points.
static labelList affectedCells(const polyMesh &, const labelList &changedFaces)
Helper function: get affected cells from faces.
const labelList & surfaces() const
label nPoints() const
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
virtual Ostream & write(const token &)=0
Write next token to stream.
Omanip< int > setw(const int i)
Definition: IOmanip.H:199
static tmp< pointVectorField > makeDisplacementField(const pointMesh &pMesh, const labelList &adaptPatchIDs)
Helper function to make a pointVectorField with correct.
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1017
virtual const faceList & faces() const
Return raw faces.
Definition: polyMesh.C:1004
fileName path() const
Return path.
Definition: Time.H:269
const labelListList & faceEdges() const
labelList meshedPatches() const
Get patchIDs for patches added in addMeshedPatch.
static bool checkCellDeterminant(const bool report, const scalar minDet, const polyMesh &, const vectorField &faceAreas, const labelList &checkFaces, const labelList &affectedCells, labelHashSet *setPtr)
Area of internal faces v.s. boundary faces.
label nTotalCells() const
Return total number of cells in decomposed mesh.
label nInternalFaces() const
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:243
word timeName() const
Replacement for Time::timeName() : return oldInstance (if.
A HashTable to objects of type <T> with a label key.
Definition: Map.H:49
static vectorField calcNearestSurface(const meshRefinement &meshRefiner, const scalarField &snapDist, const indirectPrimitivePatch &, pointField &nearestPoint, vectorField &nearestNormal)
Per patch point calculate point on nearest surface. Set as.