diff --git a/src/finiteArea/Make/files b/src/finiteArea/Make/files index 85d56a5e87c452baa3c4648fb7666f9519bdc750..2badfb337dbc4ba0167c44651b5a20094e979729 100644 --- a/src/finiteArea/Make/files +++ b/src/finiteArea/Make/files @@ -4,6 +4,7 @@ faMesh/faMeshDemandDrivenData.C faMesh/faMeshPatches.C faMesh/faMeshTopology.C faMesh/faMeshUpdate.C +faMesh/faMeshBoundaryHalo.C faMesh/faBoundaryMesh/faBoundaryMesh.C faPatches = faMesh/faPatches diff --git a/src/finiteArea/faMesh/faMesh.C b/src/finiteArea/faMesh/faMesh.C index 021b83a3270939527d78974084ca6dd6f6fe0ee8..40763bce48befa868216208c758a06289d23ed24 100644 --- a/src/finiteArea/faMesh/faMesh.C +++ b/src/finiteArea/faMesh/faMesh.C @@ -27,6 +27,7 @@ License \*---------------------------------------------------------------------------*/ #include "faMesh.H" +#include "faMeshBoundaryHalo.H" #include "faGlobalMeshData.H" #include "Time.H" #include "polyMesh.H" @@ -109,6 +110,36 @@ static labelList selectPatchFaces // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // +void Foam::faMesh::checkBoundaryEdgeLabelRange +( + const labelUList& edgeLabels +) const +{ + label nErrors = 0; + + for (const label edgei : edgeLabels) + { + if (edgei < nInternalEdges_ || edgei >= nEdges_) + { + if (!nErrors++) + { + FatalErrorInFunction + << "Boundary edge label out of range " + << nInternalEdges_ << ".." << (nEdges_-1) << nl + << " "; + } + + FatalError<< ' ' << edgei; + } + } + + if (nErrors) + { + FatalError << nl << exit(FatalError); + } +} + + void Foam::faMesh::initPatch() const { patchPtr_.reset @@ -120,6 +151,9 @@ void Foam::faMesh::initPatch() const ) ); bndConnectPtr_.reset(nullptr); + haloMapPtr_.reset(nullptr); + haloFaceCentresPtr_.reset(nullptr); + haloFaceNormalsPtr_.reset(nullptr); } @@ -168,10 +202,21 @@ void Foam::faMesh::setPrimitiveMeshData() } +void Foam::faMesh::clearHalo() const +{ + DebugInFunction << "Clearing halo information" << endl; + + haloMapPtr_.reset(nullptr); + haloFaceCentresPtr_.reset(nullptr); + haloFaceNormalsPtr_.reset(nullptr); +} + + void Foam::faMesh::clearGeomNotAreas() const { DebugInFunction << "Clearing geometry" << endl; + clearHalo(); patchPtr_.reset(nullptr); bndConnectPtr_.reset(nullptr); deleteDemandDrivenData(SPtr_); @@ -274,7 +319,11 @@ Foam::faMesh::faMesh(const polyMesh& pMesh) faceCurvaturesPtr_(nullptr), edgeTransformTensorsPtr_(nullptr), correctPatchPointNormalsPtr_(nullptr), - globalMeshDataPtr_(nullptr) + globalMeshDataPtr_(nullptr), + + haloMapPtr_(nullptr), + haloFaceCentresPtr_(nullptr), + haloFaceNormalsPtr_(nullptr) { DebugInFunction << "Creating from IOobject" << endl; @@ -368,7 +417,11 @@ Foam::faMesh::faMesh faceCurvaturesPtr_(nullptr), edgeTransformTensorsPtr_(nullptr), correctPatchPointNormalsPtr_(nullptr), - globalMeshDataPtr_(nullptr) + globalMeshDataPtr_(nullptr), + + haloMapPtr_(nullptr), + haloFaceCentresPtr_(nullptr), + haloFaceNormalsPtr_(nullptr) {} diff --git a/src/finiteArea/faMesh/faMesh.H b/src/finiteArea/faMesh/faMesh.H index 313db0517aa566dadeccb0d9826a09814dd5dc28..c18674b716c2fa78093f88833b3122495f982634 100644 --- a/src/finiteArea/faMesh/faMesh.H +++ b/src/finiteArea/faMesh/faMesh.H @@ -71,6 +71,7 @@ namespace Foam { // Forward Declarations +class faMeshBoundaryHalo; class faMeshLduAddressing; class faMeshMapper; class faPatchData; @@ -291,11 +292,20 @@ class faMesh mutable boolList* correctPatchPointNormalsPtr_; - // Other mesh-related data + // Other mesh-related data //- Parallel info mutable autoPtr<faGlobalMeshData> globalMeshDataPtr_; + //- Mapping/swapping for boundary edge halo neighbours + mutable std::unique_ptr<faMeshBoundaryHalo> haloMapPtr_; + + //- Face centres for boundary edge halo neighbours + mutable std::unique_ptr<pointField> haloFaceCentresPtr_; + + //- Face normals for boundary edge halo neighbours + mutable std::unique_ptr<vectorField> haloFaceNormalsPtr_; + // Static Private Data @@ -376,6 +386,9 @@ class faMesh //- Clear geometry but not the face areas void clearGeomNotAreas() const; + //- Clear boundary halo information + void clearHalo() const; + //- Clear geometry void clearGeom() const; @@ -386,6 +399,12 @@ class faMesh void clearOut() const; + // Halo handling + + //- Calculate halo centres/normals + void calcHaloFaceGeometry() const; + + // Helpers //- Create a single patch @@ -404,6 +423,32 @@ class faMesh ) const; + //- Fatal error if edge labels are out of range + void checkBoundaryEdgeLabelRange(const labelUList& edgeLabels) const; + + //- Extract list from contiguous (unordered) boundary data + //- to the locally sorted order. + template<class T> + List<T> boundarySubset + ( + const UList<T>& bndField, + const labelUList& edgeLabels + ) const + { + #ifdef FULLDEBUG + checkBoundaryEdgeLabelRange(edgeLabels); + #endif + // Like UIndirectList but with an offset + List<T> result(edgeLabels.size()); + forAll(edgeLabels, i) + { + result[i] = bndField[edgeLabels[i] - nInternalEdges_]; + } + return result; + } + + + public: // Public Typedefs @@ -607,6 +652,21 @@ public: //- (does not include own proc) List<labelPair> boundaryProcSizes() const; + //- Mapping/swapping for boundary halo neighbours + const faMeshBoundaryHalo& boundaryHaloMap() const; + + //- Face centres of boundary halo neighbours + const pointField& haloFaceCentres() const; + + //- Face normals of boundary halo neighbours + const vectorField& haloFaceNormals() const; + + //- Face centres of boundary halo neighbours for specified patch + tmp<pointField> haloFaceCentres(const label patchi) const; + + //- Face normals of boundary halo neighbours for specified patch + tmp<vectorField> haloFaceNormals(const label patchi) const; + // Mesh motion and morphing diff --git a/src/finiteArea/faMesh/faMeshBoundaryHalo.C b/src/finiteArea/faMesh/faMeshBoundaryHalo.C new file mode 100644 index 0000000000000000000000000000000000000000..81aa63c11752cc6f2ef48b326e3104e86bba5e4f --- /dev/null +++ b/src/finiteArea/faMesh/faMeshBoundaryHalo.C @@ -0,0 +1,194 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 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 "faMeshBoundaryHalo.H" +#include "faMesh.H" +#include "globalIndex.H" +#include "Pstream.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(faMeshBoundaryHalo, 0); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::faMeshBoundaryHalo::faMeshBoundaryHalo(const label comm) +: + mapDistributeBase(comm), + inputMeshFaces_(), + boundaryToCompact_() +{} + + +Foam::faMeshBoundaryHalo::faMeshBoundaryHalo(const faMesh& areaMesh) +: + mapDistributeBase(), + inputMeshFaces_(), + boundaryToCompact_() +{ + reset(areaMesh); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +void Foam::faMeshBoundaryHalo::clear() +{ + static_cast<mapDistributeBase&>(*this) = mapDistributeBase(); + + inputMeshFaces_.clear(); + boundaryToCompact_.clear(); +} + + +Foam::label Foam::faMeshBoundaryHalo::haloSize() const +{ + if (Pstream::parRun()) + { + return boundaryToCompact_.size(); + } + else + { + return inputMeshFaces_.size(); + } +} + + +void Foam::faMeshBoundaryHalo::reset(const faMesh& areaMesh) +{ + inputMeshFaces_.clear(); + boundaryToCompact_.clear(); + + const auto& procConnections = areaMesh.boundaryConnections(); + + if (!Pstream::parRun()) + { + // Serial - extract halo numbers directly + + inputMeshFaces_.resize(procConnections.size()); + + forAll(procConnections, connecti) + { + // Connected neighbour, non-parallel = must be local + const auto& tuple = procConnections[connecti]; + // const label nbrProci = tuple.first(); + const label nbrFacei = tuple.second(); + + inputMeshFaces_[connecti] = nbrFacei; + } + + return; + } + + const label nProcs = Pstream::nProcs(comm_); + const label myRank = Pstream::myProcNo(comm_); + + const globalIndex globalFaceNum(areaMesh.mesh().nFaces()); + + // Boundary inside faces in polyMesh face ids + const labelList insideFaces + ( + UIndirectList<label> + ( + areaMesh.faceLabels(), + areaMesh.patch().boundaryFaces() + ) + ); + + + // Slightly circuitous, but allows maximum reuse of mapDistributeBase + + // 1. Construct a connectivity map using global face numbers + + labelListList connectivity(areaMesh.nBoundaryEdges()); + List<Map<label>> compactMap(nProcs, Map<label>(0)); + + // All local mesh faces used + labelHashSet localUsed(insideFaces); + + forAll(connectivity, connecti) + { + labelList& edgeFaces = connectivity[connecti]; + edgeFaces.resize(2); + + // Owner is the boundary inside face (our side) + // Neighbour is the boundary outside face + + // Connected neighbour + const auto& tuple = procConnections[connecti]; + const label nbrProci = tuple.first(); + const label nbrFacei = tuple.second(); + + if (myRank == nbrProci) + { + // Processor-local connectivity + localUsed.insert(nbrFacei); + } + + // Global addressing for the connectivity + edgeFaces[0] = globalFaceNum.toGlobal(insideFaces[connecti]); + edgeFaces[1] = globalFaceNum.toGlobal(nbrProci, nbrFacei); + } + + // Create and replace mapping + static_cast<mapDistributeBase&>(*this) = mapDistributeBase + ( + globalFaceNum, + connectivity, + compactMap, + Pstream::msgType(), + comm_ + ); + + // List of local mesh faces referenced. + // Includes inside and locally connected outside faces + + inputMeshFaces_ = localUsed.sortedToc(); + + boundaryToCompact_.clear(); + boundaryToCompact_.resize(connectivity.size()); + + // After creating the map, connectivity is localized *and* + // uses compact numbering! + + // Extract the neighbour connection (compact numbering) + forAll(connectivity, connecti) + { + const labelList& edgeFaces = connectivity[connecti]; + // const label face0 = edgeFaces[0]; + const label face1 = edgeFaces[1]; + + boundaryToCompact_[connecti] = face1; + } +} + + +// ************************************************************************* // diff --git a/src/finiteArea/faMesh/faMeshBoundaryHalo.H b/src/finiteArea/faMesh/faMeshBoundaryHalo.H new file mode 100644 index 0000000000000000000000000000000000000000..4a8f58b33e69b8adbf416bb53c9062eba710058c --- /dev/null +++ b/src/finiteArea/faMesh/faMeshBoundaryHalo.H @@ -0,0 +1,154 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 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::faMeshBoundaryHalo + +Description + Class for obtaining halo face data for the boundary edges. + The ordering follows that natural edge ordering of the underlying + primitive patch. + +Note + The halo faces can be located on-processor or off-processor. + +SourceFiles + faMeshBoundaryHalo.C + faMeshBoundaryHaloTemplates.C + +\*---------------------------------------------------------------------------*/ + +#ifndef faMeshBoundaryHalo_H +#define faMeshBoundaryHalo_H + +#include "mapDistributeBase.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// Forward Declarations +class faMesh; + +/*---------------------------------------------------------------------------*\ + Class faMeshBoundaryHalo Declaration +\*---------------------------------------------------------------------------*/ + +class faMeshBoundaryHalo +: + public mapDistributeBase +{ + // Private Data + + //- List of local input mesh faces required + labelList inputMeshFaces_; + + //- Internal mapping from boundary index to compact + labelList boundaryToCompact_; + + +public: + + // Declare name of the class and its debug switch + ClassName("faMeshBoundaryHalo"); + + + // Constructors + + //- Default construct + explicit faMeshBoundaryHalo(const label comm = UPstream::worldComm); + + //- Construct from mesh + explicit faMeshBoundaryHalo(const faMesh& mesh); + + + // Member Functions + + //- Clear out all parameters + void clear(); + + //- Redefine map and connectivity for a mesh + void reset(const faMesh& mesh); + + //- The local data size (output) + label haloSize() const; + + //- List of local input mesh faces required. + // \note will not correspond exactly to the boundary inside faces. + // Duplicates have been removed and it also contains the + // processor-local neighbour faces, which would otherwise not be + // handled by the distribute method. + const labelList& inputMeshFaces() const noexcept + { + return inputMeshFaces_; + } + + + // Other + + //- Distribute sparse data. + // On output it is adjusted. + template<class Type> + void distributeSparse + ( + List<Type>& fld, + const labelUList& sparseInputLocations, + const labelUList& compactOutputMapping + ) const; + + //- Distribute sparse data. + // On output it is adjusted. + template<class Type> + void distributeSparse + ( + List<Type>& fld, + const labelUList& sparseInputLocations + ) const; + + //- Distribute sparse data. + // The input field one enty per sparse id (inputMeshFaces). + // On output it will have for the input sparse + // The input field contains location. + template<class Type> + void distributeSparse(List<Type>& fld) const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository + #include "faMeshBoundaryHaloTemplates.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/finiteArea/faMesh/faMeshBoundaryHaloTemplates.C b/src/finiteArea/faMesh/faMeshBoundaryHaloTemplates.C new file mode 100644 index 0000000000000000000000000000000000000000..69f2214b04fe0e94ad7c762b0d02c2c110b77e60 --- /dev/null +++ b/src/finiteArea/faMesh/faMeshBoundaryHaloTemplates.C @@ -0,0 +1,107 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | www.openfoam.com + \\/ M anipulation | +------------------------------------------------------------------------------- + Copyright (C) 2021 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 "faMeshBoundaryHalo.H" +#include "UIndirectList.H" +#include "Pstream.H" + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template<class Type> +void Foam::faMeshBoundaryHalo::distributeSparse +( + List<Type>& fld, + const labelUList& sparseInputLocations, + const labelUList& compactOutputMapping +) const +{ + if (!Pstream::parRun()) + { + return; // No-op in serial + } + + // Construct data in compact addressing + List<Type> compactFld(constructSize_, Zero); + + if (sparseInputLocations.empty()) + { + // Copy in as dense field + forAll(fld, i) + { + compactFld[i] = fld[i]; + } + } + else + { + if (fld.size() != sparseInputLocations.size()) + { + FatalErrorInFunction + << "Input field size (" << fld.size() + << " != sparse ids size (" + << sparseInputLocations.size() << ")\n" + << exit(FatalError); + } + + // Copy in sparse locations + forAll(sparseInputLocations, i) + { + const label idx = sparseInputLocations[i]; + if (idx != -1) + { + compactFld[idx] = fld[i]; + } + } + } + + + // Pull all data + mapDistributeBase::distribute<Type>(compactFld); + + // Rewrite to output + fld = UIndirectList<Type>(compactFld, compactOutputMapping); +} + + +template<class Type> +void Foam::faMeshBoundaryHalo::distributeSparse +( + List<Type>& fld, + const labelUList& sparseInputLocations +) const +{ + this->distributeSparse(fld, sparseInputLocations, boundaryToCompact_); +} + + +template<class Type> +void Foam::faMeshBoundaryHalo::distributeSparse(List<Type>& fld) const +{ + this->distributeSparse(fld, inputMeshFaces_, boundaryToCompact_); +} + + +// ************************************************************************* // diff --git a/src/finiteArea/faMesh/faMeshTopology.C b/src/finiteArea/faMesh/faMeshTopology.C index 179054ac0b6066b63a1b0240bddd889abf1dab6b..938dbe9dbdc4f18b60abbb48baa462cd83497462 100644 --- a/src/finiteArea/faMesh/faMeshTopology.C +++ b/src/finiteArea/faMesh/faMeshTopology.C @@ -26,6 +26,7 @@ License \*---------------------------------------------------------------------------*/ #include "faMesh.H" +#include "faMeshBoundaryHalo.H" #include "globalMeshData.H" #include "indirectPrimitivePatch.H" #include "edgeHashes.H" @@ -800,4 +801,130 @@ Foam::List<Foam::labelPair> Foam::faMesh::boundaryProcSizes() const } +const Foam::faMeshBoundaryHalo& Foam::faMesh::boundaryHaloMap() const +{ + if (!haloMapPtr_) + { + haloMapPtr_.reset(new faMeshBoundaryHalo(*this)); + } + + return *haloMapPtr_; +} + + +void Foam::faMesh::calcHaloFaceGeometry() const +{ + if (haloFaceCentresPtr_ || haloFaceNormalsPtr_) + { + FatalErrorInFunction + << "Halo centres/normals already calculated" + << exit(FatalError); + } + + DebugInFunction + << "Calculating halo face centres/normals" << endl; + + const faceList& faces = mesh().faces(); + const pointField& points = mesh().points(); + + const faMeshBoundaryHalo& halo = boundaryHaloMap(); + + const labelList& inputFaceIds = halo.inputMeshFaces(); + + haloFaceCentresPtr_.reset(new pointField); + haloFaceNormalsPtr_.reset(new vectorField); + + auto& centres = *haloFaceCentresPtr_; + auto& normals = *haloFaceNormalsPtr_; + + centres.resize(inputFaceIds.size()); + normals.resize(inputFaceIds.size()); + + // My values + forAll(inputFaceIds, i) + { + const face& f = faces[inputFaceIds[i]]; + + centres[i] = f.centre(points); + normals[i] = f.unitNormal(points); + } + + // Swap information and resize + halo.distributeSparse(centres); + halo.distributeSparse(normals); +} + + +const Foam::pointField& Foam::faMesh::haloFaceCentres() const +{ + if (!haloFaceCentresPtr_ || !haloFaceNormalsPtr_) + { + calcHaloFaceGeometry(); + } + + return *haloFaceCentresPtr_; +} + + +const Foam::vectorField& Foam::faMesh::haloFaceNormals() const +{ + if (!haloFaceCentresPtr_ || !haloFaceNormalsPtr_) + { + calcHaloFaceGeometry(); + } + + return *haloFaceNormalsPtr_; +} + + +Foam::tmp<Foam::pointField> +Foam::faMesh::haloFaceCentres(const label patchi) const +{ + if (patchi < 0 || patchi >= boundary().size()) + { + FatalErrorInFunction + << "Patch " << patchi << " is out-of-range 0.." + << (boundary().size()-1) << nl + << exit(FatalError); + } + + return tmp<pointField>::New + ( + List<point> + ( + boundarySubset + ( + this->haloFaceCentres(), + boundary()[patchi].edgeLabels() + ) + ) + ); +} + + +Foam::tmp<Foam::vectorField> +Foam::faMesh::haloFaceNormals(const label patchi) const +{ + if (patchi < 0 || patchi >= boundary().size()) + { + FatalErrorInFunction + << "Patch " << patchi << " is out-of-range 0.." + << (boundary().size()-1) << nl + << exit(FatalError); + } + + return tmp<vectorField>::New + ( + List<vector> + ( + boundarySubset + ( + this->haloFaceNormals(), + boundary()[patchi].edgeLabels() + ) + ) + ); +} + + // ************************************************************************* //