diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
index e8c2915f1388ac6adb33308df8383e528cec41a3..beba23a4da5f01e3bcfa84cf953ced9e9d43a3da 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMesh.C
@@ -58,6 +58,7 @@ Description
#include "globalIndex.H"
#include "IOmanip.H"
#include "decompositionModel.H"
+#include "fvMeshTools.H"
using namespace Foam;
@@ -813,6 +814,7 @@ int main(int argc, char *argv[])
readScalar(meshDict.lookup("mergeTolerance"))
);
+ const Switch keepPatches(meshDict.lookupOrDefault("keepPatches", false));
// Read decomposePar dictionary
@@ -1517,6 +1519,12 @@ int main(int argc, char *argv[])
motionDict
);
+ // Remove zero sized patches originating from faceZones
+ if (!keepPatches && !wantSnap && !wantLayers)
+ {
+ fvMeshTools::removeEmptyPatches(mesh, true);
+ }
+
writeMesh
(
"Refined mesh",
@@ -1559,6 +1567,12 @@ int main(int argc, char *argv[])
snapParams
);
+ // Remove zero sized patches originating from faceZones
+ if (!keepPatches && !wantLayers)
+ {
+ fvMeshTools::removeEmptyPatches(mesh, true);
+ }
+
writeMesh
(
"Snapped mesh",
@@ -1609,6 +1623,12 @@ int main(int argc, char *argv[])
distributor
);
+ // Remove zero sized patches originating from faceZones
+ if (!keepPatches)
+ {
+ fvMeshTools::removeEmptyPatches(mesh, true);
+ }
+
writeMesh
(
"Layer mesh",
@@ -1622,6 +1642,7 @@ int main(int argc, char *argv[])
}
+
{
// Check final mesh
Info<< "Checking final mesh ..." << endl;
diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
index 5a8508d151ec16293278ecca79cabe0339d42515..506a68786b9c4a4b55d7dff724397196dffca076 100644
--- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
+++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict
@@ -296,8 +296,18 @@ castellatedMeshControls
locationInMesh (5 0.28 0.43);
// Whether any faceZones (as specified in the refinementSurfaces)
- // are only on the boundary of corresponding cellZones or also allow
- // free-standing zone faces. Not used if there are no faceZones.
+ // are only on the boundary of corresponding cellZones.
+ // Not used if there are no faceZones. The behaviour has changed
+ // with respect to previous versions:
+ // true : all intersections with surface are put in faceZone
+ // (same behaviour as before)
+ // false : depending on the type of surface intersected:
+ // - if intersecting surface has faceZone only (so no cellZone)
+ // leave in faceZone (so behave as if set to true) (= changed
+ // behaviour)
+ // - if intersecting surface has faceZone and cellZone
+ // remove if inbetween same cellZone or if on boundary
+ // (same behaviour as before)
allowFreeStandingZoneFaces true;
diff --git a/src/dynamicMesh/fvMeshTools/fvMeshTools.C b/src/dynamicMesh/fvMeshTools/fvMeshTools.C
index ecbe00325d73bcf6032465a864b8f973195d4b3d..df8aaccdb4b027bd179017e8ec72207606cfd8b8 100644
--- a/src/dynamicMesh/fvMeshTools/fvMeshTools.C
+++ b/src/dynamicMesh/fvMeshTools/fvMeshTools.C
@@ -355,6 +355,69 @@ void Foam::fvMeshTools::reorderPatches
}
+Foam::labelList Foam::fvMeshTools::removeEmptyPatches
+(
+ fvMesh& mesh,
+ const bool validBoundary
+)
+{
+ const polyBoundaryMesh& pbm = mesh.boundaryMesh();
+
+ labelList newToOld(pbm.size());
+ labelList oldToNew(pbm.size(), -1);
+ label newI = 0;
+
+
+ // Assumes all non-coupled boundaries are on all processors!
+ forAll(pbm, patchI)
+ {
+ const polyPatch& pp = pbm[patchI];
+
+ if (!isA<processorPolyPatch>(pp))
+ {
+ label nFaces = pp.size();
+ if (validBoundary)
+ {
+ reduce(nFaces, sumOp<label>());
+ }
+
+ if (nFaces > 0)
+ {
+ newToOld[newI] = patchI;
+ oldToNew[patchI] = newI++;
+ }
+ }
+ }
+
+ // Same for processor patches (but need no reduction)
+ forAll(pbm, patchI)
+ {
+ const polyPatch& pp = pbm[patchI];
+
+ if (isA<processorPolyPatch>(pp) && pp.size())
+ {
+ newToOld[newI] = patchI;
+ oldToNew[patchI] = newI++;
+ }
+ }
+
+ newToOld.setSize(newI);
+
+ // Move all deleteable patches to the end
+ forAll(oldToNew, patchI)
+ {
+ if (oldToNew[patchI] == -1)
+ {
+ oldToNew[patchI] = newI++;
+ }
+ }
+
+ reorderPatches(mesh, oldToNew, newToOld.size(), validBoundary);
+
+ return newToOld;
+}
+
+
Foam::autoPtr<Foam::fvMesh> Foam::fvMeshTools::newMesh
(
const IOobject& io,
diff --git a/src/dynamicMesh/fvMeshTools/fvMeshTools.H b/src/dynamicMesh/fvMeshTools/fvMeshTools.H
index e458425d0b2870db96ece3e047c1b27a1e4a871f..37b368961c52ff790a4615c05e070cbf14465110 100644
--- a/src/dynamicMesh/fvMeshTools/fvMeshTools.H
+++ b/src/dynamicMesh/fvMeshTools/fvMeshTools.H
@@ -113,7 +113,7 @@ public:
static void zeroPatchFields(fvMesh& mesh, const label patchI);
//- Reorder and remove trailing patches. If validBoundary call is parallel
- // synced and all add the same patch with same settings
+ // synced
static void reorderPatches
(
fvMesh&,
@@ -122,6 +122,11 @@ public:
const bool validBoundary
);
+ //- Remove zero sized patches. All but processor patches are
+ // assumed to be present on all processors (so size will be reduced
+ // if validBoundary). Return map from new
+ // to old patches
+ static labelList removeEmptyPatches(fvMesh&, const bool validBoundary);
//- Read mesh or create dummy mesh (0 cells, >0 patches). Works in two
// modes according to masterOnlyReading:
@@ -133,7 +138,6 @@ public:
const IOobject& io,
const bool masterOnlyReading
);
-
};
diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
index f5bfcc5f7524e887ea04fe855481dce91a56614d..0a93e1e4a8a69ce9f235dd7fb0bd726007335465 100644
--- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
+++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/refinementParameters/refinementParameters.C
@@ -74,14 +74,7 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
if (dict.readIfPresent("locationInMesh", locationInMesh))
{
locationsInMesh_.append(locationInMesh);
- zonesInMesh_.append("noneIfNotSet");// special name for no cellZone
-
- if (dict.found("locationsInMesh"))
- {
- FatalIOErrorInFunction(dict)
- << "Cannot both specify 'locationInMesh' and 'locationsInMesh'"
- << exit(FatalIOError);
- }
+ zonesInMesh_.append("none"); // special name for no cellZone
}
List<Tuple2<point, word> > pointsToZone;
@@ -95,6 +88,10 @@ Foam::refinementParameters::refinementParameters(const dictionary& dict)
{
locationsInMesh_[nZones] = pointsToZone[i].first();
zonesInMesh_[nZones] = pointsToZone[i].second();
+ if (zonesInMesh_[nZones] == word::null)
+ {
+ zonesInMesh_[nZones] = "none";
+ }
nZones++;
}
}
@@ -152,12 +149,7 @@ Foam::labelList Foam::refinementParameters::addCellZonesToMesh
labelList zoneIDs(zonesInMesh_.size(), -1);
forAll(zonesInMesh_, i)
{
- if
- (
- zonesInMesh_[i] != word::null
- && zonesInMesh_[i] != "none"
- && zonesInMesh_[i] != "noneIfNotSet"
- )
+ if (zonesInMesh_[i] != word::null && zonesInMesh_[i] != "none")
{
zoneIDs[i] = surfaceZonesInfo::addCellZone
(
@@ -242,8 +234,8 @@ Foam::labelList Foam::refinementParameters::zonedLocations
{
if
(
- zonesInMesh[i] == word::null
- || zonesInMesh[i] != "noneIfNotSet"
+ zonesInMesh[i] != word::null
+ && zonesInMesh[i] != "none"
)
{
indices.append(i);
@@ -264,8 +256,8 @@ Foam::labelList Foam::refinementParameters::unzonedLocations
{
if
(
- zonesInMesh[i] != word::null
- && zonesInMesh[i] == "noneIfNotSet"
+ zonesInMesh[i] == word::null
+ || zonesInMesh[i] == "none"
)
{
indices.append(i);
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
index 706a130b1539ce0f8e4fd4f4dd9a14f37565dfd5..9b375b1b89f4364f150f9327b344ab3646f5a65d 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinement.H
@@ -507,7 +507,6 @@ private:
void getIntersections
(
const labelList& surfacesToTest,
- const labelList& neiLevel,
const pointField& neiCc,
const labelList& testFaces,
@@ -515,7 +514,19 @@ private:
labelList& globalRegion2
) const;
- //- Determine patches for baffles on all intersected unnamed faces
+ //- Calculate intersections on zoned faces. Return per face -1
+ // or the index of the surface and the orientation w.r.t. surface
+ void getIntersections
+ (
+ const labelList& surfacesToTest,
+ const pointField& neiCc,
+ const labelList& testFaces,
+
+ labelList& namedSurfaceIndex,
+ PackedBoolList& posOrientation
+ ) const;
+
+ //- Determine patches for baffles
void getBafflePatches
(
const labelList& globalToMasterPatch,
@@ -635,13 +646,13 @@ private:
labelList& cellToZone
) const;
- //- Finds zone per cell for cells inside named surfaces that have
- // an inside point specified.
+ //- Finds zone per cell for cells inside region for which name
+ // is specified.
void findCellZoneInsideWalk
(
- const labelList& locationSurfaces,
- const labelList& namedSurfaceIndex,
- const labelList& surfaceToCellZone,
+ const pointField& locationsInMesh,
+ const labelList& zonesInMesh,
+ const labelList& blockedFace, // per face -1 or some index >= 0
labelList& cellToZone
) const;
@@ -650,8 +661,8 @@ private:
void findCellZoneInsideWalk
(
const pointField& locationsInMesh,
- const wordList& regionsInMesh,
- const labelList& blockedFace, // per face -1 or some index >= 0
+ const wordList& zoneNamesInMesh,
+ const labelList& faceToZone, // per face -1 or some index >= 0
labelList& cellToZone
) const;
@@ -670,6 +681,7 @@ private:
void findCellZoneTopo
(
const pointField& locationsInMesh,
+ const labelList& allSurfaceIndex,
const labelList& namedSurfaceIndex,
const labelList& surfaceToCellZone,
labelList& cellToZone
@@ -679,10 +691,23 @@ private:
// - clear out any blocked faces inbetween same cell zone.
void makeConsistentFaceIndex
(
+ const labelList& zoneToNamedSurface,
const labelList& cellToZone,
labelList& namedSurfaceIndex
) const;
+ //- Calculate cellZone allocation
+ void zonify
+ (
+ const bool allowFreeStandingZoneFaces,
+ const pointField& locationsInMesh,
+ const wordList& zonesInMesh,
+
+ labelList& cellToZone,
+ labelList& namedSurfaceIndex,
+ PackedBoolList& posOrientation
+ ) const;
+
//- Put cells into cellZone, faces into faceZone
void zonify
(
diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
index 6fe484ac27c877787b19feeb9e348a4c8a41c4c2..7d350c4e683b2b16da68236192af67710c7be63d 100644
--- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
+++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C
@@ -51,6 +51,9 @@ License
#include "IOmanip.H"
#include "refinementParameters.H"
+#include "zeroGradientFvPatchFields.H"
+#include "volFields.H"
+
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::meshRefinement::createBaffle
@@ -169,7 +172,6 @@ Foam::label Foam::meshRefinement::createBaffle
void Foam::meshRefinement::getIntersections
(
const labelList& surfacesToTest,
- const labelList& neiLevel,
const pointField& neiCc,
const labelList& testFaces,
@@ -211,7 +213,7 @@ void Foam::meshRefinement::getIntersections
calcCellCellRays
(
neiCc,
- neiLevel,
+ labelList(neiCc.size(), -1),
testFaces,
start,
end,
@@ -222,6 +224,7 @@ void Foam::meshRefinement::getIntersections
// Do test for intersections
// ~~~~~~~~~~~~~~~~~~~~~~~~~
+
labelList surface1;
List<pointIndexHit> hit1;
labelList region1;
@@ -290,136 +293,133 @@ void Foam::meshRefinement::getBafflePatches
labelList& neiPatch
) const
{
- // Check all unnamed surfaces
+ // 1. Determine cell zones
+ // ~~~~~~~~~~~~~~~~~~~~~~~
+ // Note that this does not determine the surface that was intersected
+ // so that is done in step 2 below.
+ labelList cellToZone;
{
- const labelList testFaces(intersectedFaces());
-
- labelList globalRegion1;
- labelList globalRegion2;
- getIntersections
+ labelList namedSurfaceIndex;
+ PackedBoolList posOrientation;
+ zonify
(
- surfaceZonesInfo::getUnnamedSurfaces(surfaces_.surfZones()),
- neiLevel,
- neiCc,
- testFaces,
- globalRegion1,
- globalRegion2
+ true, // allowFreeStandingZoneFaces
+ locationsInMesh,
+ zonesInMesh,
+
+ cellToZone,
+ namedSurfaceIndex,
+ posOrientation
);
- ownPatch.setSize(mesh_.nFaces());
- ownPatch = -1;
- neiPatch.setSize(mesh_.nFaces());
- neiPatch = -1;
- forAll(testFaces, i)
+
+ // Some stats
+ if (debug)
{
- label faceI = testFaces[i];
- if (globalRegion1[faceI] != -1)
+ label nZones = gMax(cellToZone)+1;
+
+ label nUnvisited = 0;
+ label nBackgroundCells = 0;
+ labelList nZoneCells(nZones, 0);
+ forAll(cellToZone, cellI)
{
- ownPatch[faceI] = globalToMasterPatch[globalRegion1[faceI]];
+ label zoneI = cellToZone[cellI];
+ if (zoneI >= 0)
+ {
+ nZoneCells[zoneI]++;
+ }
+ else if (zoneI == -1)
+ {
+ nBackgroundCells++;
+ }
+ else if (zoneI == -2)
+ {
+ nUnvisited++;
+ }
+ else
+ {
+ FatalErrorIn("meshRefinement::getBafflePatches()")
+ << "problem" << exit(FatalError);
+ }
}
- if (globalRegion2[faceI] != -1)
+ reduce(nUnvisited, sumOp<label>());
+ reduce(nBackgroundCells, sumOp<label>());
+ forAll(nZoneCells, zoneI)
{
- neiPatch[faceI] = globalToMasterPatch[globalRegion2[faceI]];
+ reduce(nZoneCells[zoneI], sumOp<label>());
}
+ Info<< "nUnvisited :" << nUnvisited << endl;
+ Info<< "nBackgroundCells:" << nBackgroundCells << endl;
+ Info<< "nZoneCells :" << nZoneCells << endl;
}
}
- if (locationsInMesh.size() > 1)
- {
- // Now we need to filter out any possible intersections between
- // multiple regions. Only the faces on the outside of all
- // regions are candidates for baffling.
- // For this we need to go to per-cell information
-
- labelList cellToZone(mesh_.nCells(), -2);
-
-
- // Closed surfaces with cellZone specified.
- const labelList closedNamedSurfaces
- (
- surfaceZonesInfo::getClosedNamedSurfaces
- (
- surfaces_.surfZones(),
- surfaces_.geometry(),
- surfaces_.surfaces()
- )
- );
-
- if (closedNamedSurfaces.size())
- {
- Info<< "Found " << closedNamedSurfaces.size()
- << " closed, named surfaces. Assigning cells in/outside"
- << " these surfaces to the corresponding cellZone."
- << nl << endl;
+ // 2. Check intersections of all unnamed surfaces
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- findCellZoneGeometric
- (
- neiCc,
- closedNamedSurfaces, // indices of closed surfaces
- ownPatch, // per face patch
- labelList(mesh_.boundaryMesh().size(), 0),// cellZone per patch
+ const labelList testFaces(intersectedFaces());
- cellToZone
- );
- }
+ labelList globalRegion1;
+ labelList globalRegion2;
+ getIntersections
+ (
+ surfaceZonesInfo::getUnnamedSurfaces(surfaces_.surfZones()),
+ neiCc,
+ testFaces,
+ globalRegion1,
+ globalRegion2
+ );
- // Locations in mesh - walk
- findCellZoneInsideWalk
- (
- locationsInMesh,
- zonesInMesh,
- ownPatch, // per face -1 or some index >= 0
- cellToZone
- );
+ labelList neiCellZone;
+ syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
- // Now we will still have some -2 in cellToZone on parts that
- // are unreachable from the locationsInMesh or named surfaces. We only
- // want to keep the interfaces to these cellZones.
+ // 3. Baffle all boundary faces
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Baffle all boundary faces except those on outside of unvisited cells
+ // (cellToZone = -2)
+ // Use patch according to globalRegion1,2
+ // Note: the behaviour is slightly different from version3.0 and earlier
+ // in that it will not create baffles which are outside the meshed
+ // domain. On a medium testcase (motorBike tutorial geometry) this
+ // selects about 20 cells less (out of 120000). These cells are where
+ // there might e.g. be a single cell which is fully unreachable.
+
+ ownPatch.setSize(mesh_.nFaces());
+ ownPatch = -1;
+ neiPatch.setSize(mesh_.nFaces());
+ neiPatch = -1;
+ forAll(testFaces, i)
+ {
+ label faceI = testFaces[i];
+ label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
+ label neiZone = -1;
- labelList neiCellZone;
- syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
+ if (mesh_.isInternalFace(faceI))
+ {
+ neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
+ }
+ else
+ {
+ neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
+ }
- // Only keep baffles where one of the sides has an unset cellToZone (-2)
- for (label faceI=0; faceI < mesh_.nInternalFaces(); faceI++)
+ if (ownZone != -2 || neiZone != -2)
{
- if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
+ if (globalRegion1[faceI] != -1)
{
- label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
- label neiZone = cellToZone[mesh_.faceNeighbour()[faceI]];
-
- if (ownZone != -2 && neiZone != -2)
- {
- ownPatch[faceI] = -1;
- neiPatch[faceI] = -1;
- }
+ ownPatch[faceI] = globalToMasterPatch[globalRegion1[faceI]];
}
- }
- for
- (
- label faceI = mesh_.nInternalFaces();
- faceI < mesh_.nFaces();
- faceI++
- )
- {
- if (ownPatch[faceI] != -1 || neiPatch[faceI] != -1)
+ if (globalRegion2[faceI] != -1)
{
- label ownZone = cellToZone[mesh_.faceOwner()[faceI]];
- label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
-
- if (ownZone != -2 && neiZone != -2)
- {
- ownPatch[faceI] = -1;
- neiPatch[faceI] = -1;
- }
+ neiPatch[faceI] = globalToMasterPatch[globalRegion2[faceI]];
}
}
}
-
// No need to parallel sync since intersection data (surfaceIndex_ etc.)
// already guaranteed to be synced...
// However:
@@ -1417,12 +1417,19 @@ void Foam::meshRefinement::findCellZoneGeometric
forAll(insideSurfaces, cellI)
{
- if (cellToZone[cellI] == -2)
- {
- label surfI = insideSurfaces[cellI];
+ label surfI = insideSurfaces[cellI];
- if (surfI != -1)
+ if (surfI != -1)
+ {
+ if (cellToZone[cellI] == -2)
+ {
+ cellToZone[cellI] = surfaceToCellZone[surfI];
+ }
+ else if (cellToZone[cellI] == -1)
{
+ // ? Allow named surface to override background zone (-1)
+ // This is used in the multiRegionHeater tutorial where the
+ // locationInMesh is inside a named surface.
cellToZone[cellI] = surfaceToCellZone[surfI];
}
}
@@ -1548,7 +1555,9 @@ void Foam::meshRefinement::findCellZoneGeometric
if (namedSurfaceIndex[faceI] == -1 && (ownZone != neiZone))
{
- // Give face the zone of min cell zone
+ // Give face the zone of min cell zone (but only if the
+ // cellZone originated from a closed, named surface)
+
label minZone;
if (ownZone == -1)
{
@@ -1563,7 +1572,13 @@ void Foam::meshRefinement::findCellZoneGeometric
minZone = min(ownZone, neiZone);
}
- namedSurfaceIndex[faceI] = findIndex(surfaceToCellZone, minZone);
+ // Make sure the cellZone originated from a closed surface
+ label geomSurfI = findIndex(surfaceToCellZone, minZone);
+
+ if (geomSurfI != -1)
+ {
+ namedSurfaceIndex[faceI] = geomSurfI;
+ }
}
}
@@ -1601,11 +1616,13 @@ void Foam::meshRefinement::findCellZoneGeometric
minZone = min(ownZone, neiZone);
}
- namedSurfaceIndex[faceI] = findIndex
- (
- surfaceToCellZone,
- minZone
- );
+ // Make sure the cellZone originated from a closed surface
+ label geomSurfI = findIndex(surfaceToCellZone, minZone);
+
+ if (geomSurfI != -1)
+ {
+ namedSurfaceIndex[faceI] = geomSurfI;
+ }
}
}
}
@@ -1616,103 +1633,10 @@ void Foam::meshRefinement::findCellZoneGeometric
}
-void Foam::meshRefinement::findCellZoneInsideWalk
-(
- const labelList& locationSurfaces, // indices of surfaces with inside point
- const labelList& namedSurfaceIndex, // per face index of named surface
- const labelList& surfaceToCellZone, // cell zone index per surface
-
- labelList& cellToZone
-) const
-{
- // Analyse regions. Reuse regionsplit
- boolList blockedFace(mesh_.nFaces());
- //selectSeparatedCoupledFaces(blockedFace);
-
- forAll(namedSurfaceIndex, faceI)
- {
- if (namedSurfaceIndex[faceI] == -1)
- {
- blockedFace[faceI] = false;
- }
- else
- {
- blockedFace[faceI] = true;
- }
- }
- // No need to sync since namedSurfaceIndex already is synced
-
- // Set region per cell based on walking
- regionSplit cellRegion(mesh_, blockedFace);
- blockedFace.clear();
-
-
- // Force calculation of face decomposition (used in findCell)
- (void)mesh_.tetBasePtIs();
-
- const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
-
- // For all locationSurface find the cell
- forAll(locationSurfaces, i)
- {
- label surfI = locationSurfaces[i];
-
- const point& insidePoint = surfZones[surfI].zoneInsidePoint();
-
- // Find the region containing the insidePoint
- label keepRegionI = findRegion
- (
- mesh_,
- cellRegion,
- mergeDistance_*vector(1,1,1),
- insidePoint
- );
-
- Info<< "For surface " << surfaces_.names()[surfI]
- << " found point " << insidePoint
- << " in global region " << keepRegionI
- << " out of " << cellRegion.nRegions() << " regions." << endl;
-
- if (keepRegionI == -1)
- {
- FatalErrorInFunction
- << "Point " << insidePoint
- << " is not inside the mesh." << nl
- << "Bounding box of the mesh:" << mesh_.bounds()
- << exit(FatalError);
- }
-
- // Set all cells with this region
- forAll(cellRegion, cellI)
- {
- if (cellRegion[cellI] == keepRegionI)
- {
- if (cellToZone[cellI] == -2)
- {
- cellToZone[cellI] = surfaceToCellZone[surfI];
- }
- else if (cellToZone[cellI] != surfaceToCellZone[surfI])
- {
- WarningInFunction
- << "Cell " << cellI
- << " at " << mesh_.cellCentres()[cellI]
- << " is inside surface " << surfaces_.names()[surfI]
- << " but already marked as being in zone "
- << cellToZone[cellI] << endl
- << "This can happen if your surfaces are not"
- << " (sufficiently) closed."
- << endl;
- }
- }
- }
- }
-}
-
-
void Foam::meshRefinement::findCellZoneInsideWalk
(
const pointField& locationsInMesh,
- const wordList& zonesInMesh,
+ const labelList& zonesInMesh,
const labelList& faceToZone, // per face -1 or some index >= 0
labelList& cellToZone
@@ -1733,7 +1657,7 @@ void Foam::meshRefinement::findCellZoneInsideWalk
blockedFace[faceI] = true;
}
}
- // No need to sync since blockedFace already is synced
+ // No need to sync since faceToZone already is synced
// Set region per cell based on walking
regionSplit cellRegion(mesh_, blockedFace);
@@ -1752,7 +1676,7 @@ void Foam::meshRefinement::findCellZoneInsideWalk
{
// Get location and index of zone ("none" for cellZone -1)
const point& insidePoint = locationsInMesh[i];
- label zoneID = mesh_.cellZones().findZoneID(zonesInMesh[i]);
+ label zoneID = zonesInMesh[i];
// Find the region containing the insidePoint
label keepRegionI = findRegion
@@ -1763,7 +1687,12 @@ void Foam::meshRefinement::findCellZoneInsideWalk
insidePoint
);
- Info<< "For cellZone " << zonesInMesh[i]
+ Info<< "For cellZone "
+ << (
+ zoneID == -1
+ ? "none"
+ : mesh_.cellZones()[zoneID].name()
+ )
<< " found point " << insidePoint
<< " in global region " << keepRegionI
<< " out of " << cellRegion.nRegions() << " regions." << endl;
@@ -1803,7 +1732,12 @@ void Foam::meshRefinement::findCellZoneInsideWalk
WarningInFunction
<< "Cell " << cellI
<< " at " << mesh_.cellCentres()[cellI]
- << " is inside cellZone " << zonesInMesh[i]
+ << " is inside cellZone "
+ << (
+ zoneID == -1
+ ? "none"
+ : mesh_.cellZones()[zoneID].name()
+ )
<< " from locationInMesh " << insidePoint
<< " but already marked as being in zone "
<< mesh_.cellZones()[cellToZone[cellI]].name()
@@ -1817,54 +1751,32 @@ void Foam::meshRefinement::findCellZoneInsideWalk
}
}
}
+}
- // Check if any unassigned regions
- label nUnzoned = 0;
- forAll(regionToZone, regionI)
- {
- if (regionToZone[regionI] == -2)
- {
- // region that has not been assigned a cellZone
- nUnzoned++;
- }
- }
-
- if (nUnzoned > 0)
- {
- Info<< "Detected " << nUnzoned << " regions in the mesh that do"
- << " not have a locationInMesh." << nl
- << "Per unzoned region displaying a single cell centre:" << endl;
-
- // Determine single location per unzoned region
- pointField regionToLocation
- (
- regionToZone.size(),
- point(GREAT, GREAT, GREAT)
- );
- forAll(cellRegion, cellI)
- {
- label regionI = cellRegion[cellI];
- if (regionToZone[regionI] == -2)
- {
- const point& cc = mesh_.cellCentres()[cellI];
- minMagSqrEqOp<point>()(regionToLocation[regionI], cc);
- }
- }
+void Foam::meshRefinement::findCellZoneInsideWalk
+(
+ const pointField& locationsInMesh,
+ const wordList& zoneNamesInMesh,
+ const labelList& faceToZone, // per face -1 or some index >= 0
- Pstream::listCombineGather(regionToLocation, minMagSqrEqOp<point>());
- Pstream::listCombineScatter(regionToLocation);
+ labelList& cellToZone
+) const
+{
+ const cellZoneMesh& czs = mesh_.cellZones();
- forAll(regionToZone, regionI)
- {
- if (regionToZone[regionI] == -2)
- {
- Info<< '\t' << regionI
- << '\t' << regionToLocation[regionI] << endl;
- }
- }
- Info<< endl;
+ labelList zoneIDs(zoneNamesInMesh.size());
+ forAll(zoneNamesInMesh, i)
+ {
+ zoneIDs[i] = czs.findZoneID(zoneNamesInMesh[i]);
}
+ findCellZoneInsideWalk
+ (
+ locationsInMesh,
+ zoneIDs,
+ faceToZone,
+ cellToZone
+ );
}
@@ -1929,6 +1841,7 @@ bool Foam::meshRefinement::calcRegionToZone
void Foam::meshRefinement::findCellZoneTopo
(
const pointField& locationsInMesh,
+ const labelList& allSurfaceIndex,
const labelList& namedSurfaceIndex,
const labelList& surfaceToCellZone,
labelList& cellToZone
@@ -1942,9 +1855,9 @@ void Foam::meshRefinement::findCellZoneTopo
// Analyse regions. Reuse regionsplit
boolList blockedFace(mesh_.nFaces());
- forAll(namedSurfaceIndex, faceI)
+ forAll(allSurfaceIndex, faceI)
{
- if (namedSurfaceIndex[faceI] == -1)
+ if (allSurfaceIndex[faceI] == -1)
{
blockedFace[faceI] = false;
}
@@ -1970,12 +1883,21 @@ void Foam::meshRefinement::findCellZoneTopo
// Note: could be improved to count number of cells per region.
forAll(cellToZone, cellI)
{
- if (cellToZone[cellI] != -2)
+ label regionI = cellRegion[cellI];
+ if (cellToZone[cellI] != -2 && regionToCellZone[regionI] == -2)
{
- regionToCellZone[cellRegion[cellI]] = cellToZone[cellI];
+ regionToCellZone[regionI] = cellToZone[cellI];
}
}
+ // Synchronise regionToCellZone.
+ // Note:
+ // - region numbers are identical on all processors
+ // - keepRegion is identical ,,
+ // - cellZones are identical ,,
+ Pstream::listCombineGather(regionToCellZone, maxEqOp<label>());
+ Pstream::listCombineScatter(regionToCellZone);
+
// Find the region containing the keepPoint
forAll(locationsInMesh, i)
@@ -2002,7 +1924,10 @@ void Foam::meshRefinement::findCellZoneTopo
<< exit(FatalError);
}
- // Mark default region with zone -1.
+ // Mark default region with zone -1. Note that all regions should
+ // already be matched to a cellZone through the loop over cellToZone.
+ // This is just to mop up any remaining regions. Not sure whether
+ // this is needed?
if (regionToCellZone[keepRegionI] == -2)
{
regionToCellZone[keepRegionI] = -1;
@@ -2096,18 +2021,6 @@ void Foam::meshRefinement::findCellZoneTopo
}
- forAll(regionToCellZone, regionI)
- {
- label zoneI = regionToCellZone[regionI];
-
- if (zoneI == -2)
- {
- FatalErrorInFunction
- << "For region " << regionI << " haven't set cell zone."
- << exit(FatalError);
- }
- }
-
if (debug)
{
forAll(regionToCellZone, regionI)
@@ -2121,17 +2034,30 @@ void Foam::meshRefinement::findCellZoneTopo
// Rework into cellToZone
forAll(cellToZone, cellI)
{
- cellToZone[cellI] = regionToCellZone[cellRegion[cellI]];
+ label regionI = cellRegion[cellI];
+ if (cellToZone[cellI] == -2 && regionToCellZone[regionI] != -2)
+ {
+ cellToZone[cellI] = regionToCellZone[regionI];
+ }
}
}
void Foam::meshRefinement::makeConsistentFaceIndex
(
+ const labelList& surfaceMap,
const labelList& cellToZone,
- labelList& faceToZone
+ labelList& namedSurfaceIndex
) const
{
+ // Make namedSurfaceIndex consistent with cellToZone - clear out any
+ // blocked faces inbetween same cell zone (or background (=-1))
+ // Do not do this for surfaces relating to 'pure' faceZones i.e.
+ // faceZones without a cellZone. Note that we cannot check here
+ // for different cellZones on either side but no namedSurfaceIndex
+ // since cellZones can now originate from locationsInMesh as well
+ // (instead of only through named surfaces)
+
const labelList& faceOwner = mesh_.faceOwner();
const labelList& faceNeighbour = mesh_.faceNeighbour();
@@ -2140,17 +2066,14 @@ void Foam::meshRefinement::makeConsistentFaceIndex
label ownZone = cellToZone[faceOwner[faceI]];
label neiZone = cellToZone[faceNeighbour[faceI]];
- if (ownZone == neiZone && faceToZone[faceI] != -1)
- {
- faceToZone[faceI] = -1;
- }
- else if (ownZone != neiZone && faceToZone[faceI] == -1)
+ if
+ (
+ ownZone == neiZone
+ && namedSurfaceIndex[faceI] != -1
+ && surfaceMap[namedSurfaceIndex[faceI]] == -1
+ )
{
- FatalErrorInFunction
- << "Different cell zones on either side of face " << faceI
- << " at " << mesh_.faceCentres()[faceI]
- << " but face not marked with a surface."
- << abort(FatalError);
+ namedSurfaceIndex[faceI] = -1;
}
}
@@ -2174,17 +2097,14 @@ void Foam::meshRefinement::makeConsistentFaceIndex
label ownZone = cellToZone[faceOwner[faceI]];
label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
- if (ownZone == neiZone && faceToZone[faceI] != -1)
- {
- faceToZone[faceI] = -1;
- }
- else if (ownZone != neiZone && faceToZone[faceI] == -1)
+ if
+ (
+ ownZone == neiZone
+ && namedSurfaceIndex[faceI] != -1
+ && surfaceMap[namedSurfaceIndex[faceI]] == -1
+ )
{
- FatalErrorInFunction
- << "Different cell zones on either side of face "
- << faceI << " at " << mesh_.faceCentres()[faceI]
- << " but face not marked with a surface."
- << abort(FatalError);
+ namedSurfaceIndex[faceI] = -1;
}
}
}
@@ -2194,27 +2114,430 @@ void Foam::meshRefinement::makeConsistentFaceIndex
forAll(pp, i)
{
label faceI = pp.start()+i;
- faceToZone[faceI] = -1;
+
+ if
+ (
+ namedSurfaceIndex[faceI] != -1
+ && surfaceMap[namedSurfaceIndex[faceI]] == -1
+ )
+ {
+ namedSurfaceIndex[faceI] = -1;
+ }
}
}
}
}
-void Foam::meshRefinement::handleSnapProblems
+void Foam::meshRefinement::getIntersections
(
- const snapParameters& snapParams,
- const bool useTopologicalSnapDetection,
- const bool removeEdgeConnectedCells,
- const scalarField& perpendicularAngle,
- const dictionary& motionDict,
- Time& runTime,
- const labelList& globalToMasterPatch,
- const labelList& globalToSlavePatch
-)
-{
- Info<< nl
- << "Introducing baffles to block off problem cells" << nl
+ const labelList& surfacesToTest,
+ const pointField& neiCc,
+ const labelList& testFaces,
+
+ labelList& namedSurfaceIndex,
+ PackedBoolList& posOrientation
+) const
+{
+ const pointField& cellCentres = mesh_.cellCentres();
+ const labelList& faceOwner = mesh_.faceOwner();
+ const labelList& faceNeighbour = mesh_.faceNeighbour();
+
+ namedSurfaceIndex.setSize(mesh_.nFaces());
+ namedSurfaceIndex = -1;
+
+ posOrientation.setSize(mesh_.nFaces());
+ posOrientation = false;
+
+ // Statistics: number of faces per faceZone
+ labelList nSurfFaces(surfaces_.surfZones().size(), 0);
+
+ // Collect segments
+ // ~~~~~~~~~~~~~~~~
+
+ pointField start(testFaces.size());
+ pointField end(testFaces.size());
+
+ forAll(testFaces, i)
+ {
+ label faceI = testFaces[i];
+
+ if (mesh_.isInternalFace(faceI))
+ {
+ start[i] = cellCentres[faceOwner[faceI]];
+ end[i] = cellCentres[faceNeighbour[faceI]];
+ }
+ else
+ {
+ start[i] = cellCentres[faceOwner[faceI]];
+ end[i] = neiCc[faceI-mesh_.nInternalFaces()];
+ }
+ }
+
+ // Extend segments a bit
+ {
+ const vectorField smallVec(ROOTSMALL*(end-start));
+ start -= smallVec;
+ end += smallVec;
+ }
+
+
+ // Do test for intersections
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~
+ // Note that we intersect all intersected faces again. Could reuse
+ // the information already in surfaceIndex_.
+
+ labelList surface1;
+ List<pointIndexHit> hit1;
+ vectorField normal1;
+ labelList surface2;
+ List<pointIndexHit> hit2;
+ vectorField normal2;
+ {
+ labelList region1;
+ labelList region2;
+ surfaces_.findNearestIntersection
+ (
+ surfacesToTest,
+ start,
+ end,
+
+ surface1,
+ hit1,
+ region1,
+ normal1,
+
+ surface2,
+ hit2,
+ region2,
+ normal2
+ );
+ }
+
+ forAll(testFaces, i)
+ {
+ label faceI = testFaces[i];
+ const vector& area = mesh_.faceAreas()[faceI];
+
+ if (surface1[i] != -1)
+ {
+ // If both hit should probably choose 'nearest'
+ if
+ (
+ surface2[i] != -1
+ && (
+ magSqr(hit2[i].hitPoint())
+ < magSqr(hit1[i].hitPoint())
+ )
+ )
+ {
+ namedSurfaceIndex[faceI] = surface2[i];
+ posOrientation[faceI] = ((area&normal2[i]) > 0);
+ nSurfFaces[surface2[i]]++;
+ }
+ else
+ {
+ namedSurfaceIndex[faceI] = surface1[i];
+ posOrientation[faceI] = ((area&normal1[i]) > 0);
+ nSurfFaces[surface1[i]]++;
+ }
+ }
+ else if (surface2[i] != -1)
+ {
+ namedSurfaceIndex[faceI] = surface2[i];
+ posOrientation[faceI] = ((area&normal2[i]) > 0);
+ nSurfFaces[surface2[i]]++;
+ }
+ }
+
+
+ // surfaceIndex might have different surfaces on both sides if
+ // there happen to be a (obviously thin) surface with different
+ // regions between the cell centres. If one is on a named surface
+ // and the other is not this might give problems so sync.
+ syncTools::syncFaceList
+ (
+ mesh_,
+ namedSurfaceIndex,
+ maxEqOp<label>()
+ );
+
+ // Print a bit
+ if (debug)
+ {
+ forAll(nSurfFaces, surfI)
+ {
+ Pout<< "Surface:"
+ << surfaces_.names()[surfI]
+ << " nZoneFaces:" << nSurfFaces[surfI] << nl;
+ }
+ Pout<< endl;
+ }
+}
+
+
+void Foam::meshRefinement::zonify
+(
+ const bool allowFreeStandingZoneFaces,
+ const pointField& locationsInMesh,
+ const wordList& zonesInMesh,
+
+ labelList& cellToZone,
+ labelList& namedSurfaceIndex,
+ PackedBoolList& posOrientation
+) const
+{
+ // Determine zones for cells and faces
+ // cellToZone:
+ // -2 : unset
+ // -1 : not in any zone (zone 'none' or background zone)
+ // >=0 : zoneID
+ // namedSurfaceIndex, posOrientation:
+ // -1 : face not intersected by named surface
+ // >=0 : index of named surface
+ // (and posOrientation: surface normal v.s. face normal)
+
+ const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
+
+ // Swap neighbouring cell centres and cell level
+ labelList neiLevel(mesh_.nFaces()-mesh_.nInternalFaces());
+ pointField neiCc(mesh_.nFaces()-mesh_.nInternalFaces());
+ calcNeighbourData(neiLevel, neiCc);
+
+
+ labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
+
+ // Get map from surface to cellZone (or -1)
+ labelList surfaceToCellZone;
+ if (namedSurfaces.size())
+ {
+ // Get/add cellZones corresponding to surface names
+ surfaceToCellZone = surfaceZonesInfo::addCellZonesToMesh
+ (
+ surfZones,
+ namedSurfaces,
+ mesh_
+ );
+ }
+
+
+ cellToZone.setSize(mesh_.nCells());
+ cellToZone = -2;
+
+ namedSurfaceIndex.clear();
+ posOrientation.clear();
+
+
+
+ // 1. Test all (unnamed & named) surfaces
+
+ labelList globalRegion1;
+ labelList globalRegion2;
+ getIntersections
+ (
+ identity(surfaces_.surfaces().size()), // surfacesToTest,
+ neiCc,
+ intersectedFaces(), // testFaces
+ globalRegion1,
+ globalRegion2
+ );
+
+ if (namedSurfaces.size())
+ {
+ getIntersections
+ (
+ namedSurfaces,
+ neiCc,
+ intersectedFaces(),
+ namedSurfaceIndex,
+ posOrientation
+ );
+ }
+
+
+ // 2. Walk from locationsInMesh. Hard set cellZones.
+
+ if (locationsInMesh.size())
+ {
+ Info<< "Setting cellZones according to locationsInMesh:" << endl;
+
+ labelList locationsZoneIDs(zonesInMesh.size(), -1);
+ forAll(locationsInMesh, i)
+ {
+ const word& name = zonesInMesh[i];
+
+ Info<< "Location : " << locationsInMesh[i] << nl
+ << " cellZone : " << name << endl;
+
+ if (name != "none")
+ {
+ label zoneID = mesh_.cellZones().findZoneID(name);
+ if (zoneID == -1)
+ {
+ FatalErrorIn("meshRefinement::zonify(..)") << "problem"
+ << abort(FatalError);
+ }
+ locationsZoneIDs[i] = zoneID;
+ }
+ }
+ Info<< endl;
+
+
+ // Assign cellZone according to seed points
+ findCellZoneInsideWalk
+ (
+ locationsInMesh, // locations
+ locationsZoneIDs, // index of cellZone (or -1)
+ globalRegion1, // per face -1 (unblocked) or >= 0 (blocked)
+
+ cellToZone
+ );
+ }
+
+
+ // 3. Mark named-surfaces-with-insidePoint. Hard set cellZones.
+
+ labelList locationSurfaces
+ (
+ surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
+ );
+
+ if (locationSurfaces.size())
+ {
+ Info<< "Found " << locationSurfaces.size()
+ << " named surfaces with a provided inside point."
+ << " Assigning cells inside these surfaces"
+ << " to the corresponding cellZone."
+ << nl << endl;
+
+ // Collect per surface the -insidePoint -cellZone
+ pointField insidePoints(locationSurfaces.size());
+ labelList insidePointCellZoneIDs(locationSurfaces.size(), -1);
+ forAll(locationSurfaces, i)
+ {
+ label surfI = locationSurfaces[i];
+ insidePoints[i] = surfZones[surfI].zoneInsidePoint();
+
+ const word& name = surfZones[surfI].cellZoneName();
+ if (name != "none")
+ {
+ label zoneID = mesh_.cellZones().findZoneID(name);
+ if (zoneID == -1)
+ {
+ FatalErrorIn("meshRefinement::zonify(..)") << "problem"
+ << abort(FatalError);
+ }
+ insidePointCellZoneIDs[i] = zoneID;
+ }
+ }
+
+ findCellZoneInsideWalk
+ (
+ insidePoints, // locations
+ insidePointCellZoneIDs, // index of cellZone
+ globalRegion1, // per face -1 (unblocked) or >= 0 (blocked)
+ cellToZone
+ );
+ }
+
+
+ // 4. Mark named-surfaces-with-geometric faces. Do geometric test. Soft set
+ // cellZones. Correct through making consistent.
+
+ // Closed surfaces with cellZone specified.
+ labelList closedNamedSurfaces
+ (
+ surfaceZonesInfo::getClosedNamedSurfaces
+ (
+ surfZones,
+ surfaces_.geometry(),
+ surfaces_.surfaces()
+ )
+ );
+
+ if (closedNamedSurfaces.size())
+ {
+ Info<< "Found " << closedNamedSurfaces.size()
+ << " closed, named surfaces. Assigning cells in/outside"
+ << " these surfaces to the corresponding cellZone."
+ << nl << endl;
+
+ findCellZoneGeometric
+ (
+ neiCc,
+ closedNamedSurfaces, // indices of closed surfaces
+ namedSurfaceIndex, // per face index of named surface
+ surfaceToCellZone, // cell zone index per surface
+
+ cellToZone
+ );
+ }
+
+
+ // 5. Find any unassigned regions (from regionSplit)
+
+ if (namedSurfaces.size())
+ {
+ Info<< "Walking from known cellZones; crossing a faceZone "
+ << "face changes cellZone" << nl << endl;
+
+ findCellZoneTopo
+ (
+ pointField(0),
+ globalRegion1, // To split up cells
+ namedSurfaceIndex, // Step across named surfaces to propagate
+ surfaceToCellZone,
+ cellToZone
+ );
+
+ // Make sure namedSurfaceIndex is unset inbetween same cell zones.
+ if (!allowFreeStandingZoneFaces)
+ {
+ Info<< "Only keeping zone faces inbetween different cellZones."
+ << nl << endl;
+
+ // Surfaces with faceZone but no cellZone
+ labelList standaloneNamedSurfaces
+ (
+ surfaceZonesInfo::getStandaloneNamedSurfaces
+ (
+ surfZones
+ )
+ );
+
+ // Construct map from surface index to index in
+ // standaloneNamedSurfaces (or -1)
+ labelList surfaceMap(surfZones.size(), -1);
+ forAll(standaloneNamedSurfaces, i)
+ {
+ surfaceMap[standaloneNamedSurfaces[i]] = i;
+ }
+
+ makeConsistentFaceIndex
+ (
+ surfaceMap,
+ cellToZone,
+ namedSurfaceIndex
+ );
+ }
+ }
+}
+
+
+void Foam::meshRefinement::handleSnapProblems
+(
+ const snapParameters& snapParams,
+ const bool useTopologicalSnapDetection,
+ const bool removeEdgeConnectedCells,
+ const scalarField& perpendicularAngle,
+ const dictionary& motionDict,
+ Time& runTime,
+ const labelList& globalToMasterPatch,
+ const labelList& globalToSlavePatch
+)
+{
+ Info<< nl
+ << "Introducing baffles to block off problem cells" << nl
<< "----------------------------------------------" << nl
<< endl;
@@ -3796,10 +4119,13 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
wordPairHashTable& zonesToFaceZone
)
{
- const pointField& cellCentres = mesh_.cellCentres();
- const labelList& faceOwner = mesh_.faceOwner();
- const labelList& faceNeighbour = mesh_.faceNeighbour();
+ if (locationsInMesh.size() != zonesInMesh.size())
+ {
+ FatalErrorIn("zonify(..)") << "problem" << abort(FatalError);
+ }
+
const polyBoundaryMesh& patches = mesh_.boundaryMesh();
+ const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
// Swap neighbouring cell centres and cell level
@@ -3808,77 +4134,12 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
calcNeighbourData(neiLevel, neiCc);
-
- // Put the cells into the correct zone
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Zone per cell:
- // -2 : unset
- // -1 : not in any zone (zone 'none')
- // >=0: zoneID
- labelList cellToZone(mesh_.nCells(), -2);
-
-
- // Add from locationsInMesh
- // ~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Filter out keepPoints
- labelList zonedIndices(refinementParameters::zonedLocations(zonesInMesh));
-
- if (zonedIndices.size())
- {
- // Explicitly provided locations and their cellZone. Determine all
- // blocked faces and do walking. Updates cellToZone.
-
- Info<< "Setting cellZones according to locationsInMesh:" << endl;
- //forAll(locationsInMesh, regionI)
- forAll(zonedIndices, j)
- {
- label regionI = zonedIndices[j];
- Info<< "Location : " << locationsInMesh[regionI] << nl
- << " cellZone : " << zonesInMesh[regionI] << endl;
- }
- Info<< endl;
-
- // Test all (unnamed & named) surfaces
- labelList globalRegion1;
- labelList globalRegion2;
- getIntersections
- (
- identity(surfaces_.surfaces().size()), // surfacesToTest,
- neiLevel,
- neiCc,
- intersectedFaces(), // testFaces
- globalRegion1,
- globalRegion2
- );
-
- // Assign cellZone according to seed points
- findCellZoneInsideWalk
- (
- pointField(locationsInMesh, zonedIndices), // locations
- UIndirectList<word>(zonesInMesh, zonedIndices)(),// name of region
- globalRegion1, // per face -1 (unblocked) or >= 0 (blocked)
- cellToZone
- );
- }
-
-
-
- // Mark faces intersecting zoned surfaces
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- const PtrList<surfaceZonesInfo>& surfZones = surfaces_.surfZones();
+ // Add any faceZones, cellZones originating from surface to the mesh
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+ labelList surfaceToCellZone;
+ labelList surfaceToFaceZone;
labelList namedSurfaces(surfaceZonesInfo::getNamedSurfaces(surfZones));
-
-
- //- Face face orientation of zone (only from those intersecting named
- // surfaces)
- PackedBoolList posOrientation(mesh_.nFaces());
- //- Per face index of faceZone or -1
- labelList faceToZone(mesh_.nFaces(), -1);
-
if (namedSurfaces.size())
{
Info<< "Setting cellZones according to named surfaces:" << endl;
@@ -3892,312 +4153,100 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
}
Info<< endl;
-
// Add zones to mesh
- labelList surfaceToFaceZone =
- surfaceZonesInfo::addFaceZonesToMesh
- (
- surfZones,
- namedSurfaces,
- mesh_
- );
-
- labelList surfaceToCellZone =
- surfaceZonesInfo::addCellZonesToMesh
- (
- surfZones,
- namedSurfaces,
- mesh_
- );
-
-
-
-
- // Like surfaceIndex_ but only for named surfaces.
- labelList namedSurfaceIndex(mesh_.nFaces(), -1);
-
- {
- // Statistics: number of faces per faceZone
- labelList nSurfFaces(surfZones.size(), 0);
-
- // Note: for all internal faces? internal + coupled?
- // Since zonify is run after baffling the surfaceIndex_ on baffles
- // is
- // not synchronised across both baffle faces. Fortunately we don't
- // do zonify baffle faces anyway (they are normal boundary faces).
-
- // Collect candidate faces
- // ~~~~~~~~~~~~~~~~~~~~~~~
-
- labelList testFaces(intersectedFaces());
-
- // Collect segments
- // ~~~~~~~~~~~~~~~~
-
- pointField start(testFaces.size());
- pointField end(testFaces.size());
-
- forAll(testFaces, i)
- {
- label faceI = testFaces[i];
-
- if (mesh_.isInternalFace(faceI))
- {
- start[i] = cellCentres[faceOwner[faceI]];
- end[i] = cellCentres[faceNeighbour[faceI]];
- }
- else
- {
- start[i] = cellCentres[faceOwner[faceI]];
- end[i] = neiCc[faceI-mesh_.nInternalFaces()];
- }
- }
-
- // Extend segments a bit
- {
- const vectorField smallVec(ROOTSMALL*(end-start));
- start -= smallVec;
- end += smallVec;
- }
-
-
- // Do test for intersections
- // ~~~~~~~~~~~~~~~~~~~~~~~~~
- // Note that we intersect all intersected faces again. Could reuse
- // the information already in surfaceIndex_.
-
- labelList surface1;
- List<pointIndexHit> hit1;
- vectorField normal1;
- labelList surface2;
- List<pointIndexHit> hit2;
- vectorField normal2;
- {
- labelList region1;
- labelList region2;
- surfaces_.findNearestIntersection
- (
- namedSurfaces,
- start,
- end,
-
- surface1,
- hit1,
- region1,
- normal1,
-
- surface2,
- hit2,
- region2,
- normal2
- );
- }
-
- forAll(testFaces, i)
- {
- label faceI = testFaces[i];
- const vector& area = mesh_.faceAreas()[faceI];
-
- if (surface1[i] != -1)
- {
- // If both hit should probably choose 'nearest'
- if
- (
- surface2[i] != -1
- && (
- magSqr(hit2[i].hitPoint())
- < magSqr(hit1[i].hitPoint())
- )
- )
- {
- namedSurfaceIndex[faceI] = surface2[i];
- posOrientation[faceI] = ((area&normal2[i]) > 0);
- nSurfFaces[surface2[i]]++;
- }
- else
- {
- namedSurfaceIndex[faceI] = surface1[i];
- posOrientation[faceI] = ((area&normal1[i]) > 0);
- nSurfFaces[surface1[i]]++;
- }
- }
- else if (surface2[i] != -1)
- {
- namedSurfaceIndex[faceI] = surface2[i];
- posOrientation[faceI] = ((area&normal2[i]) > 0);
- nSurfFaces[surface2[i]]++;
- }
- }
-
-
- // surfaceIndex might have different surfaces on both sides if
- // there happen to be a (obviously thin) surface with different
- // regions between the cell centres. If one is on a named surface
- // and the other is not this might give problems so sync.
- syncTools::syncFaceList
- (
- mesh_,
- namedSurfaceIndex,
- maxEqOp<label>()
- );
-
- // Print a bit
- if (debug)
- {
- forAll(nSurfFaces, surfI)
- {
- Pout<< "Surface:"
- << surfaces_.names()[surfI]
- << " nZoneFaces:" << nSurfFaces[surfI] << nl;
- }
- Pout<< endl;
- }
- }
-
-
- // Now we have for all faces the intersection with a named surfaces (or
- // -1). Use this information to set zone of cells
-
- // Set using geometric test
- // ~~~~~~~~~~~~~~~~~~~~~~~~
-
- // Closed surfaces with cellZone specified.
- labelList closedNamedSurfaces
+ surfaceToCellZone = surfaceZonesInfo::addCellZonesToMesh
(
- surfaceZonesInfo::getClosedNamedSurfaces
- (
- surfZones,
- surfaces_.geometry(),
- surfaces_.surfaces()
- )
+ surfZones,
+ namedSurfaces,
+ mesh_
);
-
- if (closedNamedSurfaces.size())
- {
- Info<< "Found " << closedNamedSurfaces.size()
- << " closed, named surfaces. Assigning cells in/outside"
- << " these surfaces to the corresponding cellZone."
- << nl << endl;
-
- findCellZoneGeometric
- (
- neiCc,
- closedNamedSurfaces, // indices of closed surfaces
- namedSurfaceIndex, // per face index of named surface
- surfaceToCellZone, // cell zone index per surface
-
- cellToZone
- );
- }
-
-
- // Set using provided locations
- // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-
- labelList locationSurfaces
+ surfaceToFaceZone = surfaceZonesInfo::addFaceZonesToMesh
(
- surfaceZonesInfo::getInsidePointNamedSurfaces(surfZones)
+ surfZones,
+ namedSurfaces,
+ mesh_
);
+ }
- if (locationSurfaces.size())
- {
- Info<< "Found " << locationSurfaces.size()
- << " named surfaces with a provided inside point."
- << " Assigning cells inside these surfaces"
- << " to the corresponding cellZone."
- << nl << endl;
-
- findCellZoneInsideWalk
- (
- locationSurfaces, // indices of closed surfaces
- namedSurfaceIndex, // per face index of named surface
- surfaceToCellZone, // cell zone index per surface
-
- cellToZone
- );
- }
- // Set using walking
- // ~~~~~~~~~~~~~~~~~
+ // Put the cells into the correct zone
+ // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
- Info<< "Walking from " << locationsInMesh.size()
- << " locations-in-mesh "
- << " to assign cellZones "
- << "- crossing a faceZone face changes cellZone" << nl << endl;
+ // Zone per cell:
+ // -2 : unset
+ // -1 : not in any zone (zone 'none')
+ // >=0: zoneID
+ // namedSurfaceIndex:
+ // -1 : face not intersecting a named surface
+ // >=0 : index of named surface
+ labelList cellToZone;
+ labelList namedSurfaceIndex;
+ PackedBoolList posOrientation;
+ zonify
+ (
+ allowFreeStandingZoneFaces,
+ locationsInMesh,
+ zonesInMesh,
- labelList zoneIDs(refinementParameters::unzonedLocations(zonesInMesh));
+ cellToZone,
+ namedSurfaceIndex,
+ posOrientation
+ );
- findCellZoneTopo
- (
- pointField(locationsInMesh, zoneIDs),
- namedSurfaceIndex,
- surfaceToCellZone,
- cellToZone
- );
+ // Convert namedSurfaceIndex (index of named surfaces) to
+ // actual faceZone index
+ //- Per face index of faceZone or -1
+ labelList faceToZone(mesh_.nFaces(), -1);
- // Make sure namedSurfaceIndex is unset inbetween same cell zones.
- if (!allowFreeStandingZoneFaces)
+ forAll(namedSurfaceIndex, faceI)
+ {
+ label surfI = namedSurfaceIndex[faceI];
+ if (surfI != -1)
{
- Info<< "Only keeping zone faces inbetween different cellZones."
- << nl << endl;
-
- makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
+ faceToZone[faceI] = surfaceToFaceZone[surfI];
}
+ }
- // Convert namedSurfaceIndex (index of named surfaces) to
- // actual faceZone index
+ if (debug&MESH)
+ {
+ const_cast<Time&>(mesh_.time())++;
+ Pout<< "Writing cell zone allocation on mesh to time "
+ << timeName() << endl;
+ write
+ (
+ debugType(debug),
+ writeType(writeLevel() | WRITEMESH),
+ mesh_.time().path()/"cellToZone"
+ );
+ volScalarField volCellToZone
+ (
+ IOobject
+ (
+ "cellToZone",
+ mesh_.time().timeName(),
+ mesh_,
+ IOobject::NO_READ,
+ IOobject::AUTO_WRITE,
+ false
+ ),
+ mesh_,
+ dimensionedScalar("zero", dimless, 0),
+ zeroGradientFvPatchScalarField::typeName
+ );
- forAll(namedSurfaceIndex, faceI)
+ forAll(cellToZone, cellI)
{
- label surfI = namedSurfaceIndex[faceI];
- if (surfI != -1)
- {
- faceToZone[faceI] = surfaceToFaceZone[surfI];
- }
+ volCellToZone[cellI] = cellToZone[cellI];
}
+ volCellToZone.write();
}
- //if (debug&MESH)
- //{
- // const_cast<Time&>(mesh_.time())++;
- // Pout<< "Writing mesh to time " << timeName() << endl;
- // write
- // (
- // debugType(debug),
- // writeType(writeLevel() | WRITEMESH),
- // mesh_.time().path()/"cellToZone"
- // );
- // volScalarField volCellToZone
- // (
- // IOobject
- // (
- // "cellToZone",
- // mesh_.time().timeName(),
- // mesh_,
- // IOobject::NO_READ,
- // IOobject::AUTO_WRITE,
- // false
- // ),
- // mesh_,
- // dimensionedScalar("zero", dimless, 0),
- // zeroGradientFvPatchScalarField::typeName
- // );
- //
- // forAll(cellToZone, cellI)
- // {
- // volCellToZone[cellI] = cellToZone[cellI];
- // }
- // volCellToZone.write();
- //}
-
-
-
// Allocate and assign faceZones from cellZones
// ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -4208,7 +4257,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
- if (faceToZone[faceI] == -1)
+ if (faceToZone[faceI] == -1) // Not named surface
{
// Face not yet in a faceZone. (it might already have been
// done so by a 'named' surface). Check if inbetween different
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
index 1962ced8b3ba68157aa126800b4f7b1608b35dad..7021d0cf4522ab1871260b144bb3c9c2dbf9da36 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.C
@@ -226,6 +226,32 @@ Foam::labelList Foam::surfaceZonesInfo::getNamedSurfaces
}
+Foam::labelList Foam::surfaceZonesInfo::getStandaloneNamedSurfaces
+(
+ const PtrList<surfaceZonesInfo>& surfList
+)
+{
+ labelList namedSurfaces(surfList.size());
+
+ label namedI = 0;
+ forAll(surfList, surfI)
+ {
+ if
+ (
+ surfList.set(surfI)
+ && surfList[surfI].faceZoneName().size()
+ && !surfList[surfI].cellZoneName().size()
+ )
+ {
+ namedSurfaces[namedI++] = surfI;
+ }
+ }
+ namedSurfaces.setSize(namedI);
+
+ return namedSurfaces;
+}
+
+
Foam::labelList Foam::surfaceZonesInfo::getClosedNamedSurfaces
(
const PtrList<surfaceZonesInfo>& surfList,
diff --git a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H
index 3b5a73c3f5c3c5778b2fa17565eaeca5f34f2444..a6af0f24113f74154cbf7b0dd73f9833a77c9e61 100644
--- a/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H
+++ b/src/mesh/autoMesh/autoHexMesh/refinementSurfaces/surfaceZonesInfo.H
@@ -191,6 +191,12 @@ public:
const PtrList<surfaceZonesInfo>& surfList
);
+ //- Get indices of named surfaces without a cellZone
+ static labelList getStandaloneNamedSurfaces
+ (
+ const PtrList<surfaceZonesInfo>& surfList
+ );
+
//- Get indices of surfaces with a cellZone that are closed and
// have 'inside' or 'outside' selection.
static labelList getClosedNamedSurfaces