From c605a62f76b55765df7183bdabe1c397e6233956 Mon Sep 17 00:00:00 2001
From: Kutalmis Bercin <kutalmis.bercin@esi-group.com>
Date: Thu, 4 Jun 2020 15:40:14 +0100
Subject: [PATCH] INT: integrate planeToFaceZone topoSetSource
---
src/meshTools/Make/files | 1 +
.../planeToFaceZone/planeToFaceZone.C | 413 ++++++++++++++++++
.../planeToFaceZone/planeToFaceZone.H | 208 +++++++++
3 files changed, 622 insertions(+)
create mode 100644 src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.C
create mode 100644 src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.H
diff --git a/src/meshTools/Make/files b/src/meshTools/Make/files
index 1d46037d254..31fe65adff5 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 00000000000..53de3eb39f8
--- /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 <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#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::plainToFaceZone::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, notEqOp<bool>());
+
+ // 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<face>(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>(),
+ labelList()
+ );
+
+ // Combine edge regions to create a list of what regions a given
+ // region is connected to
+ List<labelHashSet> 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<labelHashSet>());
+ 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<label>(regionMap, newSetFaceRegions);
+
+ // Report the final number and size of the regions
+ regionNFaces = labelList(nRegions, 0);
+ forAll(newSetFaces, newSetFacei)
+ {
+ regionNFaces[newSetFaceRegions[newSetFacei]] ++;
+ }
+ Pstream::listCombineGather(regionNFaces, plusEqOp<label>());
+ Pstream::listCombineScatter(regionNFaces);
+ Info<< " Found " << nRegions << " contiguous regions with "
+ << regionNFaces << " faces" << endl;
+ }
+
+ // Step 4: Choose the closest region to output
+ label selectedRegioni = -1;
+ {
+ // Compute the region centres
+ scalarField regionMagAreas(nRegions, 0);
+ pointField regionCentres(nRegions, Zero);
+ forAll(newSetFaces, newSetFacei)
+ {
+ const label facei = newSetFaces[newSetFacei];
+ const label regioni = newSetFaceRegions[newSetFacei];
+
+ const vector& a = mesh_.faceAreas()[facei];
+ const point& c = mesh_.faceCentres()[facei];
+
+ regionMagAreas[regioni] += mag(a);
+ regionCentres[regioni] += mag(a)*c;
+ }
+ Pstream::listCombineGather(regionMagAreas, plusEqOp<scalar>());
+ Pstream::listCombineGather(regionCentres, plusEqOp<point>());
+ Pstream::listCombineScatter(regionMagAreas);
+ Pstream::listCombineScatter(regionCentres);
+ regionCentres /= regionMagAreas;
+
+ // Find the region centroid closest to the reference point
+ selectedRegioni =
+ returnReduce
+ (
+ findMin(mag(regionCentres - point_)()),
+ minOp<label>()
+ );
+
+ // Report the selection
+ Info<< " Selecting region " << selectedRegioni << " with "
+ << regionNFaces[selectedRegioni]
+ << " faces as the closest to point " << point_ << endl;
+ }
+
+ // Step 5: Remove any faces from the set list not in the selected region
+ {
+ // Remove faces from the list by shuffling up and resizing
+ label newSetFacei0 = 0;
+ forAll(newSetFaces, newSetFacei)
+ {
+ newSetFaces[newSetFacei0] = newSetFaces[newSetFacei];
+
+ if (newSetFaceRegions[newSetFacei] == selectedRegioni)
+ {
+ ++ newSetFacei0;
+ }
+ }
+ newSetFaces.resize(newSetFacei0);
+ }
+ }
+
+ // Modify the face zone set
+ DynamicList<label> newAddressing;
+ DynamicList<bool> newFlipMap;
+ if (add)
+ {
+ // Start from copy
+ newAddressing = DynamicList<label>(fzSet.addressing());
+ newFlipMap = DynamicList<bool>(fzSet.flipMap());
+
+ // Add anything from the new set that is not already in the zone set
+ forAll(newSetFaces, newSetFacei)
+ {
+ const label facei = newSetFaces[newSetFacei];
+
+ if (!fzSet.found(facei))
+ {
+ newAddressing.append(facei);
+ newFlipMap.append(cellIsAbovePlane[mesh_.faceOwner()[facei]]);
+ }
+ }
+ }
+ else
+ {
+ // Start from empty
+ newAddressing = DynamicList<label>(fzSet.addressing().size());
+ newFlipMap = DynamicList<bool>(fzSet.flipMap().size());
+
+ // Add everything from the zone set that is not also in the new set
+ labelHashSet newSet(newSetFaces);
+ forAll(fzSet.addressing(), i)
+ {
+ const label facei = fzSet.addressing()[i];
+
+ if (!newSet.found(facei))
+ {
+ newAddressing.append(facei);
+ newFlipMap.append(cellIsAbovePlane[mesh_.faceOwner()[facei]]);
+ }
+ }
+ }
+ fzSet.addressing().transfer(newAddressing);
+ fzSet.flipMap().transfer(newFlipMap);
+ fzSet.updateSet();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::planeToFaceZone::planeToFaceZone
+(
+ const polyMesh& mesh,
+ const dictionary& dict
+)
+:
+ topoSetSource(mesh),
+ point_(dict.get<vector>("point")),
+ normal_(dict.get<vector>("normal")),
+ option_
+ (
+ faceZoneActionNames_.getOrDefault("option", dict, faceZoneAction::ALL)
+ )
+{}
+
+
+Foam::planeToFaceZone::planeToFaceZone
+(
+ const polyMesh& mesh,
+ Istream& is
+)
+:
+ topoSetSource(mesh),
+ point_(checkIs(is)),
+ normal_(checkIs(is)),
+ option_(faceZoneActionNames_.read(checkIs(is)))
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::planeToFaceZone::applyToSet
+(
+ const topoSetSource::setAction action,
+ topoSet& set
+) const
+{
+ if (!isA<faceZoneSet>(set))
+ {
+ WarningInFunction
+ << "Operation only allowed on a faceZoneSet." << endl;
+ }
+ else
+ {
+ faceZoneSet& fzSet = refCast<faceZoneSet>(set);
+
+ if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+ {
+ Info<< " Adding faces which form a plane at " << point_
+ << " with normal " << normal_ << endl;
+
+ combine(fzSet, true);
+ }
+ else if (action == topoSetSource::DELETE)
+ {
+ Info<< " Removing faces which form a plane at " << point_
+ << " with normal " << normal_ << endl;
+
+ combine(fzSet, false);
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.H b/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.H
new file mode 100644
index 00000000000..04ecc718570
--- /dev/null
+++ b/src/meshTools/sets/faceZoneSources/planeToFaceZone/planeToFaceZone.H
@@ -0,0 +1,208 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / 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 <http://www.gnu.org/licenses/>.
+
+Class
+ Foam::planeToFaceZone
+
+Description
+ A \c topoSetSource to select faces based on the adjacent cell centres
+ spanning a given plane. The plane is defined by a point and normal vector.
+
+ Operands:
+ \table
+ Operand | Type | Location
+ output 1 | faceSet | $FOAM_CASE/constant/polyMesh/sets/\<set\>
+ output 2 | faceZone | $FOAM_CASE/constant/polyMesh/faceZones
+ \endtable
+
+Usage
+ Minimal example by using \c system/topoSetDict.actions:
+ \verbatim
+ {
+ // Mandatory (inherited) entries
+ name <name>;
+ type faceZoneSet;
+ action <action>;
+
+ // Mandatory entries
+ source planeToFaceZone;
+ point (<px> <py> <pz>);
+ normal (<nx> <ny> <nz>);
+
+ // Optional entries
+ option <option>;
+ }
+ \endverbatim
+
+ where the entries mean:
+ \table
+ Property | Description | Type | Req'd | Dflt
+ name | Name of faceZone | word | yes | -
+ type | Type name: faceZoneSet | word | yes | -
+ action | Action applied on faces - see below | word | yes | -
+ source | Source name: planeToFaceZone | word | yes | -
+ point | A point on the plane | vector | yes | -
+ normal | The plane normal for orientation | vector | yes | -
+ option | Selection method - see below | word | no | all
+ \endtable
+
+ Options for the \c action entry:
+ \verbatim
+ new | Create a new faceZone from selected faces of a faceZone
+ add | Add selected faces of a faceZone into this faceZone
+ subtract | Remove selected faces of a faceZone from this faceZone
+ \endverbatim
+
+ Options for the \c option entry:
+ \verbatim
+ all |
+ closest |
+ \endverbatim
+
+Note
+ Additionally, an include entry can be specified. When omitted or set to
+ "all", then all faces that meet the criteria are included in the set. When
+ set to "closest", just the faces that belong to the closest contiguous
+ region to the plane point are included. This latter setting is useful when
+ defining face zones through channels on which the flow rate is to be
+ computed, as it keeps the set local to a single channel.
+
+See also
+ - Foam::topoSetSource
+
+SourceFiles
+ planeToFaceZone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef planeToFaceZone_H
+#define planeToFaceZone_H
+
+#include "topoSetSource.H"
+#include "Enum.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+class faceZoneSet;
+
+/*---------------------------------------------------------------------------*\
+ Class planeToFaceZone Declaration
+\*---------------------------------------------------------------------------*/
+
+class planeToFaceZone
+:
+ public topoSetSource
+{
+public:
+
+ //- Enumeration defining the valid options
+ enum faceZoneAction
+ {
+ ALL,
+ CLOSEST
+ };
+
+
+private:
+
+ // Private Data
+
+ //- Add usage string
+ static addToUsageTable usage_;
+
+ //- Included region names
+ //static const NamedEnum<include, 2> includeNames_;
+ static const Enum<faceZoneAction> faceZoneActionNames_;
+
+ //- Point on the plane
+ const vector point_;
+
+ //- Normal to the plane
+ const vector normal_;
+
+ //- Included regions
+ const faceZoneAction option_;
+
+
+ // Private Member Functions
+
+ void combine(faceZoneSet& set, const bool add) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("planeToFaceZone");
+
+
+ // Constructors
+
+ //- Construct from dictionary
+ planeToFaceZone
+ (
+ const polyMesh& mesh,
+ const dictionary& dict
+ );
+
+ //- Construct from Istream
+ planeToFaceZone
+ (
+ const polyMesh& mesh,
+ Istream&
+ );
+
+
+ //- Destructor
+ virtual ~planeToFaceZone() = default;
+
+
+ // Member Functions
+
+ virtual sourceType setType() const
+ {
+ return FACESETSOURCE;
+ }
+
+ virtual void applyToSet
+ (
+ const topoSetSource::setAction action,
+ topoSet&
+ ) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
--
GitLab