diff --git a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C index 5cdab910ac74b0cc8a6bc406c25d81c3a37dcca9..d9ec65959ca7c4bebfa9fda4430b9dba3b39275a 100644 --- a/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C +++ b/src/mesh/autoMesh/autoHexMesh/meshRefinement/meshRefinementBaffles.C @@ -51,6 +51,9 @@ 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,6 +131,46 @@ 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, @@ -342,11 +385,11 @@ Foam::Map<Foam::labelPair> Foam::meshRefinement::getZoneBafflePatches } } } - return bafflePatch; } +// Create baffle for every face where ownPatch != -1 Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles ( const labelList& ownPatch, @@ -466,7 +509,6 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles } } } - // Pick up neighbour side of baffle (added faces) forAll(faceMap, faceI) { @@ -493,6 +535,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createBaffles void Foam::meshRefinement::checkZoneFaces() const { const faceZoneMesh& fZones = mesh_.faceZones(); + const polyBoundaryMesh& pbm = mesh_.boundaryMesh(); forAll(pbm, patchI) @@ -621,26 +664,27 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::createZoneBaffles Info<< "Created " << nZoneFaces << " baffles in = " << mesh_.time().cpuTimeIncrement() << " s\n" << nl << endl; } - return map; } +// 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. labelList nBafflesPerEdge(mesh_.nEdges(), 0); + // This algorithm is quite tricky. We don't want to use edgeFaces and // also want it to run in parallel so it is now an algorithm over // all (boundary) faces instead. @@ -654,8 +698,10 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles // So now any edge that is used by baffle faces only will have the // value 2*1000000+2*1. + const label baffleValue = 1000000; + // Count number of boundary faces per edge // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -683,9 +729,11 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles } } + DynamicList<label> fe0; DynamicList<label> fe1; + // Count number of duplicate boundary faces per edge // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -879,6 +927,7 @@ Foam::List<Foam::labelPair> Foam::meshRefinement::freeStandingBaffles } +// Merge baffles. Gets pairs of faces. Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles ( const List<labelPair>& couples @@ -1010,6 +1059,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::mergeBaffles } +// Finds region per cell for cells inside closed named surfaces void Foam::meshRefinement::findCellZoneGeometric ( const pointField& neiCc, @@ -1229,6 +1279,7 @@ 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 { @@ -1247,7 +1298,6 @@ void Foam::meshRefinement::findCellZoneInsideWalk blockedFace[faceI] = true; } } - // No need to sync since namedSurfaceIndex already is synced // Set region per cell based on walking @@ -1383,11 +1433,14 @@ bool Foam::meshRefinement::calcRegionToZone } } } - return changed; } +// 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, @@ -1396,11 +1449,6 @@ 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()); @@ -1604,6 +1652,8 @@ 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, @@ -1793,21 +1843,7 @@ Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces const labelList& faceNeighbour = mesh_.faceNeighbour(); - // 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); + DynamicList<label> faceLabels(mesh_.nFaces()/20); for (label faceI = 0; faceI < mesh_.nInternalFaces(); faceI++) { @@ -1817,7 +1853,7 @@ Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces // Free standing baffle? label ownZone = cellToZone[faceOwner[faceI]]; label neiZone = cellToZone[faceNeighbour[faceI]]; - if (ownZone == neiZone) + if (max(ownZone, neiZone) == -1) { faceLabels.append(faceI); } @@ -1836,14 +1872,13 @@ Foam::labelList Foam::meshRefinement::freeStandingBaffleFaces // Free standing baffle? label ownZone = cellToZone[faceOwner[faceI]]; label neiZone = neiCellZone[faceI-mesh_.nInternalFaces()]; - if (ownZone == neiZone) + if (max(ownZone, neiZone) == -1) { faceLabels.append(faceI); } } } } - return faceLabels.shrink(); } @@ -2100,6 +2135,7 @@ void Foam::meshRefinement::consistentOrientation } + DynamicList<label> changedEdges; DynamicList<patchFaceOrientation> changedInfo; @@ -2188,6 +2224,7 @@ void Foam::meshRefinement::consistentOrientation } + // Walk PatchEdgeFaceWave < @@ -2289,6 +2326,7 @@ void Foam::meshRefinement::consistentOrientation // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // +// Split off unreachable areas of mesh. void Foam::meshRefinement::baffleAndSplitMesh ( const bool doHandleSnapProblems, @@ -2479,6 +2517,7 @@ void Foam::meshRefinement::baffleAndSplitMesh } +// Split off (with optional buffer layers) unreachable areas of mesh. Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::splitMesh ( const label nBufferLayers, @@ -2765,6 +2804,8 @@ 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 @@ -2818,6 +2859,8 @@ 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 @@ -2827,6 +2870,7 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::dupNonManifoldPoints() } +// Zoning Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify ( const point& keepPoint, @@ -3128,23 +3172,22 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify // Get coupled neighbour cellZone. Set to -1 on non-coupled patches. - labelList neiCellZone; - syncTools::swapBoundaryCellList(mesh_, cellToZone, neiCellZone); + labelList neiCellZone(mesh_.nFaces()-mesh_.nInternalFaces(), -1); forAll(patches, patchI) { const polyPatch& pp = patches[patchI]; - if (!pp.coupled()) + if (pp.coupled()) { - label bFaceI = pp.start()-mesh_.nInternalFaces(); forAll(pp, i) { - neiCellZone[bFaceI++] = -1; + label faceI = pp.start()+i; + neiCellZone[faceI-mesh_.nInternalFaces()] = + cellToZone[mesh_.faceOwner()[faceI]]; } } } - - + syncTools::swapBoundaryFaceList(mesh_, neiCellZone); // Get per face whether is it master (of a coupled set of faces) const PackedBoolList isMasterFace(syncTools::getMasterFaces(mesh_)); @@ -3185,13 +3228,6 @@ 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); @@ -3290,19 +3326,22 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify if (surfI != -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; - if (ownZone == neiZone) + label maxZone = max(ownZone, neiZone); + + if (maxZone == -1) { // free-standing face. Use geometrically derived orientation flip = meshFlipMap[faceI]; } + else if (ownZone == maxZone) + { + flip = false; + } else { flip = true; @@ -3345,18 +3384,26 @@ Foam::autoPtr<Foam::mapPolyMesh> Foam::meshRefinement::zonify bool flip; - if (ownZone == neiZone) + label maxZone = max(ownZone, neiZone); + + if (maxZone == -1) { // 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 = - ( - ownZone == -1 - || (neiZone != -1 && ownZone > neiZone) - ); + flip = false; } meshMod.setAction