Commit c7642d0c authored by Henry Weller's avatar Henry Weller
Browse files

meshRefinementBaffles: Correct faceZone orientation

Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1479
Patches provided by Bruno Santos based on the work of Mattijs Janssens
parent 62945ab1
......@@ -369,7 +369,7 @@ private:
const labelList& globalToSlavePatch
) const;
//- Determine patches for baffles
//- Determine patches for baffles on all intersected unnamed faces
void getBafflePatches
(
const labelList& globalToMasterPatch,
......@@ -514,6 +514,8 @@ private:
labelList& cellToZone
) const;
//- Make namedSurfaceIndex consistent with cellToZone
// - clear out any blocked faces inbetween same cell zone.
void makeConsistentFaceIndex
(
const labelList& cellToZone,
......@@ -536,11 +538,11 @@ private:
// Some patch utilities
//- Get all faces in namedSurfaceIndex that have no cellZone on
//- Get all faces in faceToZone that have no cellZone on
// either side.
labelList freeStandingBaffleFaces
(
const labelList& namedSurfaceIndex,
const labelList& faceToZone,
const labelList& cellToZone,
const labelList& neiCellZone
) const;
......
......@@ -51,9 +51,6 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// Repatches external face or creates baffle for internal face
// with user specified patches (might be different for both sides).
// Returns label of added face.
Foam::label Foam::meshRefinement::createBaffle
(
const label faceI,
......@@ -128,52 +125,11 @@ Foam::label Foam::meshRefinement::createBaffle
}
//// Check if we are a boundary face and normal of surface does
//// not align with test vector. In this case there'd probably be
//// a freestanding 'baffle' so we might as well not create it.
//// Note that since it is not a proper baffle we cannot detect it
//// afterwards so this code cannot be merged with the
//// filterDuplicateFaces code.
//bool Foam::meshRefinement::validBaffleTopology
//(
// const label faceI,
// const vector& n1,
// const vector& n2,
// const vector& testDir
//) const
//{
//
// label patchI = mesh_.boundaryMesh().whichPatch(faceI);
// if (patchI == -1 || mesh_.boundaryMesh()[patchI].coupled())
// {
// return true;
// }
// else if (mag(n1&n2) > cos(degToRad(30)))
// {
// // Both normals aligned. Check that test vector perpendicularish to
// // surface normal
// scalar magTestDir = mag(testDir);
// if (magTestDir > VSMALL)
// {
// if (mag(n1&(testDir/magTestDir)) < cos(degToRad(45)))
// {
// //Pout<< "** disabling baffling face "
// // << mesh_.faceCentres()[faceI] << endl;
// return false;
// }
// }
// }
// return true;
//}
// Determine patches for baffles on all intersected unnamed faces
void Foam::meshRefinement::getBafflePatches
(
const labelList& globalToMasterPatch,
const labelList& neiLevel,
const pointField& neiCc,
labelList& ownPatch,
labelList& neiPatch
) const
......@@ -386,7 +342,6 @@ Foam::Map<Foam::labelPair> Foam::meshRefinement::getZoneBafflePatches
}
// Create baffle for every face where ownPatch != -1
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles
(
const labelList& ownPatch,
......@@ -659,17 +614,16 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles
}
// Extract those baffles (duplicate) faces that are on the edge of a baffle
// region. These are candidates for merging.
// Done by counting the number of baffles faces per mesh edge. If edge
// has 2 boundary faces and both are baffle faces it is the edge of a baffle
// region.
Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
(
const List<labelPair>& couples,
const scalar planarAngle
) const
{
// Done by counting the number of baffles faces per mesh edge. If edge
// has 2 boundary faces and both are baffle faces it is the edge of a baffle
// region.
// All duplicate faces on edge of the patch are to be merged.
// So we count for all edges of duplicate faces how many duplicate
// faces use them.
......@@ -918,7 +872,6 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles
}
// Merge baffles. Gets pairs of faces.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles
(
const List<labelPair>& couples
......@@ -1420,10 +1373,6 @@ bool Foam::meshRefinement::calcRegionToZone
}
// Finds region per cell. Assumes:
// - region containing keepPoint does not go into a cellZone
// - all other regions can be found by crossing faces marked in
// namedSurfaceIndex.
void Foam::meshRefinement::findCellZoneTopo
(
const point& keepPoint,
......@@ -1432,6 +1381,11 @@ void Foam::meshRefinement::findCellZoneTopo
labelList& cellToZone
) const
{
// Assumes:
// - region containing keepPoint does not go into a cellZone
// - all other regions can be found by crossing faces marked in
// namedSurfaceIndex.
// Analyse regions. Reuse regionsplit
boolList blockedFace(mesh_.nFaces());
......@@ -1629,8 +1583,6 @@ void Foam::meshRefinement::findCellZoneTopo
}
// Make namedSurfaceIndex consistent with cellToZone - clear out any blocked
// faces inbetween same cell zone.
void Foam::meshRefinement::makeConsistentFaceIndex
(
const labelList& cellToZone,
......@@ -1810,7 +1762,7 @@ void Foam::meshRefinement::handleSnapProblems
Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
(
const labelList& namedSurfaceIndex,
const labelList& faceToZone,
const labelList& cellToZone,
const labelList& neiCellZone
) const
......@@ -1820,17 +1772,30 @@ Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
const labelList& faceNeighbour = mesh_.faceNeighbour();
DynamicList<label> faceLabels(mesh_.nFaces()/20);
// We want to pick up the faces to orient. These faces come in
// two variants:
// - faces originating from stand-alone faceZones
// (these will most likely have no cellZone on either side so
// ownZone and neiZone both -1)
// - sticky-up faces originating from a 'bulge' in a outside of
// a cellZone. These will have the same cellZone on either side.
// How to orient these is not really clearly defined so do them
// same as stand-alone faceZone faces for now. (Normally these will
// already have been removed by the 'allowFreeStandingZoneFaces=false'
// default setting)
// Note that argument neiCellZone will have -1 on uncoupled boundaries.
DynamicList<label> faceLabels(mesh_.nFaces()/100);
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
if (faceToZone[faceI] != -1)
{
// Free standing baffle?
label ownZone = cellToZone[faceOwner[faceI]];
label neiZone = cellToZone[faceNeighbour[faceI]];
if (max(ownZone, neiZone) == -1)
if (ownZone == neiZone)
{
faceLabels.append(faceI);
}
......@@ -1843,13 +1808,12 @@ Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces
forAll(pp, i)
{
label faceI = pp.start()+i;
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
if (faceToZone[faceI] != -1)
{
// Free standing baffle?
label ownZone = cellToZone[faceOwner[faceI]];
label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
if (max(ownZone, neiZone) == -1)
if (ownZone == neiZone)
{
faceLabels.append(faceI);
}
......@@ -2303,7 +2267,6 @@ void Foam::meshRefinement::consistentOrientation
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// Split off unreachable areas of mesh.
void Foam::meshRefinement::baffleAndSplitMesh
(
const bool doHandleSnapProblems,
......@@ -2494,7 +2457,6 @@ void Foam::meshRefinement::baffleAndSplitMesh
}
// Split off (with optional buffer layers) unreachable areas of mesh.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
(
const label nBufferLayers,
......@@ -2778,8 +2740,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh
}
// Find boundary points that connect to more than one cell region and
// split them.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints
(
const localPointRegion& regionSide
......@@ -2833,8 +2793,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints
}
// Find boundary points that connect to more than one cell region and
// split them.
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
{
// Analyse which points need to be duplicated
......@@ -2844,7 +2802,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints()
}
// Zoning
Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
(
const point& keepPoint,
......@@ -3140,28 +3097,46 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
makeConsistentFaceIndex(cellToZone, namedSurfaceIndex);
}
//- Per face index of faceZone or -1
labelList faceToZone(mesh_.nFaces(), -1);
// Convert namedSurfaceIndex (index of named surfaces) to
// actual faceZone index
forAll(namedSurfaceIndex, faceI)
{
label surfI = namedSurfaceIndex[faceI];
if (surfI != -1)
{
faceToZone[faceI] = surfaceToFaceZone[surfI];
}
}
// Topochange container
polyTopoChange meshMod(mesh_);
// Get coupled neighbour cellZone. Set to -1 on non-coupled patches.
labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1);
labelList neiCellZone;
syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone);
forAll(patches, patchI)
{
const polyPatch& pp = patches[patchI];
if (pp.coupled())
if (!pp.coupled())
{
label bFaceI = pp.start()-mesh_.nInternalFaces();
forAll(pp, i)
{
label faceI = pp.start()+i;
neiCellZone[faceI-mesh_.nInternalFaces()] =
cellToZone[mesh_.faceOwner()[faceI]];
neiCellZone[bFaceI++] = -1;
}
}
}
syncTools::swapBoundaryFaceList(mesh_, neiCellZone);
// Get per face whether is it master (of a coupled set of faces)
const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_));
......@@ -3188,7 +3163,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
mesh_.faces(),
freeStandingBaffleFaces
(
namedSurfaceIndex,
faceToZone,
cellToZone,
neiCellZone
)
......@@ -3202,18 +3177,25 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
Info<< "Detected " << nFreeStanding << " free-standing zone faces"
<< endl;
if (debug)
{
OBJstream str(mesh_.time().path()/"freeStanding.obj");
str.write(patch.localFaces(), patch.localPoints(), false);
}
// Detect non-manifold edges
labelList nMasterFacesPerEdge;
calcPatchNumMasterFaces(isMasterFace, patch, nMasterFacesPerEdge);
// Mark zones. Even a single original surface might create multiple
// disconnected/non-manifold-connected zones
labelList faceToZone;
labelList faceToConnectedZone;
const label nZones = markPatchZones
(
patch,
nMasterFacesPerEdge,
faceToZone
faceToConnectedZone
);
Map<label> nPosOrientation(2*nZones);
......@@ -3230,7 +3212,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
isMasterFace,
patch,
nMasterFacesPerEdge,
faceToZone,
faceToConnectedZone,
nPosOrientation,
meshFlipMap
......@@ -3254,7 +3236,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
n = -1;
}
nPosOrientation.find(faceToZone[faceI])() += n;
nPosOrientation.find(faceToConnectedZone[faceI])() += n;
}
}
Pstream::mapCombineGather(nPosOrientation, plusEqOp<label>());
......@@ -3281,7 +3263,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
isMasterFace,
patch,
nMasterFacesPerEdge,
faceToZone,
faceToConnectedZone,
nPosOrientation,
meshFlipMap
......@@ -3295,30 +3277,31 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++)
{
label surfI = namedSurfaceIndex[faceI];
label faceZoneI = faceToZone[faceI];
if (surfI != -1)
if (faceZoneI != -1)
{
// Orient face zone to have slave cells in max cell zone.
// Note: logic to use flipMap should be consistent with logic
// to pick up the freeStandingBaffleFaces!
label ownZone = cellToZone[faceOwner[faceI]];
label neiZone = cellToZone[faceNeighbour[faceI]];
bool flip;
label maxZone = max(ownZone, neiZone);
if (maxZone == -1)
if (ownZone == neiZone)
{
// free-standing face. Use geometrically derived orientation
flip = meshFlipMap[faceI];
}
else if (ownZone == maxZone)
{
flip = false;
}
else
{
flip = true;
flip =
(
ownZone == -1
|| (neiZone != -1 && ownZone > neiZone)
);
}
meshMod.setAction
......@@ -3332,7 +3315,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
false, // face flip
-1, // patch for face
false, // remove from zone
surfaceToFaceZone[surfI], // zone for face
faceZoneI, // zone for face
flip // face flip in zone
)
);
......@@ -3349,35 +3332,27 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
forAll(pp, i)
{
label surfI = namedSurfaceIndex[faceI];
label faceZoneI = faceToZone[faceI];
if (surfI != -1)
if (faceZoneI != -1)
{
label ownZone = cellToZone[faceOwner[faceI]];
label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()];
bool flip;
label maxZone = max(ownZone, neiZone);
if (maxZone == -1)
if (ownZone == neiZone)
{
// free-standing face. Use geometrically derived orientation
flip = meshFlipMap[faceI];
}
else if (ownZone == neiZone)
{
// Free-standing zone face or coupled boundary. Keep master
// face unflipped.
flip = !isMasterFace[faceI];
}
else if (neiZone == maxZone)
{
flip = true;
}
else
{
flip = false;
flip =
(
ownZone == -1
|| (neiZone != -1 && ownZone > neiZone)
);
}
meshMod.setAction
......@@ -3391,7 +3366,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify
false, // face flip
patchI, // patch for face
false, // remove from zone
surfaceToFaceZone[surfI], // zone for face
faceZoneI, // zone for face
flip // face flip in zone
)
);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment