FreeFOAM The Cross-Platform CFD Toolkit
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) 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"
27 #include <finiteVolume/fvMesh.H>
28 #include <OpenFOAM/syncTools.H>
29 #include <OpenFOAM/Time.H>
31 #include <meshTools/pointSet.H>
32 #include <meshTools/faceSet.H>
34 #include <OpenFOAM/OFstream.H>
35 #include <meshTools/cellSet.H>
38 #include <OpenFOAM/IOmanip.H>
39 
40 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
41 
42 void Foam::meshRefinement::markBoundaryFace
43 (
44  const label faceI,
45  boolList& isBoundaryFace,
46  boolList& isBoundaryEdge,
47  boolList& isBoundaryPoint
48 ) const
49 {
50  isBoundaryFace[faceI] = true;
51 
52  const labelList& fEdges = mesh_.faceEdges(faceI);
53 
54  forAll(fEdges, fp)
55  {
56  isBoundaryEdge[fEdges[fp]] = true;
57  }
58 
59  const face& f = mesh_.faces()[faceI];
60 
61  forAll(f, fp)
62  {
63  isBoundaryPoint[f[fp]] = true;
64  }
65 }
66 
67 
68 void Foam::meshRefinement::findNearest
69 (
70  const labelList& meshFaces,
71  List<pointIndexHit>& nearestInfo,
72  labelList& nearestSurface,
73  labelList& nearestRegion,
74  vectorField& nearestNormal
75 ) const
76 {
77  pointField fc(meshFaces.size());
78  forAll(meshFaces, i)
79  {
80  fc[i] = mesh_.faceCentres()[meshFaces[i]];
81  }
82 
83  const labelList allSurfaces(identity(surfaces_.surfaces().size()));
84 
85  surfaces_.findNearest
86  (
87  allSurfaces,
88  fc,
89  scalarField(fc.size(), sqr(GREAT)), // sqr of attraction
90  nearestSurface,
91  nearestInfo
92  );
93 
94  // Do normal testing per surface.
95  nearestNormal.setSize(nearestInfo.size());
96  nearestRegion.setSize(nearestInfo.size());
97 
98  forAll(allSurfaces, surfI)
99  {
100  DynamicList<pointIndexHit> localHits;
101 
102  forAll(nearestSurface, i)
103  {
104  if (nearestSurface[i] == surfI)
105  {
106  localHits.append(nearestInfo[i]);
107  }
108  }
109 
110  label geomI = surfaces_.surfaces()[surfI];
111 
112  pointField localNormals;
113  surfaces_.geometry()[geomI].getNormal(localHits, localNormals);
114 
115  labelList localRegion;
116  surfaces_.geometry()[geomI].getRegion(localHits, localRegion);
117 
118  label localI = 0;
119  forAll(nearestSurface, i)
120  {
121  if (nearestSurface[i] == surfI)
122  {
123  nearestNormal[i] = localNormals[localI];
124  nearestRegion[i] = localRegion[localI];
125  localI++;
126  }
127  }
128  }
129 }
130 
131 
132 Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
133 (
134  const scalarField& perpendicularAngle,
135  const labelList& globalToPatch
136 ) const
137 {
138  // Construct addressing engine from all patches added for meshing.
139  autoPtr<indirectPrimitivePatch> ppPtr
140  (
142  (
143  mesh_,
144  meshedPatches()
145  )
146  );
147  const indirectPrimitivePatch& pp = ppPtr();
148 
149 
150  // 1. Collect faces to test
151  // ~~~~~~~~~~~~~~~~~~~~~~~~
152 
153  DynamicList<label> candidateFaces(pp.size()/20);
154 
155  const labelListList& edgeFaces = pp.edgeFaces();
156 
157  const labelList& cellLevel = meshCutter_.cellLevel();
158 
159  forAll(edgeFaces, edgeI)
160  {
161  const labelList& eFaces = edgeFaces[edgeI];
162 
163  if (eFaces.size() == 2)
164  {
165  label face0 = pp.addressing()[eFaces[0]];
166  label face1 = pp.addressing()[eFaces[1]];
167 
168  label cell0 = mesh_.faceOwner()[face0];
169  label cell1 = mesh_.faceOwner()[face1];
170 
171  if (cellLevel[cell0] > cellLevel[cell1])
172  {
173  // cell0 smaller.
174  const vector& n0 = pp.faceNormals()[eFaces[0]];
175  const vector& n1 = pp.faceNormals()[eFaces[1]];
176 
177  if (mag(n0 & n1) < 0.1)
178  {
179  candidateFaces.append(face0);
180  }
181  }
182  else if (cellLevel[cell1] > cellLevel[cell0])
183  {
184  // cell1 smaller.
185  const vector& n0 = pp.faceNormals()[eFaces[0]];
186  const vector& n1 = pp.faceNormals()[eFaces[1]];
187 
188  if (mag(n0 & n1) < 0.1)
189  {
190  candidateFaces.append(face1);
191  }
192  }
193  }
194  }
195  candidateFaces.shrink();
196 
197  Info<< "Testing " << returnReduce(candidateFaces.size(), sumOp<label>())
198  << " faces on edge-connected cells of differing level."
199  << endl;
200 
201  if (debug)
202  {
203  faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
204  Pout<< "Writing " << fSet.size()
205  << " with problematic topology to faceSet "
206  << fSet.objectPath() << endl;
207  fSet.write();
208  }
209 
210 
211  // 2. Find nearest surface on candidate faces
212  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
213 
214  List<pointIndexHit> nearestInfo;
215  labelList nearestSurface;
216  labelList nearestRegion;
217  vectorField nearestNormal;
218  findNearest
219  (
220  candidateFaces,
221  nearestInfo,
222  nearestSurface,
223  nearestRegion,
224  nearestNormal
225  );
226 
227 
228  // 3. Test angle to surface
229  // ~~~~~~~~~~~~~~~~~~~~~~~~
230 
231  Map<label> candidateCells(candidateFaces.size());
232 
233  faceSet perpFaces(mesh_, "perpendicularFaces", pp.size()/100);
234 
235  forAll(candidateFaces, i)
236  {
237  label faceI = candidateFaces[i];
238 
239  vector n = mesh_.faceAreas()[faceI];
240  n /= mag(n);
241 
242  label region = surfaces_.globalRegion
243  (
244  nearestSurface[i],
245  nearestRegion[i]
246  );
247 
248  scalar angle =
249  perpendicularAngle[region]
250  / 180.0
252 
253  if (angle >= 0)
254  {
255  if (mag(n & nearestNormal[i]) < Foam::sin(angle))
256  {
257  perpFaces.insert(faceI);
258  candidateCells.insert
259  (
260  mesh_.faceOwner()[faceI],
261  globalToPatch[region]
262  );
263  }
264  }
265  }
266 
267  if (debug)
268  {
269  Pout<< "Writing " << perpFaces.size()
270  << " faces that are perpendicular to the surface to set "
271  << perpFaces.objectPath() << endl;
272  perpFaces.write();
273  }
274  return candidateCells;
275 }
276 
277 
278 // Check if moving face to new points causes a 'collapsed' face.
279 // Uses new point position only for the face, not the neighbouring
280 // cell centres
281 bool Foam::meshRefinement::isCollapsedFace
282 (
283  const pointField& points,
284  const pointField& neiCc,
285  const scalar minFaceArea,
286  const scalar maxNonOrtho,
287  const label faceI
288 ) const
289 {
290  vector s = mesh_.faces()[faceI].normal(points);
291  scalar magS = mag(s);
292 
293  // Check face area
294  if (magS < minFaceArea)
295  {
296  return true;
297  }
298 
299  // Check orthogonality
300  const point& ownCc = mesh_.cellCentres()[mesh_.faceOwner()[faceI]];
301 
302  if (mesh_.isInternalFace(faceI))
303  {
304  label nei = mesh_.faceNeighbour()[faceI];
305  vector d = ownCc - mesh_.cellCentres()[nei];
306 
307  scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
308 
309  if (dDotS < maxNonOrtho)
310  {
311  return true;
312  }
313  else
314  {
315  return false;
316  }
317  }
318  else
319  {
320  label patchI = mesh_.boundaryMesh().whichPatch(faceI);
321 
322  if (mesh_.boundaryMesh()[patchI].coupled())
323  {
324  vector d = ownCc - neiCc[faceI-mesh_.nInternalFaces()];
325 
326  scalar dDotS = (d & s)/(mag(d)*magS + VSMALL);
327 
328  if (dDotS < maxNonOrtho)
329  {
330  return true;
331  }
332  else
333  {
334  return false;
335  }
336  }
337  else
338  {
339  // Collapsing normal boundary face does not cause problems with
340  // non-orthogonality
341  return true;
342  }
343  }
344 }
345 
346 
347 // Check if moving cell to new points causes it to collapse.
348 bool Foam::meshRefinement::isCollapsedCell
349 (
350  const pointField& points,
351  const scalar volFraction,
352  const label cellI
353 ) const
354 {
355  scalar vol = mesh_.cells()[cellI].mag(points, mesh_.faces());
356 
357  if (vol/mesh_.cellVolumes()[cellI] < volFraction)
358  {
359  return true;
360  }
361  else
362  {
363  return false;
364  }
365 }
366 
367 
368 // Returns list with for every internal face -1 or the patch they should
369 // be baffled into. Gets run after createBaffles so all the unzoned surface
370 // intersections have already been turned into baffles. (Note: zoned surfaces
371 // are not baffled at this stage)
372 // Used to remove cells by baffling all their faces and have the
373 // splitMeshRegions chuck away non used regions.
374 Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
375 (
376  const dictionary& motionDict,
377  const bool removeEdgeConnectedCells,
378  const scalarField& perpendicularAngle,
379  const labelList& globalToPatch
380 ) const
381 {
382  const labelList& cellLevel = meshCutter_.cellLevel();
383  const labelList& pointLevel = meshCutter_.pointLevel();
384  const polyBoundaryMesh& patches = mesh_.boundaryMesh();
385 
386  // Per internal face (boundary faces not used) the patch that the
387  // baffle should get (or -1)
388  labelList facePatch(mesh_.nFaces(), -1);
389 
390  // Mark all points and edges on baffle patches (so not on any inlets,
391  // outlets etc.)
392  boolList isBoundaryPoint(mesh_.nPoints(), false);
393  boolList isBoundaryEdge(mesh_.nEdges(), false);
394  boolList isBoundaryFace(mesh_.nFaces(), false);
395 
396  // Fill boundary data. All elements on meshed patches get marked.
397  // Get the labels of added patches.
398  labelList adaptPatchIDs(meshedPatches());
399 
400  forAll(adaptPatchIDs, i)
401  {
402  label patchI = adaptPatchIDs[i];
403 
404  const polyPatch& pp = patches[patchI];
405 
406  label faceI = pp.start();
407 
408  forAll(pp, j)
409  {
410  markBoundaryFace
411  (
412  faceI,
413  isBoundaryFace,
414  isBoundaryEdge,
415  isBoundaryPoint
416  );
417 
418  faceI++;
419  }
420  }
421 
422  // Swap neighbouring cell centres and cell level
423  labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
424  pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
425  calcNeighbourData(neiLevel, neiCc);
426 
427 
428  // Count of faces marked for baffling
429  label nBaffleFaces = 0;
430 
431  // Count of faces not baffled since would not cause a collapse
432  label nPrevented = 0;
433 
434  if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
435  {
436  Info<< "markFacesOnProblemCells :"
437  << " Checking for edge-connected cells of highly differing sizes."
438  << endl;
439 
440  // Pick up the cells that need to be removed and (a guess for)
441  // the patch they should be patched with.
442  Map<label> problemCells
443  (
444  findEdgeConnectedProblemCells
445  (
446  perpendicularAngle,
447  globalToPatch
448  )
449  );
450 
451  // Baffle all faces of cells that need to be removed
452  forAllConstIter(Map<label>, problemCells, iter)
453  {
454  const cell& cFaces = mesh_.cells()[iter.key()];
455 
456  forAll(cFaces, i)
457  {
458  label faceI = cFaces[i];
459 
460  if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
461  {
462  facePatch[faceI] = getBafflePatch(facePatch, faceI);
463  nBaffleFaces++;
464 
465  // Mark face as a 'boundary'
466  markBoundaryFace
467  (
468  faceI,
469  isBoundaryFace,
470  isBoundaryEdge,
471  isBoundaryPoint
472  );
473  }
474  }
475  }
476  Info<< "markFacesOnProblemCells : Marked "
477  << returnReduce(nBaffleFaces, sumOp<label>())
478  << " additional internal faces to be converted into baffles"
479  << " due to "
480  << returnReduce(problemCells.size(), sumOp<label>())
481  << " cells edge-connected to lower level cells." << endl;
482 
483  if (debug)
484  {
485  cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
486  Pout<< "Writing " << problemCellSet.size()
487  << " cells that are edge connected to coarser cell to set "
488  << problemCellSet.objectPath() << endl;
489  problemCellSet.write();
490  }
491  }
492 
494  (
495  mesh_,
496  isBoundaryPoint,
497  orEqOp<bool>(),
498  false, // null value
499  false // no separation
500  );
501 
503  (
504  mesh_,
505  isBoundaryEdge,
506  orEqOp<bool>(),
507  false, // null value
508  false // no separation
509  );
510 
512  (
513  mesh_,
514  isBoundaryFace,
515  orEqOp<bool>(),
516  false // no separation
517  );
518 
519 
520  // See if checking for collapse
521  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
522 
523  // Collapse checking parameters
524  scalar volFraction = -1;
525  if (motionDict.found("minVolCollapseRatio"))
526  {
527  volFraction = readScalar(motionDict.lookup("minVolCollapseRatio"));
528  }
529  const bool checkCollapse = (volFraction > 0);
530  scalar minArea = -1;
531  scalar maxNonOrtho = -1;
532 
533 
534  // Find nearest (non-baffle) surface
535  pointField newPoints;
536 
537  if (checkCollapse)
538  {
539  minArea = readScalar(motionDict.lookup("minArea"));
540  maxNonOrtho = readScalar(motionDict.lookup("maxNonOrtho"));
541 
542  Info<< "markFacesOnProblemCells :"
543  << " Deleting all-anchor surface cells only if"
544  << "snapping them violates mesh quality constraints:" << nl
545  << " snapped/original cell volume < " << volFraction << nl
546  << " face area < " << minArea << nl
547  << " non-orthogonality > " << maxNonOrtho << nl
548  << endl;
549 
550  // Construct addressing engine.
551  autoPtr<indirectPrimitivePatch> ppPtr
552  (
554  (
555  mesh_,
556  adaptPatchIDs
557  )
558  );
559  const indirectPrimitivePatch& pp = ppPtr();
560  const pointField& localPoints = pp.localPoints();
561  const labelList& meshPoints = pp.meshPoints();
562 
563  List<pointIndexHit> hitInfo;
564  labelList hitSurface;
565  surfaces_.findNearest
566  (
567  surfaces_.getUnnamedSurfaces(),
568  localPoints,
569  scalarField(localPoints.size(), sqr(GREAT)), // sqr of attraction
570  hitSurface,
571  hitInfo
572  );
573 
574  // Start of from current points
575  newPoints = mesh_.points();
576 
577  forAll(hitInfo, i)
578  {
579  if (hitInfo[i].hit())
580  {
581  newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
582  }
583  }
584  }
585 
586 
587 
588  // For each cell count the number of anchor points that are on
589  // the boundary:
590  // 8 : check the number of (baffle) boundary faces. If 3 or more block
591  // off the cell since the cell would get squeezed down to a diamond
592  // (probably; if the 3 or more faces are unrefined (only use the
593  // anchor points))
594  // 7 : store. Used to check later on whether there are points with
595  // 3 or more of these cells. (note that on a flat surface a boundary
596  // point will only have 4 cells connected to it)
597 
598 
599  // Does cell have exactly 7 of its 8 anchor points on the boundary?
600  PackedBoolList hasSevenBoundaryAnchorPoints(mesh_.nCells());
601  // If so what is the remaining non-boundary anchor point?
602  labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
603 
604  // On-the-fly addressing storage.
605  DynamicList<label> dynFEdges;
606  DynamicList<label> dynCPoints;
607 
608  forAll(cellLevel, cellI)
609  {
610  const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
611 
612  // Get number of anchor points (pointLevel <= cellLevel)
613 
614  label nBoundaryAnchors = 0;
615  label nNonAnchorBoundary = 0;
616  label nonBoundaryAnchor = -1;
617 
618  forAll(cPoints, i)
619  {
620  label pointI = cPoints[i];
621 
622  if (pointLevel[pointI] <= cellLevel[cellI])
623  {
624  // Anchor point
625  if (isBoundaryPoint[pointI])
626  {
627  nBoundaryAnchors++;
628  }
629  else
630  {
631  // Anchor point which is not on the surface
632  nonBoundaryAnchor = pointI;
633  }
634  }
635  else if (isBoundaryPoint[pointI])
636  {
637  nNonAnchorBoundary++;
638  }
639  }
640 
641  if (nBoundaryAnchors == 8)
642  {
643  const cell& cFaces = mesh_.cells()[cellI];
644 
645  // Count boundary faces.
646  label nBfaces = 0;
647 
648  forAll(cFaces, cFaceI)
649  {
650  if (isBoundaryFace[cFaces[cFaceI]])
651  {
652  nBfaces++;
653  }
654  }
655 
656  // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
657  // We assume that this situation is where there is a single
658  // cell sticking out which would get flattened.
659 
660  // Eugene: delete cell no matter what.
661  //if (nBfaces > 1)
662  {
663  if
664  (
665  checkCollapse
666  && !isCollapsedCell(newPoints, volFraction, cellI)
667  )
668  {
669  nPrevented++;
670  //Pout<< "Preventing collapse of 8 anchor point cell "
671  // << cellI << " at " << mesh_.cellCentres()[cellI]
672  // << " since new volume "
673  // << mesh_.cells()[cellI].mag(newPoints, mesh_.faces())
674  // << " old volume " << mesh_.cellVolumes()[cellI]
675  // << endl;
676  }
677  else
678  {
679  // Block all faces of cell
680  forAll(cFaces, cf)
681  {
682  label faceI = cFaces[cf];
683 
684  if
685  (
686  facePatch[faceI] == -1
687  && mesh_.isInternalFace(faceI)
688  )
689  {
690  facePatch[faceI] = getBafflePatch(facePatch, faceI);
691  nBaffleFaces++;
692 
693  // Mark face as a 'boundary'
694  markBoundaryFace
695  (
696  faceI,
697  isBoundaryFace,
698  isBoundaryEdge,
699  isBoundaryPoint
700  );
701  }
702  }
703  }
704  }
705  }
706  else if (nBoundaryAnchors == 7)
707  {
708  // Mark the cell. Store the (single!) non-boundary anchor point.
709  hasSevenBoundaryAnchorPoints.set(cellI, 1u);
710  nonBoundaryAnchors.insert(nonBoundaryAnchor);
711  }
712  }
713 
714 
715  // Loop over all points. If a point is connected to 4 or more cells
716  // with 7 anchor points on the boundary set those cell's non-boundary faces
717  // to baffles
718 
719  DynamicList<label> dynPCells;
720 
721  forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
722  {
723  label pointI = iter.key();
724 
725  const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
726 
727  // Count number of 'hasSevenBoundaryAnchorPoints' cells.
728  label n = 0;
729 
730  forAll(pCells, i)
731  {
732  if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
733  {
734  n++;
735  }
736  }
737 
738  if (n > 3)
739  {
740  // Point in danger of being what? Remove all 7-cells.
741  forAll(pCells, i)
742  {
743  label cellI = pCells[i];
744 
745  if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
746  {
747  if
748  (
749  checkCollapse
750  && !isCollapsedCell(newPoints, volFraction, cellI)
751  )
752  {
753  nPrevented++;
754  //Pout<< "Preventing collapse of 7 anchor cell "
755  // << cellI
756  // << " at " << mesh_.cellCentres()[cellI]
757  // << " since new volume "
758  // << mesh_.cells()[cellI].mag
759  // (newPoints, mesh_.faces())
760  // << " old volume " << mesh_.cellVolumes()[cellI]
761  // << endl;
762  }
763  else
764  {
765  const cell& cFaces = mesh_.cells()[cellI];
766 
767  forAll(cFaces, cf)
768  {
769  label faceI = cFaces[cf];
770 
771  if
772  (
773  facePatch[faceI] == -1
774  && mesh_.isInternalFace(faceI)
775  )
776  {
777  facePatch[faceI] = getBafflePatch
778  (
779  facePatch,
780  faceI
781  );
782  nBaffleFaces++;
783 
784  // Mark face as a 'boundary'
785  markBoundaryFace
786  (
787  faceI,
788  isBoundaryFace,
789  isBoundaryEdge,
790  isBoundaryPoint
791  );
792  }
793  }
794  }
795  }
796  }
797  }
798  }
799 
800 
801  // Sync all. (note that pointdata and facedata not used anymore but sync
802  // anyway)
803 
805  (
806  mesh_,
807  isBoundaryPoint,
808  orEqOp<bool>(),
809  false, // null value
810  false // no separation
811  );
812 
814  (
815  mesh_,
816  isBoundaryEdge,
817  orEqOp<bool>(),
818  false, // null value
819  false // no separation
820  );
821 
823  (
824  mesh_,
825  isBoundaryFace,
826  orEqOp<bool>(),
827  false // no separation
828  );
829 
830 
831  // Find faces with all edges on the boundary and make them baffles
832  for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
833  {
834  if (facePatch[faceI] == -1)
835  {
836  const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
837  label nFaceBoundaryEdges = 0;
838 
839  forAll(fEdges, fe)
840  {
841  if (isBoundaryEdge[fEdges[fe]])
842  {
843  nFaceBoundaryEdges++;
844  }
845  }
846 
847  if (nFaceBoundaryEdges == fEdges.size())
848  {
849  if
850  (
851  checkCollapse
852  && !isCollapsedFace
853  (
854  newPoints,
855  neiCc,
856  minArea,
857  maxNonOrtho,
858  faceI
859  )
860  )
861  {
862  nPrevented++;
863  //Pout<< "Preventing collapse of face " << faceI
864  // << " with all boundary edges "
865  // << " at " << mesh_.faceCentres()[faceI]
866  // << endl;
867  }
868  else
869  {
870  facePatch[faceI] = getBafflePatch(facePatch, faceI);
871  nBaffleFaces++;
872 
873  // Do NOT update boundary data since this would grow blocked
874  // faces across gaps.
875  }
876  }
877  }
878  }
879 
880  forAll(patches, patchI)
881  {
882  const polyPatch& pp = patches[patchI];
883 
884  if (pp.coupled())
885  {
886  label faceI = pp.start();
887 
888  forAll(pp, i)
889  {
890  if (facePatch[faceI] == -1)
891  {
892  const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
893  label nFaceBoundaryEdges = 0;
894 
895  forAll(fEdges, fe)
896  {
897  if (isBoundaryEdge[fEdges[fe]])
898  {
899  nFaceBoundaryEdges++;
900  }
901  }
902 
903  if (nFaceBoundaryEdges == fEdges.size())
904  {
905  if
906  (
907  checkCollapse
908  && !isCollapsedFace
909  (
910  newPoints,
911  neiCc,
912  minArea,
913  maxNonOrtho,
914  faceI
915  )
916  )
917  {
918  nPrevented++;
919  //Pout<< "Preventing collapse of coupled face "
920  // << faceI
921  // << " with all boundary edges "
922  // << " at " << mesh_.faceCentres()[faceI]
923  // << endl;
924  }
925  else
926  {
927  facePatch[faceI] = getBafflePatch(facePatch, faceI);
928  nBaffleFaces++;
929 
930  // Do NOT update boundary data since this would grow
931  // blocked faces across gaps.
932  }
933  }
934  }
935 
936  faceI++;
937  }
938  }
939  }
940 
941  Info<< "markFacesOnProblemCells : marked "
942  << returnReduce(nBaffleFaces, sumOp<label>())
943  << " additional internal faces to be converted into baffles."
944  << endl;
945 
946  if (checkCollapse)
947  {
948  Info<< "markFacesOnProblemCells : prevented "
949  << returnReduce(nPrevented, sumOp<label>())
950  << " internal faces fom getting converted into baffles."
951  << endl;
952  }
953 
954  return facePatch;
955 }
956 
957 
960 //Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
961 //(
962 // const dictionary& motionDict
963 //) const
964 //{
965 // // Construct addressing engine.
966 // autoPtr<indirectPrimitivePatch> ppPtr
967 // (
968 // meshRefinement::makePatch
969 // (
970 // mesh_,
971 // meshedPatches()
972 // )
973 // );
974 // const indirectPrimitivePatch& pp = ppPtr();
975 // const pointField& localPoints = pp.localPoints();
976 // const labelList& meshPoints = pp.meshPoints();
977 //
978 // // Find nearest (non-baffle) surface
979 // pointField newPoints(mesh_.points());
980 // {
981 // List<pointIndexHit> hitInfo;
982 // labelList hitSurface;
983 // surfaces_.findNearest
984 // (
985 // surfaces_.getUnnamedSurfaces(),
986 // localPoints,
987 // scalarField(localPoints.size(), sqr(GREAT)),// sqr of attraction
988 // hitSurface,
989 // hitInfo
990 // );
991 //
992 // forAll(hitInfo, i)
993 // {
994 // if (hitInfo[i].hit())
995 // {
996 // //label pointI = meshPoints[i];
997 // //Pout<< " " << pointI << " moved from "
998 // // << mesh_.points()[pointI] << " by "
999 // // << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
1000 // // << endl;
1001 // newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
1002 // }
1003 // }
1004 // }
1005 //
1006 // // Per face (internal or coupled!) the patch that the
1007 // // baffle should get (or -1).
1008 // labelList facePatch(mesh_.nFaces(), -1);
1009 // // Count of baffled faces
1010 // label nBaffleFaces = 0;
1011 //
1012 // {
1013 // pointField oldPoints(mesh_.points());
1014 // mesh_.movePoints(newPoints);
1015 // faceSet wrongFaces(mesh_, "wrongFaces", 100);
1016 // {
1017 // //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
1018 //
1019 // // Just check the errors from squashing
1020 // // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1021 //
1022 // const labelList allFaces(identity(mesh_.nFaces()));
1023 // label nWrongFaces = 0;
1024 //
1025 // scalar minArea(readScalar(motionDict.lookup("minArea")));
1026 // if (minArea > -SMALL)
1027 // {
1028 // polyMeshGeometry::checkFaceArea
1029 // (
1030 // false,
1031 // minArea,
1032 // mesh_,
1033 // mesh_.faceAreas(),
1034 // allFaces,
1035 // &wrongFaces
1036 // );
1037 //
1038 // label nNewWrongFaces = returnReduce
1039 // (
1040 // wrongFaces.size(),
1041 // sumOp<label>()
1042 // );
1043 //
1044 // Info<< " faces with area < "
1045 // << setw(5) << minArea
1046 // << " m^2 : "
1047 // << nNewWrongFaces-nWrongFaces << endl;
1048 //
1049 // nWrongFaces = nNewWrongFaces;
1050 // }
1051 //
1053 // scalar minDet = 0.01;
1054 // if (minDet > -1)
1055 // {
1056 // polyMeshGeometry::checkCellDeterminant
1057 // (
1058 // false,
1059 // minDet,
1060 // mesh_,
1061 // mesh_.faceAreas(),
1062 // allFaces,
1063 // polyMeshGeometry::affectedCells(mesh_, allFaces),
1064 // &wrongFaces
1065 // );
1066 //
1067 // label nNewWrongFaces = returnReduce
1068 // (
1069 // wrongFaces.size(),
1070 // sumOp<label>()
1071 // );
1072 //
1073 // Info<< " faces on cells with determinant < "
1074 // << setw(5) << minDet << " : "
1075 // << nNewWrongFaces-nWrongFaces << endl;
1076 //
1077 // nWrongFaces = nNewWrongFaces;
1078 // }
1079 // }
1080 //
1081 //
1082 // forAllConstIter(faceSet, wrongFaces, iter)
1083 // {
1084 // label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
1085 //
1086 // if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
1087 // {
1088 // facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
1089 // nBaffleFaces++;
1090 //
1091 // //Pout<< " " << iter.key()
1092 // // //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
1093 // // << " is destined for patch " << facePatch[iter.key()]
1094 // // << endl;
1095 // }
1096 // }
1097 // // Restore points.
1098 // mesh_.movePoints(oldPoints);
1099 // }
1100 //
1101 //
1102 // Info<< "markFacesOnProblemCellsGeometric : marked "
1103 // << returnReduce(nBaffleFaces, sumOp<label>())
1104 // << " additional internal and coupled faces"
1105 // << " to be converted into baffles." << endl;
1106 //
1107 // syncTools::syncFaceList
1108 // (
1109 // mesh_,
1110 // facePatch,
1111 // maxEqOp<label>(),
1112 // false // no separation
1113 // );
1114 //
1115 // return facePatch;
1116 //}
1117 
1118 
1119 // ************************ vim: set sw=4 sts=4 et: ************************ //