diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files index e428fc63dc..c383e52fda 100644 --- a/src/meshTools/Make/files +++ b/src/meshTools/Make/files @@ -211,6 +211,7 @@ $(faceZoneSources)/setsToFaceZone/setsToFaceZone.C $(faceZoneSources)/setToFaceZone/setToFaceZone.C $(faceZoneSources)/setAndNormalToFaceZone/setAndNormalToFaceZone.C $(faceZoneSources)/searchableSurfaceToFaceZone/searchableSurfaceToFaceZone.C +$(faceZoneSources)/planeToFaceZone/planeToFaceZone.C cellZoneSources = sets/cellZoneSources $(cellZoneSources)/setToCellZone/setToCellZone.C diff --git a/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.C b/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.C new file mode 100644 index 0000000000..2fbd45576d --- /dev/null +++ b/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.C @@ -0,0 +1,413 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2020 OpenFOAM Foundation + Copyright (C) 2020 OpenCFD Ltd. +------------------------------------------------------------------------------- +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 3 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, see . + +\*---------------------------------------------------------------------------*/ + +#include "planeToFaceZone.H" +#include "polyMesh.H" +#include "faceZoneSet.H" +#include "uindirectPrimitivePatch.H" +#include "PatchTools.H" +#include "syncTools.H" +#include "addToRunTimeSelectionTable.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(planeToFaceZone, 0); + addToRunTimeSelectionTable(topoSetSource, planeToFaceZone, word); + addToRunTimeSelectionTable(topoSetSource, planeToFaceZone, istream); +} + + +Foam::topoSetSource::addToUsageTable Foam::planeToFaceZone::usage_ +( + planeToFaceZone::typeName, + "\n Usage: planeToFaceZone (px py pz) (nx ny nz) include\n\n" + " Select faces for which the adjacent cell centres lie on opposite " + " of a plane\n\n" +); + + +const Foam::Enum +< + Foam::planeToFaceZone::faceZoneAction +> +Foam::planeToFaceZone::faceZoneActionNames_ +({ + { faceZoneAction::ALL, "all" }, + { faceZoneAction::CLOSEST, "closest" }, +}); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::planeToFaceZone::combine(faceZoneSet& fzSet, const bool add) const +{ + // Mark all cells with centres above the plane + boolList cellIsAbovePlane(mesh_.nCells()); + forAll(mesh_.cells(), celli) + { + cellIsAbovePlane[celli] = + ((mesh_.cellCentres()[celli] - point_) & normal_) > 0; + } + + // Mark all faces that sit between cells above and below the plane + boolList faceIsOnPlane(mesh_.nFaces()); + forAll(mesh_.faceNeighbour(), facei) + { + faceIsOnPlane[facei] = + cellIsAbovePlane[mesh_.faceOwner()[facei]] + != cellIsAbovePlane[mesh_.faceNeighbour()[facei]]; + } + forAll(mesh_.boundaryMesh(), patchi) + { + const polyPatch& patch = mesh_.boundaryMesh()[patchi]; + forAll(patch, patchFacei) + { + const label facei = patch.start() + patchFacei; + faceIsOnPlane[facei] = + patch.coupled() && cellIsAbovePlane[mesh_.faceOwner()[facei]]; + } + } + syncTools::syncFaceList(mesh_, faceIsOnPlane, notEqualOp()); + + // Convert marked faces to a list of indices + labelList newSetFaces(findIndices(faceIsOnPlane, true)); + + // If constructing a single contiguous set, remove all faces except those + // connected to the contiguous region closest to the specified point + if (option_ == faceZoneAction::CLOSEST) + { + // Step 1: Get locally contiguous regions for the new face set and the + // total number of regions across all processors. + labelList newSetFaceRegions(newSetFaces.size(), -1); + label nRegions = -1; + { + // Create a patch of the set faces + const uindirectPrimitivePatch newSetPatch + ( + UIndirectList(mesh_.faces(), newSetFaces), + mesh_.points() + ); + + // Get the region ID-s and store the total number of regions on + // each processor + labelList procNRegions(Pstream::nProcs(), -1); + procNRegions[Pstream::myProcNo()] = + PatchTools::markZones + ( + newSetPatch, + boolList(newSetPatch.nEdges(), false), + newSetFaceRegions + ); + Pstream::gatherList(procNRegions); + Pstream::scatterList(procNRegions); + + // Cumulative sum the number of regions on each processor to get an + // offset which makes the local region ID-s globally unique + labelList procRegionOffset(Pstream::nProcs(), 0); + for (label proci = 1; proci < Pstream::nProcs(); ++proci) + { + procRegionOffset[proci] += + procRegionOffset[proci - 1] + + procNRegions[proci - 1]; + } + + // Apply the offset + forAll(newSetFaces, newSetFacei) + { + newSetFaceRegions[newSetFacei] += + procRegionOffset[Pstream::myProcNo()]; + } + + // Store the total number of regions across all processors + nRegions = procRegionOffset.last() + procNRegions.last(); + } + + // Step 2: Create a region map which combines regions which are + // connected across coupled interfaces + labelList regionMap(identity(nRegions)); + { + // Put region labels on connected boundary edges and synchronise to + // create a list of all regions connected to a given edge + labelListList meshEdgeRegions(mesh_.nEdges(), labelList()); + forAll(newSetFaces, newSetFacei) + { + const label facei = newSetFaces[newSetFacei]; + const label regioni = newSetFaceRegions[newSetFacei]; + forAll(mesh_.faceEdges()[facei], faceEdgei) + { + const label edgei = mesh_.faceEdges()[facei][faceEdgei]; + meshEdgeRegions[edgei] = labelList(1, regioni); + } + } + syncTools::syncEdgeList + ( + mesh_, + meshEdgeRegions, + globalMeshData::ListPlusEqOp(), + labelList() + ); + + // Combine edge regions to create a list of what regions a given + // region is connected to + List regionRegions(nRegions); + forAll(newSetFaces, newSetFacei) + { + const label facei = newSetFaces[newSetFacei]; + const label regioni = newSetFaceRegions[newSetFacei]; + forAll(mesh_.faceEdges()[facei], faceEdgei) + { + const label edgei = mesh_.faceEdges()[facei][faceEdgei]; + forAll(meshEdgeRegions[edgei], edgeRegioni) + { + if (meshEdgeRegions[edgei][edgeRegioni] != regioni) + { + regionRegions[regioni].insert + ( + meshEdgeRegions[edgei][edgeRegioni] + ); + } + } + } + } + Pstream::listCombineGather(regionRegions, plusEqOp()); + Pstream::listCombineScatter(regionRegions); + + // Collapse the region connections into a map between each region + // and the lowest numbered region that it connects to + forAll(regionRegions, regioni) + { + forAllConstIter(labelHashSet, regionRegions[regioni], iter) + { + regionMap[iter.key()] = + min(regionMap[iter.key()], regionMap[regioni]); + } + } + } + + // Step 3: Combine connected regions + labelList regionNFaces; + { + // Remove duplicates from the region map + label regioni0 = 0; + forAll(regionMap, regioni) + { + if (regionMap[regioni] > regioni0) + { + ++ regioni0; + regionMap[regioni] = regioni0; + } + } + + // Recompute the number of regions + nRegions = regioni0 + 1; + + // Renumber the face region ID-s + newSetFaceRegions = + IndirectList