From 421a14488ea9c7488602ece1b1afd3cdcff7a683 Mon Sep 17 00:00:00 2001
From: mattijs <mattijs@hunt.opencfd.co.uk>
Date: Tue, 6 Jan 2009 20:54:08 +0000
Subject: [PATCH] problem cell removal
---
src/autoMesh/Make/files | 1 +
.../autoHexMeshDriver/autoRefineDriver.C | 10 +-
.../meshRefinement/meshRefinement.H | 29 +-
.../meshRefinement/meshRefinementBaffles.C | 664 +-----------
.../meshRefinementProblemCells.C | 952 ++++++++++++++++++
.../refinementSurfaces/refinementSurfaces.C | 44 +-
.../refinementSurfaces/refinementSurfaces.H | 9 +
7 files changed, 1053 insertions(+), 656 deletions(-)
create mode 100644 src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
diff --git a/src/autoMesh/Make/files b/src/autoMesh/Make/files
index 96141445bea..fecec5db5db 100644
--- a/src/autoMesh/Make/files
+++ b/src/autoMesh/Make/files
@@ -15,6 +15,7 @@ $(autoHexMeshDriver)/pointData/pointData.C
$(autoHexMesh)/meshRefinement/meshRefinementBaffles.C
$(autoHexMesh)/meshRefinement/meshRefinement.C
$(autoHexMesh)/meshRefinement/meshRefinementMerge.C
+$(autoHexMesh)/meshRefinement/meshRefinementProblemCells.C
$(autoHexMesh)/meshRefinement/meshRefinementRefine.C
$(autoHexMesh)/refinementSurfaces/refinementSurfaces.C
$(autoHexMesh)/shellSurfaces/shellSurfaces.C
diff --git a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
index ed9fd058086..fc70f98ccc6 100644
--- a/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
+++ b/src/autoMesh/autoHexMesh/autoHexMeshDriver/autoRefineDriver.C
@@ -518,7 +518,9 @@ void Foam::autoRefineDriver::baffleAndSplitMesh
// be like boundary face from now on so not coupled anymore.
meshRefiner_.baffleAndSplitMesh
(
- handleSnapProblems,
+ handleSnapProblems, // detect&remove potential snap problem
+ false, // perpendicular edge connected cells
+ scalarField(0), // per region perpendicular angle
!handleSnapProblems, // merge free standing baffles?
motionDict,
const_cast<Time&>(mesh.time()),
@@ -592,10 +594,14 @@ void Foam::autoRefineDriver::splitAndMergeBaffles
const_cast<Time&>(mesh.time())++;
}
+ const scalarField& perpAngle = meshRefiner_.surfaces().perpendicularAngle();
+
meshRefiner_.baffleAndSplitMesh
(
handleSnapProblems,
- false, // merge free standing baffles?
+ handleSnapProblems, // remove perp edge connected cells
+ perpAngle, // perp angle
+ false, // merge free standing baffles?
motionDict,
const_cast<Time&>(mesh.time()),
globalToPatch_,
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index da533c4d136..5895f5ae3ad 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -36,6 +36,7 @@ SourceFiles
meshRefinement.C
meshRefinementBaffles.C
meshRefinementMerge.C
+ meshRefinementProblemCells.C
meshRefinementRefine.C
\*---------------------------------------------------------------------------*/
@@ -51,6 +52,7 @@ SourceFiles
#include "indirectPrimitivePatch.H"
#include "pointFieldsFwd.H"
#include "Tuple2.H"
+#include "pointIndexHit.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -345,6 +347,8 @@ private:
polyTopoChange& meshMod
) const;
+ // Problem cell handling
+
//- Helper function to mark face as being on 'boundary'. Used by
// markFacesOnProblemCells
void markBoundaryFace
@@ -355,15 +359,32 @@ private:
boolList& isBoundaryPoint
) const;
+ void findNearest
+ (
+ const labelList& meshFaces,
+ List<pointIndexHit>& nearestInfo,
+ labelList& nearestSurface,
+ labelList& nearestRegion,
+ vectorField& nearestNormal
+ ) const;
+
+ Map<label> findEdgeConnectedProblemCells
+ (
+ const scalarField& perpendicularAngle,
+ const labelList&
+ ) const;
+
//- Returns list with for every internal face -1 or the patch
- // they should be baffled into.
+ // they should be baffled into. If removeEdgeConnectedCells is set
+ // removes cells based on perpendicularAngle.
labelList markFacesOnProblemCells
(
+ const bool removeEdgeConnectedCells,
+ const scalarField& perpendicularAngle,
const labelList& globalToPatch
) const;
- //- Returns list with for every internal face -1 or the patch
- // they should be baffled into.
+ //- Initial test of marking faces using geometric information.
labelList markFacesOnProblemCellsGeometric
(
const dictionary& motionDict,
@@ -589,6 +610,8 @@ public:
void baffleAndSplitMesh
(
const bool handleSnapProblems,
+ const bool removeEdgeConnectedCells,
+ const scalarField& perpendicularAngle,
const bool mergeFreeStanding,
const dictionary& motionDict,
Time& runTime,
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 950cc7251b4..6c5315a12a0 100644
--- a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -44,10 +44,6 @@ License
#include "OFstream.H"
#include "regionSplit.H"
#include "removeCells.H"
-#include "motionSmoother.H"
-#include "polyMeshGeometry.H"
-#include "IOmanip.H"
-#include "cellSet.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -230,7 +226,7 @@ void Foam::meshRefinement::getBafflePatches
label vertI = 0;
if (debug&OBJINTERSECTIONS)
{
- str.reset(new OFstream(mesh_.time().path()/"intersections.obj"));
+ str.reset(new OFstream(mesh_.time().timePath()/"intersections.obj"));
Pout<< "getBafflePatches : Writing surface intersections to file "
<< str().name() << nl << endl;
@@ -510,642 +506,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
}
-void Foam::meshRefinement::markBoundaryFace
-(
- const label faceI,
- boolList& isBoundaryFace,
- boolList& isBoundaryEdge,
- boolList& isBoundaryPoint
-) const
-{
- isBoundaryFace[faceI] = true;
-
- const labelList& fEdges = mesh_.faceEdges(faceI);
-
- forAll(fEdges, fp)
- {
- isBoundaryEdge[fEdges[fp]] = true;
- }
-
- const face& f = mesh_.faces()[faceI];
-
- forAll(f, fp)
- {
- isBoundaryPoint[f[fp]] = true;
- }
-}
-
-
-// Returns list with for every internal face -1 or the patch they should
-// be baffled into. Gets run after createBaffles so all the surface
-// intersections have already been turned into baffles. Used to remove cells
-// by baffling all their faces and have the splitMeshRegions chuck away non
-// used regions.
-Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
-(
- const labelList& globalToPatch
-) const
-{
- const labelList& cellLevel = meshCutter_.cellLevel();
- const labelList& pointLevel = meshCutter_.pointLevel();
- const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-
-
- // Per internal face (boundary faces not used) the patch that the
- // baffle should get (or -1)
- labelList facePatch(mesh_.nFaces(), -1);
-
- // Mark all points and edges on baffle patches (so not on any inlets,
- // outlets etc.)
- boolList isBoundaryPoint(mesh_.nPoints(), false);
- boolList isBoundaryEdge(mesh_.nEdges(), false);
- boolList isBoundaryFace(mesh_.nFaces(), false);
-
- // Fill boundary data. All elements on meshed patches get marked.
- // Get the labels of added patches.
- labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
-
- forAll(adaptPatchIDs, i)
- {
- label patchI = adaptPatchIDs[i];
-
- const polyPatch& pp = patches[patchI];
-
- label faceI = pp.start();
-
- forAll(pp, j)
- {
- markBoundaryFace
- (
- faceI,
- isBoundaryFace,
- isBoundaryEdge,
- isBoundaryPoint
- );
-
- faceI++;
- }
- }
-
- syncTools::syncPointList
- (
- mesh_,
- isBoundaryPoint,
- orEqOp<bool>(),
- false, // null value
- false // no separation
- );
-
- syncTools::syncEdgeList
- (
- mesh_,
- isBoundaryEdge,
- orEqOp<bool>(),
- false, // null value
- false // no separation
- );
-
- syncTools::syncFaceList
- (
- mesh_,
- isBoundaryFace,
- orEqOp<bool>(),
- false // no separation
- );
-
-
- // For each cell count the number of anchor points that are on
- // the boundary:
- // 8 : check the number of (baffle) boundary faces. If 3 or more block
- // off the cell since the cell would get squeezed down to a diamond
- // (probably; if the 3 or more faces are unrefined (only use the
- // anchor points))
- // 7 : store. Used to check later on whether there are points with
- // 3 or more of these cells. (note that on a flat surface a boundary
- // point will only have 4 cells connected to it)
-
- // Does cell have exactly 7 of its 8 anchor points on the boundary?
- PackedList<1> hasSevenBoundaryAnchorPoints(mesh_.nCells(), 0u);
- // If so what is the remaining non-boundary anchor point?
- labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
-
- // On-the-fly addressing storage.
- DynamicList<label> dynFEdges;
- DynamicList<label> dynCPoints;
-
- // Count of faces marked for baffling
- label nBaffleFaces = 0;
-
- forAll(cellLevel, cellI)
- {
- const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
-
- // Get number of anchor points (pointLevel == cellLevel)
-
- label nBoundaryAnchors = 0;
- label nNonAnchorBoundary = 0;
- label nonBoundaryAnchor = -1;
-
- forAll(cPoints, i)
- {
- label pointI = cPoints[i];
-
- if (pointLevel[pointI] <= cellLevel[cellI])
- {
- // Anchor point
- if (isBoundaryPoint[pointI])
- {
- nBoundaryAnchors++;
- }
- else
- {
- // Anchor point which is not on the surface
- nonBoundaryAnchor = pointI;
- }
- }
- else if (isBoundaryPoint[pointI])
- {
- nNonAnchorBoundary++;
- }
- }
-
- if (nBoundaryAnchors == 8)
- {
- const cell& cFaces = mesh_.cells()[cellI];
-
- // Count boundary faces.
- label nBfaces = 0;
-
- forAll(cFaces, cFaceI)
- {
- if (isBoundaryFace[cFaces[cFaceI]])
- {
- nBfaces++;
- }
- }
-
- // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
- // We assume that this situation is where there is a single
- // cell sticking out which would get flattened.
-
- // Eugene: delete cell no matter what.
- //if (nBfaces > 1)
- {
- forAll(cFaces, cf)
- {
- label faceI = cFaces[cf];
-
- if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
- {
- facePatch[faceI] = getBafflePatch(facePatch, faceI);
- nBaffleFaces++;
-
- // Mark face as a 'boundary'
- markBoundaryFace
- (
- faceI,
- isBoundaryFace,
- isBoundaryEdge,
- isBoundaryPoint
- );
- }
- }
- }
- }
- else if (nBoundaryAnchors == 7)
- {
- // Mark the cell. Store the (single!) non-boundary anchor point.
- hasSevenBoundaryAnchorPoints.set(cellI, 1u);
- nonBoundaryAnchors.insert(nonBoundaryAnchor);
- }
- }
-
-
- // Loop over all points. If a point is connected to 4 or more cells
- // with 7 anchor points on the boundary set those cell's non-boundary faces
- // to baffles
-
- DynamicList<label> dynPCells;
-
- forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
- {
- label pointI = iter.key();
-
- const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
-
- // Count number of 'hasSevenBoundaryAnchorPoints' cells.
- label n = 0;
-
- forAll(pCells, i)
- {
- if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
- {
- n++;
- }
- }
-
- if (n > 3)
- {
- // Point in danger of being what? Remove all 7-cells.
- forAll(pCells, i)
- {
- label cellI = pCells[i];
-
- if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
- {
- const cell& cFaces = mesh_.cells()[cellI];
-
- forAll(cFaces, cf)
- {
- label faceI = cFaces[cf];
-
- if
- (
- facePatch[faceI] == -1
- && mesh_.isInternalFace(faceI)
- )
- {
- facePatch[faceI] = getBafflePatch(facePatch, faceI);
- nBaffleFaces++;
-
- // Mark face as a 'boundary'
- markBoundaryFace
- (
- faceI,
- isBoundaryFace,
- isBoundaryEdge,
- isBoundaryPoint
- );
- }
- }
- }
- }
- }
- }
-
-
- // Sync all. (note that pointdata and facedata not used anymore but sync
- // anyway)
-
- syncTools::syncPointList
- (
- mesh_,
- isBoundaryPoint,
- orEqOp<bool>(),
- false, // null value
- false // no separation
- );
-
- syncTools::syncEdgeList
- (
- mesh_,
- isBoundaryEdge,
- orEqOp<bool>(),
- false, // null value
- false // no separation
- );
-
- syncTools::syncFaceList
- (
- mesh_,
- isBoundaryFace,
- orEqOp<bool>(),
- false // no separation
- );
-
-
- // Find faces with all edges on the boundary and make them baffles
- for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
- {
- if (facePatch[faceI] == -1)
- {
- const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
- label nFaceBoundaryEdges = 0;
-
- forAll(fEdges, fe)
- {
- if (isBoundaryEdge[fEdges[fe]])
- {
- nFaceBoundaryEdges++;
- }
- }
-
- if (nFaceBoundaryEdges == fEdges.size())
- {
- facePatch[faceI] = getBafflePatch(facePatch, faceI);
- nBaffleFaces++;
-
- // Do NOT update boundary data since this would grow blocked
- // faces across gaps.
- }
- }
- }
-
- forAll(patches, patchI)
- {
- const polyPatch& pp = patches[patchI];
-
- if (pp.coupled())
- {
- label faceI = pp.start();
-
- forAll(pp, i)
- {
- if (facePatch[faceI] == -1)
- {
- const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
- label nFaceBoundaryEdges = 0;
-
- forAll(fEdges, fe)
- {
- if (isBoundaryEdge[fEdges[fe]])
- {
- nFaceBoundaryEdges++;
- }
- }
-
- if (nFaceBoundaryEdges == fEdges.size())
- {
- facePatch[faceI] = getBafflePatch(facePatch, faceI);
- nBaffleFaces++;
-
- // Do NOT update boundary data since this would grow
- // blocked faces across gaps.
- }
- }
-
- faceI++;
- }
- }
- }
-
- Info<< "markFacesOnProblemCells : marked "
- << returnReduce(nBaffleFaces, sumOp<label>())
- << " additional internal faces to be converted into baffles."
- << endl;
-
- return facePatch;
-}
-//XXXXXXXXXXXXXX
-// Mark faces to be baffled to prevent snapping problems. Does
-// test to find nearest surface and checks which faces would get squashed.
-Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
-(
- const dictionary& motionDict,
- const labelList& globalToPatch
-) const
-{
- // Get the labels of added patches.
- labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
-
- // Construct addressing engine.
- autoPtr<indirectPrimitivePatch> ppPtr
- (
- meshRefinement::makePatch
- (
- mesh_,
- adaptPatchIDs
- )
- );
- const indirectPrimitivePatch& pp = ppPtr();
- const pointField& localPoints = pp.localPoints();
- const labelList& meshPoints = pp.meshPoints();
-
- // Find nearest (non-baffle) surface
- pointField newPoints(mesh_.points());
- {
- List<pointIndexHit> hitInfo;
- labelList hitSurface;
- surfaces_.findNearest
- (
- surfaces_.getUnnamedSurfaces(),
- localPoints,
- scalarField(localPoints.size(), sqr(GREAT)), // sqr of attraction
- hitSurface,
- hitInfo
- );
-
- forAll(hitInfo, i)
- {
- if (hitInfo[i].hit())
- {
- //label pointI = meshPoints[i];
- //Pout<< " " << pointI << " moved from "
- // << mesh_.points()[pointI] << " by "
- // << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
- // << endl;
- newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
- }
- }
- }
-
- // Per face (internal or coupled!) the patch that the
- // baffle should get (or -1).
- labelList facePatch(mesh_.nFaces(), -1);
- // Count of baffled faces
- label nBaffleFaces = 0;
-
-
-// // Sync position? Or not since same face on both side so just sync
-// // result of baffle.
-//
-// const scalar minArea(readScalar(motionDict.lookup("minArea")));
-//
-// Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
-// << minArea << endl;
-//
-// pointField facePoints;
-// for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
-// {
-// const face& f = mesh_.faces()[faceI];
-//
-// bool usesPatchPoint = false;
-//
-// facePoints.setSize(f.size());
-// forAll(f, fp)
-// {
-// Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
-//
-// if (iter != pp.meshPointMap().end())
-// {
-// facePoints[fp] = newPosition[iter()];
-// usesPatchPoint = true;
-// }
-// else
-// {
-// facePoints[fp] = mesh_.points()[f[fp]];
-// }
-// }
-//
-// if (usesPatchPoint)
-// {
-// // Check area of face wrt original area
-// face identFace(identity(f.size()));
-//
-// if (identFace.mag(facePoints) < minArea)
-// {
-// facePatch[faceI] = getBafflePatch(facePatch, faceI);
-// nBaffleFaces++;
-// }
-// }
-// }
-//
-//
-// const polyBoundaryMesh& patches = mesh_.boundaryMesh();
-// forAll(patches, patchI)
-// {
-// const polyPatch& pp = patches[patchI];
-//
-// if (pp.coupled())
-// {
-// forAll(pp, i)
-// {
-// label faceI = pp.start()+i;
-//
-// const face& f = mesh_.faces()[faceI];
-//
-// bool usesPatchPoint = false;
-//
-// facePoints.setSize(f.size());
-// forAll(f, fp)
-// {
-// Map<label>::const_iterator iter =
-// pp.meshPointMap().find(f[fp]);
-//
-// if (iter != pp.meshPointMap().end())
-// {
-// facePoints[fp] = newPosition[iter()];
-// usesPatchPoint = true;
-// }
-// else
-// {
-// facePoints[fp] = mesh_.points()[f[fp]];
-// }
-// }
-//
-// if (usesPatchPoint)
-// {
-// // Check area of face wrt original area
-// face identFace(identity(f.size()));
-//
-// if (identFace.mag(facePoints) < minArea)
-// {
-// facePatch[faceI] = getBafflePatch(facePatch, faceI);
-// nBaffleFaces++;
-// }
-// }
-// }
-// }
-// }
-
- {
- pointField oldPoints(mesh_.points());
- mesh_.movePoints(newPoints);
- faceSet wrongFaces(mesh_, "wrongFaces", 100);
- {
- //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
-
- // Just check the errors from squashing
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- const labelList allFaces(identity(mesh_.nFaces()));
- label nWrongFaces = 0;
-
- scalar minArea(readScalar(motionDict.lookup("minArea")));
- if (minArea > -SMALL)
- {
- polyMeshGeometry::checkFaceArea
- (
- false,
- minArea,
- mesh_,
- mesh_.faceAreas(),
- allFaces,
- &wrongFaces
- );
-
- label nNewWrongFaces = returnReduce
- (
- wrongFaces.size(),
- sumOp<label>()
- );
-
- Info<< " faces with area < "
- << setw(5) << minArea
- << " m^2 : "
- << nNewWrongFaces-nWrongFaces << endl;
-
- nWrongFaces = nNewWrongFaces;
- }
-
-// scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
- scalar minDet = 0.01;
- if (minDet > -1)
- {
- polyMeshGeometry::checkCellDeterminant
- (
- false,
- minDet,
- mesh_,
- mesh_.faceAreas(),
- allFaces,
- polyMeshGeometry::affectedCells(mesh_, allFaces),
- &wrongFaces
- );
-
- label nNewWrongFaces = returnReduce
- (
- wrongFaces.size(),
- sumOp<label>()
- );
-
- Info<< " faces on cells with determinant < "
- << setw(5) << minDet << " : "
- << nNewWrongFaces-nWrongFaces << endl;
-
- nWrongFaces = nNewWrongFaces;
- }
- }
-
-
- forAllConstIter(faceSet, wrongFaces, iter)
- {
- label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
-
- if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
- {
- facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
- nBaffleFaces++;
-
- //Pout<< " " << iter.key()
- // //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
- // << " is destined for patch " << facePatch[iter.key()]
- // << endl;
- }
- }
- // Restore points.
- mesh_.movePoints(oldPoints);
- }
-
-
- Info<< "markFacesOnProblemCellsGeometric : marked "
- << returnReduce(nBaffleFaces, sumOp<label>())
- << " additional internal and coupled faces"
- << " to be converted into baffles." << endl;
-
- syncTools::syncFaceList
- (
- mesh_,
- facePatch,
- maxEqOp<label>(),
- false // no separation
- );
-
- return facePatch;
-}
-//XXXXXXXX
-
-
// Return a list of coupled face pairs, i.e. faces that use the same vertices.
// (this information is recalculated instead of maintained since would be too
// hard across splitMeshRegions).
@@ -1855,13 +1215,13 @@ void Foam::meshRefinement::findCellZoneTopo
break;
}
- // Synchronise regionToCellZone.
- // Note:
- // - region numbers are identical on all processors
- // - keepRegion is identical ,,
- // - cellZones are identical ,,
- Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
- Pstream::listCombineScatter(regionToCellZone);
+ // Synchronise regionToCellZone.
+ // Note:
+ // - region numbers are identical on all processors
+ // - keepRegion is identical ,,
+ // - cellZones are identical ,,
+ Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
+ Pstream::listCombineScatter(regionToCellZone);
}
@@ -1904,6 +1264,8 @@ void Foam::meshRefinement::findCellZoneTopo
void Foam::meshRefinement::baffleAndSplitMesh
(
const bool handleSnapProblems,
+ const bool removeEdgeConnectedCells,
+ const scalarField& perpendicularAngle,
const bool mergeFreeStanding,
const dictionary& motionDict,
Time& runTime,
@@ -1981,6 +1343,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
(
markFacesOnProblemCells
(
+ removeEdgeConnectedCells,
+ perpendicularAngle,
globalToPatch
)
//markFacesOnProblemCellsGeometric
@@ -2024,6 +1388,8 @@ void Foam::meshRefinement::baffleAndSplitMesh
(
markFacesOnProblemCells
(
+ removeEdgeConnectedCells,
+ perpendicularAngle,
globalToPatch
)
);
@@ -2099,7 +1465,7 @@ void Foam::meshRefinement::baffleAndSplitMesh
{
Pout<< "Writing subsetted mesh to time " << mesh_.time().timeName()
<< endl;
- write(debug, runTime.path()/runTime.timeName());
+ write(debug, runTime.timePath());
Pout<< "Dumped debug data in = "
<< runTime.cpuTimeIncrement() << " s\n" << nl << endl;
}
diff --git a/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
new file mode 100644
index 00000000000..8a027fa654b
--- /dev/null
+++ b/src/autoMesh/autoHexMesh/meshRefinement/meshRefinementProblemCells.C
@@ -0,0 +1,952 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 1991-2008 OpenCFD Ltd.
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+License
+ This file is part of OpenFOAM.
+
+ OpenFOAM is free software; you can redistribute it and/or modify it
+ under the terms of the GNU General Public License as published by the
+ Free Software Foundation; either version 2 of the License, or (at your
+ option) any later version.
+
+ OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+ ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+ FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+ for more details.
+
+ You should have received a copy of the GNU General Public License
+ along with OpenFOAM; if not, write to the Free Software Foundation,
+ Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*----------------------------------------------------------------------------*/
+
+#include "meshRefinement.H"
+#include "fvMesh.H"
+#include "syncTools.H"
+#include "Time.H"
+#include "refinementSurfaces.H"
+#include "pointSet.H"
+#include "faceSet.H"
+#include "indirectPrimitivePatch.H"
+#include "OFstream.H"
+#include "cellSet.H"
+#include "searchableSurfaces.H"
+#include "polyMeshGeometry.H"
+#include "IOmanip.H"
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+void Foam::meshRefinement::markBoundaryFace
+(
+ const label faceI,
+ boolList& isBoundaryFace,
+ boolList& isBoundaryEdge,
+ boolList& isBoundaryPoint
+) const
+{
+ isBoundaryFace[faceI] = true;
+
+ const labelList& fEdges = mesh_.faceEdges(faceI);
+
+ forAll(fEdges, fp)
+ {
+ isBoundaryEdge[fEdges[fp]] = true;
+ }
+
+ const face& f = mesh_.faces()[faceI];
+
+ forAll(f, fp)
+ {
+ isBoundaryPoint[f[fp]] = true;
+ }
+}
+
+
+void Foam::meshRefinement::findNearest
+(
+ const labelList& meshFaces,
+ List<pointIndexHit>& nearestInfo,
+ labelList& nearestSurface,
+ labelList& nearestRegion,
+ vectorField& nearestNormal
+) const
+{
+ pointField fc(meshFaces.size());
+ forAll(meshFaces, i)
+ {
+ fc[i] = mesh_.faceCentres()[meshFaces[i]];
+ }
+
+ const labelList allSurfaces(identity(surfaces_.surfaces().size()));
+
+ surfaces_.findNearest
+ (
+ allSurfaces,
+ fc,
+ scalarField(fc.size(), sqr(GREAT)), // sqr of attraction
+ nearestSurface,
+ nearestInfo
+ );
+
+ // Do normal testing per surface.
+ nearestNormal.setSize(nearestInfo.size());
+ nearestRegion.setSize(nearestInfo.size());
+
+ forAll(allSurfaces, surfI)
+ {
+ DynamicList<pointIndexHit> localHits;
+
+ forAll(nearestSurface, i)
+ {
+ if (nearestSurface[i] == surfI)
+ {
+ localHits.append(nearestInfo[i]);
+ }
+ }
+
+ label geomI = surfaces_.surfaces()[surfI];
+
+ pointField localNormals;
+ surfaces_.geometry()[geomI].getNormal(localHits, localNormals);
+
+ labelList localRegion;
+ surfaces_.geometry()[geomI].getRegion(localHits, localRegion);
+
+ label localI = 0;
+ forAll(nearestSurface, i)
+ {
+ if (nearestSurface[i] == surfI)
+ {
+ nearestNormal[i] = localNormals[localI];
+ nearestRegion[i] = localRegion[localI];
+ localI++;
+ }
+ }
+ }
+}
+
+
+Foam::Map<Foam::label> Foam::meshRefinement::findEdgeConnectedProblemCells
+(
+ const scalarField& perpendicularAngle,
+ const labelList& globalToPatch
+) const
+{
+ labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+ // Construct addressing engine.
+ autoPtr<indirectPrimitivePatch> ppPtr
+ (
+ meshRefinement::makePatch
+ (
+ mesh_,
+ adaptPatchIDs
+ )
+ );
+ const indirectPrimitivePatch& pp = ppPtr();
+
+
+ // 1. Collect faces to test
+ // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+ DynamicList<label> candidateFaces(pp.size()/20);
+
+ const labelListList& edgeFaces = pp.edgeFaces();
+
+ const labelList& cellLevel = meshCutter_.cellLevel();
+
+ forAll(edgeFaces, edgeI)
+ {
+ const labelList& eFaces = edgeFaces[edgeI];
+
+ if (eFaces.size() == 2)
+ {
+ label face0 = pp.addressing()[eFaces[0]];
+ label face1 = pp.addressing()[eFaces[1]];
+
+ label cell0 = mesh_.faceOwner()[face0];
+ label cell1 = mesh_.faceOwner()[face1];
+
+ if (cellLevel[cell0] > cellLevel[cell1])
+ {
+ // cell0 smaller.
+ const vector& n0 = pp.faceNormals()[eFaces[0]];
+ const vector& n1 = pp.faceNormals()[eFaces[1]];
+
+ if (mag(n0 & n1) < 0.1)
+ {
+ candidateFaces.append(face0);
+ }
+ }
+ else if (cellLevel[cell1] > cellLevel[cell0])
+ {
+ // cell1 smaller.
+ const vector& n0 = pp.faceNormals()[eFaces[0]];
+ const vector& n1 = pp.faceNormals()[eFaces[1]];
+
+ if (mag(n0 & n1) < 0.1)
+ {
+ candidateFaces.append(face1);
+ }
+ }
+ }
+ }
+ candidateFaces.shrink();
+
+ Info<< "Testing " << returnReduce(candidateFaces.size(), sumOp<label>())
+ << " faces on edge-connected cells of differing level."
+ << endl;
+
+ if (debug)
+ {
+ faceSet fSet(mesh_, "edgeConnectedFaces", candidateFaces);
+ Pout<< "Writing " << fSet.size()
+ << " with problematic topology to faceSet "
+ << fSet.objectPath() << endl;
+ fSet.write();
+ }
+
+
+ // 2. Find nearest surface on candidate faces
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ List<pointIndexHit> nearestInfo;
+ labelList nearestSurface;
+ labelList nearestRegion;
+ vectorField nearestNormal;
+ findNearest
+ (
+ candidateFaces,
+ nearestInfo,
+ nearestSurface,
+ nearestRegion,
+ nearestNormal
+ );
+
+
+ // 3. Test angle to surface
+ // ~~~~~~~~~~~~~~~~~~~~~~~~
+
+ Map<label> candidateCells(candidateFaces.size());
+
+ faceSet perpFaces(mesh_, "perpendicularFaces", pp.size()/100);
+
+ forAll(candidateFaces, i)
+ {
+ label faceI = candidateFaces[i];
+
+ vector n = mesh_.faceAreas()[faceI];
+ n /= mag(n);
+
+ label region = surfaces_.globalRegion
+ (
+ nearestSurface[i],
+ nearestRegion[i]
+ );
+
+ scalar angle =
+ perpendicularAngle[region]
+ / 180.0
+ * mathematicalConstant::pi;
+
+ if (angle >= 0)
+ {
+ if (mag(n & nearestNormal[i]) < Foam::sin(angle))
+ {
+ perpFaces.insert(faceI);
+ candidateCells.insert
+ (
+ mesh_.faceOwner()[faceI],
+ globalToPatch[region]
+ );
+ }
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "Writing " << perpFaces.size()
+ << " faces that are perpendicular to the surface to set "
+ << perpFaces.objectPath() << endl;
+ perpFaces.write();
+ }
+ return candidateCells;
+}
+
+
+// Returns list with for every internal face -1 or the patch they should
+// be baffled into. Gets run after createBaffles so all the surface
+// intersections have already been turned into baffles. Used to remove cells
+// by baffling all their faces and have the splitMeshRegions chuck away non
+// used regions.
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCells
+(
+ const bool removeEdgeConnectedCells,
+ const scalarField& perpendicularAngle,
+ const labelList& globalToPatch
+) const
+{
+ const labelList& cellLevel = meshCutter_.cellLevel();
+ const labelList& pointLevel = meshCutter_.pointLevel();
+ const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+
+
+ // Per internal face (boundary faces not used) the patch that the
+ // baffle should get (or -1)
+ labelList facePatch(mesh_.nFaces(), -1);
+
+ // Mark all points and edges on baffle patches (so not on any inlets,
+ // outlets etc.)
+ boolList isBoundaryPoint(mesh_.nPoints(), false);
+ boolList isBoundaryEdge(mesh_.nEdges(), false);
+ boolList isBoundaryFace(mesh_.nFaces(), false);
+
+ // Fill boundary data. All elements on meshed patches get marked.
+ // Get the labels of added patches.
+ labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+ forAll(adaptPatchIDs, i)
+ {
+ label patchI = adaptPatchIDs[i];
+
+ const polyPatch& pp = patches[patchI];
+
+ label faceI = pp.start();
+
+ forAll(pp, j)
+ {
+ markBoundaryFace
+ (
+ faceI,
+ isBoundaryFace,
+ isBoundaryEdge,
+ isBoundaryPoint
+ );
+
+ faceI++;
+ }
+ }
+
+ // Count of faces marked for baffling
+ label nBaffleFaces = 0;
+
+ if (removeEdgeConnectedCells && max(perpendicularAngle) >= 0)
+ {
+ Info<< "markFacesOnProblemCells :"
+ << " Checking for edge-connected cells of highly differing sizes."
+ << endl;
+
+ // Pick up the cells that need to be removed and (a guess for)
+ // the patch they should be patched with.
+ Map<label> problemCells
+ (
+ findEdgeConnectedProblemCells
+ (
+ perpendicularAngle,
+ globalToPatch
+ )
+ );
+
+ // Baffle all faces of cells that need to be removed
+ forAllConstIter(Map<label>, problemCells, iter)
+ {
+ const cell& cFaces = mesh_.cells()[iter.key()];
+
+ forAll(cFaces, i)
+ {
+ label faceI = cFaces[i];
+
+ if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
+ {
+ facePatch[faceI] = getBafflePatch(facePatch, faceI);
+ nBaffleFaces++;
+
+ // Mark face as a 'boundary'
+ markBoundaryFace
+ (
+ faceI,
+ isBoundaryFace,
+ isBoundaryEdge,
+ isBoundaryPoint
+ );
+ }
+ }
+ }
+ Info<< "markFacesOnProblemCells : Marked "
+ << returnReduce(nBaffleFaces, sumOp<label>())
+ << " additional internal faces to be converted into baffles"
+ << " due to "
+ << returnReduce(problemCells.size(), sumOp<label>())
+ << " cells edge-connected to lower level cells." << endl;
+
+ if (debug)
+ {
+ cellSet problemCellSet(mesh_, "problemCells", problemCells.toc());
+ Pout<< "Writing " << problemCellSet.size()
+ << " cells that are edge connected to coarser cell to set "
+ << problemCellSet.objectPath() << endl;
+ problemCellSet.write();
+ }
+ }
+
+ syncTools::syncPointList
+ (
+ mesh_,
+ isBoundaryPoint,
+ orEqOp<bool>(),
+ false, // null value
+ false // no separation
+ );
+
+ syncTools::syncEdgeList
+ (
+ mesh_,
+ isBoundaryEdge,
+ orEqOp<bool>(),
+ false, // null value
+ false // no separation
+ );
+
+ syncTools::syncFaceList
+ (
+ mesh_,
+ isBoundaryFace,
+ orEqOp<bool>(),
+ false // no separation
+ );
+
+
+ // For each cell count the number of anchor points that are on
+ // the boundary:
+ // 8 : check the number of (baffle) boundary faces. If 3 or more block
+ // off the cell since the cell would get squeezed down to a diamond
+ // (probably; if the 3 or more faces are unrefined (only use the
+ // anchor points))
+ // 7 : store. Used to check later on whether there are points with
+ // 3 or more of these cells. (note that on a flat surface a boundary
+ // point will only have 4 cells connected to it)
+
+ // Does cell have exactly 7 of its 8 anchor points on the boundary?
+ PackedList<1> hasSevenBoundaryAnchorPoints(mesh_.nCells(), 0u);
+ // If so what is the remaining non-boundary anchor point?
+ labelHashSet nonBoundaryAnchors(mesh_.nCells()/10000);
+
+ // On-the-fly addressing storage.
+ DynamicList<label> dynFEdges;
+ DynamicList<label> dynCPoints;
+
+ forAll(cellLevel, cellI)
+ {
+ const labelList& cPoints = mesh_.cellPoints(cellI, dynCPoints);
+
+ // Get number of anchor points (pointLevel == cellLevel)
+
+ label nBoundaryAnchors = 0;
+ label nNonAnchorBoundary = 0;
+ label nonBoundaryAnchor = -1;
+
+ forAll(cPoints, i)
+ {
+ label pointI = cPoints[i];
+
+ if (pointLevel[pointI] <= cellLevel[cellI])
+ {
+ // Anchor point
+ if (isBoundaryPoint[pointI])
+ {
+ nBoundaryAnchors++;
+ }
+ else
+ {
+ // Anchor point which is not on the surface
+ nonBoundaryAnchor = pointI;
+ }
+ }
+ else if (isBoundaryPoint[pointI])
+ {
+ nNonAnchorBoundary++;
+ }
+ }
+
+ if (nBoundaryAnchors == 8)
+ {
+ const cell& cFaces = mesh_.cells()[cellI];
+
+ // Count boundary faces.
+ label nBfaces = 0;
+
+ forAll(cFaces, cFaceI)
+ {
+ if (isBoundaryFace[cFaces[cFaceI]])
+ {
+ nBfaces++;
+ }
+ }
+
+ // If nBfaces > 1 make all non-boundary non-baffle faces baffles.
+ // We assume that this situation is where there is a single
+ // cell sticking out which would get flattened.
+
+ // Eugene: delete cell no matter what.
+ //if (nBfaces > 1)
+ {
+ forAll(cFaces, cf)
+ {
+ label faceI = cFaces[cf];
+
+ if (facePatch[faceI] == -1 && mesh_.isInternalFace(faceI))
+ {
+ facePatch[faceI] = getBafflePatch(facePatch, faceI);
+ nBaffleFaces++;
+
+ // Mark face as a 'boundary'
+ markBoundaryFace
+ (
+ faceI,
+ isBoundaryFace,
+ isBoundaryEdge,
+ isBoundaryPoint
+ );
+ }
+ }
+ }
+ }
+ else if (nBoundaryAnchors == 7)
+ {
+ // Mark the cell. Store the (single!) non-boundary anchor point.
+ hasSevenBoundaryAnchorPoints.set(cellI, 1u);
+ nonBoundaryAnchors.insert(nonBoundaryAnchor);
+ }
+ }
+
+
+ // Loop over all points. If a point is connected to 4 or more cells
+ // with 7 anchor points on the boundary set those cell's non-boundary faces
+ // to baffles
+
+ DynamicList<label> dynPCells;
+
+ forAllConstIter(labelHashSet, nonBoundaryAnchors, iter)
+ {
+ label pointI = iter.key();
+
+ const labelList& pCells = mesh_.pointCells(pointI, dynPCells);
+
+ // Count number of 'hasSevenBoundaryAnchorPoints' cells.
+ label n = 0;
+
+ forAll(pCells, i)
+ {
+ if (hasSevenBoundaryAnchorPoints.get(pCells[i]) == 1u)
+ {
+ n++;
+ }
+ }
+
+ if (n > 3)
+ {
+ // Point in danger of being what? Remove all 7-cells.
+ forAll(pCells, i)
+ {
+ label cellI = pCells[i];
+
+ if (hasSevenBoundaryAnchorPoints.get(cellI) == 1u)
+ {
+ const cell& cFaces = mesh_.cells()[cellI];
+
+ forAll(cFaces, cf)
+ {
+ label faceI = cFaces[cf];
+
+ if
+ (
+ facePatch[faceI] == -1
+ && mesh_.isInternalFace(faceI)
+ )
+ {
+ facePatch[faceI] = getBafflePatch(facePatch, faceI);
+ nBaffleFaces++;
+
+ // Mark face as a 'boundary'
+ markBoundaryFace
+ (
+ faceI,
+ isBoundaryFace,
+ isBoundaryEdge,
+ isBoundaryPoint
+ );
+ }
+ }
+ }
+ }
+ }
+ }
+
+
+ // Sync all. (note that pointdata and facedata not used anymore but sync
+ // anyway)
+
+ syncTools::syncPointList
+ (
+ mesh_,
+ isBoundaryPoint,
+ orEqOp<bool>(),
+ false, // null value
+ false // no separation
+ );
+
+ syncTools::syncEdgeList
+ (
+ mesh_,
+ isBoundaryEdge,
+ orEqOp<bool>(),
+ false, // null value
+ false // no separation
+ );
+
+ syncTools::syncFaceList
+ (
+ mesh_,
+ isBoundaryFace,
+ orEqOp<bool>(),
+ false // no separation
+ );
+
+
+ // Find faces with all edges on the boundary and make them baffles
+ for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
+ {
+ if (facePatch[faceI] == -1)
+ {
+ const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
+ label nFaceBoundaryEdges = 0;
+
+ forAll(fEdges, fe)
+ {
+ if (isBoundaryEdge[fEdges[fe]])
+ {
+ nFaceBoundaryEdges++;
+ }
+ }
+
+ if (nFaceBoundaryEdges == fEdges.size())
+ {
+ facePatch[faceI] = getBafflePatch(facePatch, faceI);
+ nBaffleFaces++;
+
+ // Do NOT update boundary data since this would grow blocked
+ // faces across gaps.
+ }
+ }
+ }
+
+ forAll(patches, patchI)
+ {
+ const polyPatch& pp = patches[patchI];
+
+ if (pp.coupled())
+ {
+ label faceI = pp.start();
+
+ forAll(pp, i)
+ {
+ if (facePatch[faceI] == -1)
+ {
+ const labelList& fEdges = mesh_.faceEdges(faceI, dynFEdges);
+ label nFaceBoundaryEdges = 0;
+
+ forAll(fEdges, fe)
+ {
+ if (isBoundaryEdge[fEdges[fe]])
+ {
+ nFaceBoundaryEdges++;
+ }
+ }
+
+ if (nFaceBoundaryEdges == fEdges.size())
+ {
+ facePatch[faceI] = getBafflePatch(facePatch, faceI);
+ nBaffleFaces++;
+
+ // Do NOT update boundary data since this would grow
+ // blocked faces across gaps.
+ }
+ }
+
+ faceI++;
+ }
+ }
+ }
+
+ Info<< "markFacesOnProblemCells : marked "
+ << returnReduce(nBaffleFaces, sumOp<label>())
+ << " additional internal faces to be converted into baffles."
+ << endl;
+
+ return facePatch;
+}
+//XXXXXXXXXXXXXX
+// Mark faces to be baffled to prevent snapping problems. Does
+// test to find nearest surface and checks which faces would get squashed.
+Foam::labelList Foam::meshRefinement::markFacesOnProblemCellsGeometric
+(
+ const dictionary& motionDict,
+ const labelList& globalToPatch
+) const
+{
+ // Get the labels of added patches.
+ labelList adaptPatchIDs(meshRefinement::addedPatches(globalToPatch));
+
+ // Construct addressing engine.
+ autoPtr<indirectPrimitivePatch> ppPtr
+ (
+ meshRefinement::makePatch
+ (
+ mesh_,
+ adaptPatchIDs
+ )
+ );
+ const indirectPrimitivePatch& pp = ppPtr();
+ const pointField& localPoints = pp.localPoints();
+ const labelList& meshPoints = pp.meshPoints();
+
+ // Find nearest (non-baffle) surface
+ pointField newPoints(mesh_.points());
+ {
+ List<pointIndexHit> hitInfo;
+ labelList hitSurface;
+ surfaces_.findNearest
+ (
+ surfaces_.getUnnamedSurfaces(),
+ localPoints,
+ scalarField(localPoints.size(), sqr(GREAT)), // sqr of attraction
+ hitSurface,
+ hitInfo
+ );
+
+ forAll(hitInfo, i)
+ {
+ if (hitInfo[i].hit())
+ {
+ //label pointI = meshPoints[i];
+ //Pout<< " " << pointI << " moved from "
+ // << mesh_.points()[pointI] << " by "
+ // << mag(hitInfo[i].hitPoint()-mesh_.points()[pointI])
+ // << endl;
+ newPoints[meshPoints[i]] = hitInfo[i].hitPoint();
+ }
+ }
+ }
+
+ // Per face (internal or coupled!) the patch that the
+ // baffle should get (or -1).
+ labelList facePatch(mesh_.nFaces(), -1);
+ // Count of baffled faces
+ label nBaffleFaces = 0;
+
+
+// // Sync position? Or not since same face on both side so just sync
+// // result of baffle.
+//
+// const scalar minArea(readScalar(motionDict.lookup("minArea")));
+//
+// Pout<< "markFacesOnProblemCellsGeometric : Comparing to minArea:"
+// << minArea << endl;
+//
+// pointField facePoints;
+// for (label faceI = mesh_.nInternalFaces(); faceI < mesh_.nFaces(); faceI++)
+// {
+// const face& f = mesh_.faces()[faceI];
+//
+// bool usesPatchPoint = false;
+//
+// facePoints.setSize(f.size());
+// forAll(f, fp)
+// {
+// Map<label>::const_iterator iter = pp.meshPointMap().find(f[fp]);
+//
+// if (iter != pp.meshPointMap().end())
+// {
+// facePoints[fp] = newPosition[iter()];
+// usesPatchPoint = true;
+// }
+// else
+// {
+// facePoints[fp] = mesh_.points()[f[fp]];
+// }
+// }
+//
+// if (usesPatchPoint)
+// {
+// // Check area of face wrt original area
+// face identFace(identity(f.size()));
+//
+// if (identFace.mag(facePoints) < minArea)
+// {
+// facePatch[faceI] = getBafflePatch(facePatch, faceI);
+// nBaffleFaces++;
+// }
+// }
+// }
+//
+//
+// const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+// forAll(patches, patchI)
+// {
+// const polyPatch& pp = patches[patchI];
+//
+// if (pp.coupled())
+// {
+// forAll(pp, i)
+// {
+// label faceI = pp.start()+i;
+//
+// const face& f = mesh_.faces()[faceI];
+//
+// bool usesPatchPoint = false;
+//
+// facePoints.setSize(f.size());
+// forAll(f, fp)
+// {
+// Map<label>::const_iterator iter =
+// pp.meshPointMap().find(f[fp]);
+//
+// if (iter != pp.meshPointMap().end())
+// {
+// facePoints[fp] = newPosition[iter()];
+// usesPatchPoint = true;
+// }
+// else
+// {
+// facePoints[fp] = mesh_.points()[f[fp]];
+// }
+// }
+//
+// if (usesPatchPoint)
+// {
+// // Check area of face wrt original area
+// face identFace(identity(f.size()));
+//
+// if (identFace.mag(facePoints) < minArea)
+// {
+// facePatch[faceI] = getBafflePatch(facePatch, faceI);
+// nBaffleFaces++;
+// }
+// }
+// }
+// }
+// }
+
+ {
+ pointField oldPoints(mesh_.points());
+ mesh_.movePoints(newPoints);
+ faceSet wrongFaces(mesh_, "wrongFaces", 100);
+ {
+ //motionSmoother::checkMesh(false, mesh_, motionDict, wrongFaces);
+
+ // Just check the errors from squashing
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+ const labelList allFaces(identity(mesh_.nFaces()));
+ label nWrongFaces = 0;
+
+ scalar minArea(readScalar(motionDict.lookup("minArea")));
+ if (minArea > -SMALL)
+ {
+ polyMeshGeometry::checkFaceArea
+ (
+ false,
+ minArea,
+ mesh_,
+ mesh_.faceAreas(),
+ allFaces,
+ &wrongFaces
+ );
+
+ label nNewWrongFaces = returnReduce
+ (
+ wrongFaces.size(),
+ sumOp<label>()
+ );
+
+ Info<< " faces with area < "
+ << setw(5) << minArea
+ << " m^2 : "
+ << nNewWrongFaces-nWrongFaces << endl;
+
+ nWrongFaces = nNewWrongFaces;
+ }
+
+// scalar minDet(readScalar(motionDict.lookup("minDeterminant")));
+ scalar minDet = 0.01;
+ if (minDet > -1)
+ {
+ polyMeshGeometry::checkCellDeterminant
+ (
+ false,
+ minDet,
+ mesh_,
+ mesh_.faceAreas(),
+ allFaces,
+ polyMeshGeometry::affectedCells(mesh_, allFaces),
+ &wrongFaces
+ );
+
+ label nNewWrongFaces = returnReduce
+ (
+ wrongFaces.size(),
+ sumOp<label>()
+ );
+
+ Info<< " faces on cells with determinant < "
+ << setw(5) << minDet << " : "
+ << nNewWrongFaces-nWrongFaces << endl;
+
+ nWrongFaces = nNewWrongFaces;
+ }
+ }
+
+
+ forAllConstIter(faceSet, wrongFaces, iter)
+ {
+ label patchI = mesh_.boundaryMesh().whichPatch(iter.key());
+
+ if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
+ {
+ facePatch[iter.key()] = getBafflePatch(facePatch, iter.key());
+ nBaffleFaces++;
+
+ //Pout<< " " << iter.key()
+ // //<< " on patch " << mesh_.boundaryMesh()[patchI].name()
+ // << " is destined for patch " << facePatch[iter.key()]
+ // << endl;
+ }
+ }
+ // Restore points.
+ mesh_.movePoints(oldPoints);
+ }
+
+
+ Info<< "markFacesOnProblemCellsGeometric : marked "
+ << returnReduce(nBaffleFaces, sumOp<label>())
+ << " additional internal and coupled faces"
+ << " to be converted into baffles." << endl;
+
+ syncTools::syncFaceList
+ (
+ mesh_,
+ facePatch,
+ maxEqOp<label>(),
+ false // no separation
+ );
+
+ return facePatch;
+}
+//XXXXXXXX
+
+
+// ************************************************************************* //
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
index 2f23e9706cd..276aa2fd014 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.C
@@ -53,8 +53,10 @@ Foam::refinementSurfaces::refinementSurfaces
{
labelList globalMinLevel(surfaceDicts.size(), 0);
labelList globalMaxLevel(surfaceDicts.size(), 0);
+ scalarField globalAngle(surfaceDicts.size(), -GREAT);
List<Map<label> > regionMinLevel(surfaceDicts.size());
List<Map<label> > regionMaxLevel(surfaceDicts.size());
+ List<Map<scalar> > regionAngle(surfaceDicts.size());
//wordList globalPatchType(surfaceDicts.size());
//List<HashTable<word> > regionPatchType(surfaceDicts.size());
@@ -80,6 +82,12 @@ Foam::refinementSurfaces::refinementSurfaces
dict.lookup("zoneInside") >> zoneInside_[surfI];
}
+ // Global perpendicular angle
+ if (dict.found("perpendicularAngle"))
+ {
+ globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
+ }
+
//// Global patch name per surface
//if (dict.found("patchType"))
//{
@@ -130,6 +138,15 @@ Foam::refinementSurfaces::refinementSurfaces
<< exit(FatalError);
}
regionMaxLevel[surfI].insert(regionI, max);
+
+ if (regionDict.found("perpendicularAngle"))
+ {
+ regionAngle[surfI].insert
+ (
+ regionI,
+ readScalar(regionDict.lookup("perpendicularAngle"))
+ );
+ }
}
}
}
@@ -170,6 +187,8 @@ Foam::refinementSurfaces::refinementSurfaces
minLevel_ = 0;
maxLevel_.setSize(nRegions);
maxLevel_ = 0;
+ perpendicularAngle_.setSize(nRegions);
+ perpendicularAngle_ = -GREAT;
//patchName_.setSize(nRegions);
//patchType_.setSize(nRegions);
@@ -182,6 +201,7 @@ Foam::refinementSurfaces::refinementSurfaces
{
minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+ perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
}
// Overwrite with region specific information
@@ -191,6 +211,7 @@ Foam::refinementSurfaces::refinementSurfaces
minLevel_[globalRegionI] = iter();
maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
+ perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
// Check validity
if
@@ -240,8 +261,10 @@ Foam::refinementSurfaces::refinementSurfaces
{
labelList globalMinLevel(surfacesDict.size(), 0);
labelList globalMaxLevel(surfacesDict.size(), 0);
+ scalarField globalAngle(surfacesDict.size(), -GREAT);
List<Map<label> > regionMinLevel(surfacesDict.size());
List<Map<label> > regionMaxLevel(surfacesDict.size());
+ List<Map<scalar> > regionAngle(surfacesDict.size());
label surfI = 0;
forAllConstIter(dictionary, surfacesDict, iter)
@@ -274,6 +297,12 @@ Foam::refinementSurfaces::refinementSurfaces
dict.lookup("zoneInside") >> zoneInside_[surfI];
}
+ // Global perpendicular angle
+ if (dict.found("perpendicularAngle"))
+ {
+ globalAngle[surfI] = readScalar(dict.lookup("perpendicularAngle"));
+ }
+
if (dict.found("regions"))
{
const dictionary& regionsDict = dict.subDict("regions");
@@ -306,6 +335,15 @@ Foam::refinementSurfaces::refinementSurfaces
regionMinLevel[surfI].insert(regionI, refLevel[0]);
regionMaxLevel[surfI].insert(regionI, refLevel[1]);
+
+ if (regionDict.found("perpendicularAngle"))
+ {
+ regionAngle[surfI].insert
+ (
+ regionI,
+ readScalar(regionDict.lookup("perpendicularAngle"))
+ );
+ }
}
}
}
@@ -326,6 +364,8 @@ Foam::refinementSurfaces::refinementSurfaces
minLevel_ = 0;
maxLevel_.setSize(nRegions);
maxLevel_ = 0;
+ perpendicularAngle_.setSize(nRegions);
+ perpendicularAngle_ = -GREAT;
forAll(globalMinLevel, surfI)
@@ -337,6 +377,7 @@ Foam::refinementSurfaces::refinementSurfaces
{
minLevel_[regionOffset_[surfI] + i] = globalMinLevel[surfI];
maxLevel_[regionOffset_[surfI] + i] = globalMaxLevel[surfI];
+ perpendicularAngle_[regionOffset_[surfI] + i] = globalAngle[surfI];
}
// Overwrite with region specific information
@@ -346,6 +387,7 @@ Foam::refinementSurfaces::refinementSurfaces
minLevel_[globalRegionI] = iter();
maxLevel_[globalRegionI] = regionMaxLevel[surfI][iter.key()];
+ perpendicularAngle_[globalRegionI] = regionAngle[surfI][iter.key()];
// Check validity
if
@@ -454,8 +496,6 @@ void Foam::refinementSurfaces::setMinLevelFields
const shellSurfaces& shells
)
{
- //minLevelFields_.setSize(surfaces_.size());
-
forAll(surfaces_, surfI)
{
const searchableSurface& geom = allGeometry_[surfaces_[surfI]];
diff --git a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
index a20df19f29f..c727303173c 100644
--- a/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
+++ b/src/autoMesh/autoHexMesh/refinementSurfaces/refinementSurfaces.H
@@ -90,6 +90,9 @@ class refinementSurfaces
//- From global region number to refinement level
labelList maxLevel_;
+ //- From global region number to perpendicular angle
+ scalarField perpendicularAngle_;
+
// Private Member Functions
@@ -178,6 +181,12 @@ public:
return maxLevel_;
}
+ //- From global region number to perpendicular angle
+ const scalarField& perpendicularAngle() const
+ {
+ return perpendicularAngle_;
+ }
+
// Helper
--
GitLab