diff --git a/src/overset/oversetMesh/Make/files b/src/overset/oversetMesh/Make/files index 1f27b60ff..d21982f92 100644 --- a/src/overset/oversetMesh/Make/files +++ b/src/overset/oversetMesh/Make/files @@ -4,6 +4,8 @@ oversetFringe/oversetFringe/oversetFringe.C oversetFringe/oversetFringe/newOversetFringe.C oversetFringe/manualFringe/manualFringe.C oversetFringe/faceCellsFringe/faceCellsFringe.C +oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.C +oversetFringe/cuttingPatchFringe/cuttingPatchFringe.C oversetFringe/overlapFringe/overlapFringe/overlapFringe.C oversetFringe/overlapFringe/layeredOverlapFringe/layeredOverlapFringe.C oversetFringe/overlapFringe/adaptiveOverlapFringe/adaptiveOverlapFringe.C diff --git a/src/overset/oversetMesh/Make/options b/src/overset/oversetMesh/Make/options index d07b7e112..9444c1444 100644 --- a/src/overset/oversetMesh/Make/options +++ b/src/overset/oversetMesh/Make/options @@ -2,12 +2,10 @@ EXE_INC = \ -I$(LIB_SRC)/finiteVolume/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ -I$(LIB_SRC)/surfMesh/lnInclude \ - -I$(LIB_SRC)/sampling/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude + -I$(LIB_SRC)/sampling/lnInclude EXE_LIBS = \ -lfiniteVolume \ -lmeshTools \ -lsurfMesh \ - -lsampling \ - -ldynamicMesh + -lsampling diff --git a/src/overset/oversetMesh/oversetAdjustPhi/oversetAdjustPhi.C b/src/overset/oversetMesh/oversetAdjustPhi/oversetAdjustPhi.C index 31086ad42..2f40298be 100644 --- a/src/overset/oversetMesh/oversetAdjustPhi/oversetAdjustPhi.C +++ b/src/overset/oversetMesh/oversetAdjustPhi/oversetAdjustPhi.C @@ -59,7 +59,7 @@ void Foam::oversetAdjustPhi const unallocLabelList& neighbour = mesh.neighbour(); // Get region split to identify separate mesh components - const labelList& regionID = om.regionID(); + const scalarField& regionID = om.regionID().internalField(); // Sum up incoming and outgoing flux scalar fringeIn = 0; @@ -80,8 +80,12 @@ void Foam::oversetAdjustPhi { // Internal face - // Debug check - if (regionID[owner[curFace]] != regionID[neighbour[curFace]]) + // Debug check. Note: use notEqual since we are comparing two + // scalars + if + ( + notEqual(regionID[owner[curFace]], regionID[neighbour[curFace]]) + ) { FatalErrorIn ( diff --git a/src/overset/oversetMesh/oversetAdjustPhi/regionWiseOversetAdjustPhi.C b/src/overset/oversetMesh/oversetAdjustPhi/regionWiseOversetAdjustPhi.C index 21d25cb47..73e682e49 100644 --- a/src/overset/oversetMesh/oversetAdjustPhi/regionWiseOversetAdjustPhi.C +++ b/src/overset/oversetMesh/oversetAdjustPhi/regionWiseOversetAdjustPhi.C @@ -60,7 +60,7 @@ void Foam::regionWiseOversetAdjustPhi const unallocLabelList& neighbour = mesh.neighbour(); // Get region split to identify separate mesh components - const labelList& regionID = om.regionID(); + const scalarField& regionID = om.regionID().internalField(); // Incoming and outgoing region fluxes scalarField regionIn(om.regions().size(), 0); @@ -94,7 +94,7 @@ void Foam::regionWiseOversetAdjustPhi // scaled forAll (phip, i) { - // Get current region index + // Get current region index. Note: conversion to label const label curRegion = regionID[fc[i]]; if (phip[i] < 0.0) @@ -124,11 +124,12 @@ void Foam::regionWiseOversetAdjustPhi { // Internal face - // Get region index + // Get region index. Note: conversion to label const label curRegion = regionID[owner[curFace]]; - // Check whether owner and neighbour belong to the same region - if (curRegion != regionID[neighbour[curFace]]) + // Check whether owner and neighbour belong to the same region. + // Note: use notEqual function since we are comparing two scalars + if (notEqual(curRegion, regionID[neighbour[curFace]])) { FatalErrorIn ( @@ -210,7 +211,7 @@ void Foam::regionWiseOversetAdjustPhi { // Processor patch, master side - // Get region index + // Get region index. Note: conversion to label const label curRegion = regionID[mesh.boundary()[patchI].faceCells()[faceI]]; @@ -321,7 +322,7 @@ void Foam::regionWiseOversetAdjustPhi { // Internal face - // Get region index + // Get region index. Note: conversion to label const label curRegion = regionID[owner[curFace]]; // Get reference to the flux for scaling @@ -358,7 +359,7 @@ void Foam::regionWiseOversetAdjustPhi if (procPatch.master()) // Owner side { - // Get region index + // Get region index. Note: conversion to label const label curRegion = regionID[mesh.boundary()[patchI].faceCells()[faceI]]; @@ -374,7 +375,7 @@ void Foam::regionWiseOversetAdjustPhi } else // Neighbouring processor side { - // Get region index + // Get region index. Note: conversion to label const label curRegion = regionID[mesh.boundary()[patchI].faceCells()[faceI]]; diff --git a/src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.C b/src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.C new file mode 100644 index 000000000..bf4bf01f6 --- /dev/null +++ b/src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.C @@ -0,0 +1,649 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.1 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "cuttingPatchFringe.H" +#include "oversetMesh.H" +#include "oversetRegion.H" +#include "polyPatchID.H" +#include "addToRunTimeSelectionTable.H" +#include "syncTools.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(cuttingPatchFringe, 0); + addToRunTimeSelectionTable + ( + oversetFringe, + cuttingPatchFringe, + dictionary + ); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::cuttingPatchFringe::calcAddressing() const +{ + // Make sure that either acceptorsPtr is unnalocated or if it is allocated, + // that it is empty + if (acceptorsPtr_ && !acceptorsPtr_->empty()) + { + FatalErrorIn + ( + "void Foam::cuttingPatchFringe::calcAddressing() const" + ) << "Addressing already calculated" + << abort(FatalError); + } + + // Get polyMesh + const polyMesh& mesh = this->mesh(); + + // Collect all cutting patches + labelHashSet patchIDs(cuttingPatchNames_.size()); + + forAll (cuttingPatchNames_, nameI) + { + // Get polyPatchID and check if valid + const polyPatchID cutPatch + ( + cuttingPatchNames_[nameI], + mesh.boundaryMesh() + ); + + if (!cutPatch.active()) + { + FatalErrorIn + ( + "void cuttingPatchFringe::calcAddressing const" + ) << "Cutting patch " << cuttingPatchNames_[nameI] + << " cannot be found." + << abort(FatalError); + } + + // Store patch ID in the set + patchIDs.insert(cutPatch.index()); + } + + if (debug) + { + Info<< "Starting cutting patch fringe assembly..." << endl; + } + + // Note: similar code as in oversetRegion::calcHoleTriMesh. Consider + // refactoring. VV, 20/May/2019 + + // Make and invert local triSurface + triFaceList triFaces; + pointField triPoints; + + // Memory management + { + triSurface ts = triSurfaceTools::triangulate + ( + mesh.boundaryMesh(), + patchIDs + ); + + // Clean mutiple points and zero-sized triangles + ts.cleanup(false); + + triFaces.setSize(ts.size()); + triPoints = ts.points(); + + forAll (ts, tsI) + { + // Bugfix: no need to reverse the face because the normals point in + // the correct direction already. VV, 20/May/2019. + triFaces[tsI] = ts[tsI]; + } + } + + if (Pstream::parRun()) + { + // Combine all faces and points into a single list + + List allTriFaces(Pstream::nProcs()); + List allTriPoints(Pstream::nProcs()); + + allTriFaces[Pstream::myProcNo()] = triFaces; + allTriPoints[Pstream::myProcNo()] = triPoints; + + Pstream::gatherList(allTriFaces); + Pstream::scatterList(allTriFaces); + + Pstream::gatherList(allTriPoints); + Pstream::scatterList(allTriPoints); + + // Re-pack points and faces + + label nTris = 0; + label nPoints = 0; + + forAll (allTriFaces, procI) + { + nTris += allTriFaces[procI].size(); + nPoints += allTriPoints[procI].size(); + } + + // Pack points + triPoints.setSize(nPoints); + + // Prepare point renumbering array + labelListList renumberPoints(Pstream::nProcs()); + + nPoints = 0; + + forAll (allTriPoints, procI) + { + const pointField& ptp = allTriPoints[procI]; + + renumberPoints[procI].setSize(ptp.size()); + + labelList& procRenumberPoints = renumberPoints[procI]; + + forAll (ptp, ptpI) + { + triPoints[nPoints] = ptp[ptpI]; + procRenumberPoints[ptpI] = nPoints; + + nPoints++; + } + } + + // Pack triangles and renumber into complete points on the fly + triFaces.setSize(nTris); + + nTris = 0; + + forAll (allTriFaces, procI) + { + const triFaceList& ptf = allTriFaces[procI]; + + const labelList& procRenumberPoints = renumberPoints[procI]; + + forAll (ptf, ptfI) + { + const triFace& procFace = ptf[ptfI]; + + triFace& renumberFace = triFaces[nTris]; + + forAll (renumberFace, rfI) + { + renumberFace[rfI] = + procRenumberPoints[procFace[rfI]]; + } + + nTris++; + } + } + } + + // Make a complete triSurface from local data + triSurface patchTriMesh(triFaces, triPoints); + + // Clean up duplicate points and zero sized triangles + patchTriMesh.cleanup(false); + + // Get this region + const oversetRegion& myRegion = this->region(); + + // Debug: write down the tri mesh + if (Pstream::master()) + { + patchTriMesh.write + ( + word + ( + "patchTriSurface_region" + myRegion.name() + ".vtk" + ) + ); + } + + // Create the tri surface search object + triSurfaceSearch patchTriSearch(patchTriMesh); + + // Get cells in this region + const labelList& myRegionCells = myRegion.regionCells(); + + // Get cell centres for inside-outside search using search object + vectorField myCC(myRegionCells.size()); + + // Cell centres from polyMesh + const vectorField& cc = mesh.cellCentres(); + + forAll (myCC, i) + { + myCC[i] = cc[myRegionCells[i]]; + } + + // Inside mask: all cells within search object will be marked + boolList insideMask(mesh.nCells(), false); + + // Get inside cells for cells in my region only + boolList myRegionInsideMask = patchTriSearch.calcInside(myCC); + + // Note: insideMask has the size of all mesh cells and + // myRegionInsideMask has the size of cells in this region + forAll (myRegionInsideMask, i) + { + insideMask[myRegionCells[i]] = myRegionInsideMask[i]; + } + + // Make sure that the cut holes for this region are also properly marked as + // "inside". This may not be the case automatically for e.g. simulations + // with appendages + const labelList& cutRegionHoles = myRegion.cutHoles(); + forAll (cutRegionHoles, i) + { + insideMask[cutRegionHoles[i]] = true; + } + + // Get necessary mesh data (from polyMesh/primitiveMesh) + const cellList& meshCells = mesh.cells(); + const unallocLabelList& owner = mesh.faceOwner(); + const unallocLabelList& neighbour = mesh.faceNeighbour(); + + // Bool list for collecting faces with at least one unmarked + // cell (to determine the acceptors for the first iteration) + boolList hasUnmarkedCell(mesh.nFaces(), false); + + // Loop through all cells + forAll (insideMask, cellI) + { + if (!insideMask[cellI]) + { + // This cell is not inside (it is unmarked). Loop through + // its faces and set the flag + const cell& cFaces = meshCells[cellI]; + + forAll (cFaces, i) + { + // Set the mark for this global face + hasUnmarkedCell[cFaces[i]] = true; + } + } + } + + // Sync the face list across processor boundaries + syncTools::syncFaceList + ( + mesh, + hasUnmarkedCell, + orEqOp(), + true + ); + + // Mark-up for all inside faces + boolList insideFaceMask(mesh.nFaces(), false); + + // Collect all acceptors for the first iteration (the cells that + // have at least one neighbour cell that is not marked) + labelHashSet acceptors(myRegionCells.size()/10); + + // Loop again through all cells and collect marked ones into + // acceptors or holes, depending on whether they have unmarked cell + // as a neighbour (indicating an acceptor) + forAll (insideMask, cellI) + { + if (insideMask[cellI]) + { + // This cell is inside the covered region + const cell& cFaces = meshCells[cellI]; + + forAll (cFaces, i) + { + // Get global face index + const label& faceI = cFaces[i]; + + // Check whether this neighbour is unmarked + if (hasUnmarkedCell[faceI]) + { + // This cell has unmarked neighbour, collect it into + // the acceptor list + acceptors.insert(cellI); + + // This cell is no longer "inside cell" + insideMask[cellI] = false; + } + else + { + // This is an "inside" face, mark it + insideFaceMask[faceI] = true; + } + } // End for all faces + } // End if cell is inside + } // End for all cells + + // Note: insideFaceMask already synced across processors because it relies + // on hasUnmarkedCell list, which has been synced just above + + // Hash set containing new acceptors (for successive iterations) + labelHashSet newAcceptors(acceptors.size()); + + // Now that we have the initial set of acceptors (and holes), loop + // nLayers away from initial donors + for (label i = 0; i < nLayers_; ++i) + { + // Face markup for propagation + boolList propagateFace(mesh.nFaces(), false); + + // Loop through all acceptors and mark faces that are around hole + // cells. This way, we make sure that we go towards the correct, + // inside direction + forAllConstIter (labelHashSet, acceptors, iter) + { + // Get the cell index and the cell + const label& cellI = iter.key(); + const cell& cFaces = meshCells[cellI]; + + // Loop through all faces of the cell + forAll (cFaces, i) + { + // Get face index (global) + const label& faceI = cFaces[i]; + + if (insideFaceMask[faceI]) + { + // This is a hole face, we are moving in the right + // direction. Mark the face for propagation + propagateFace[faceI] = true; + } + } // End for all faces of the cell + } // End for all donor cells + + // Sync the face list across processor boundaries + syncTools::syncFaceList + ( + mesh, + propagateFace, + orEqOp(), + false + ); + + // Loop through all internal faces and append acceptors + for (label faceI = 0; faceI < mesh.nInternalFaces(); ++faceI) + { + if (propagateFace[faceI]) + { + // Face is marked, select owner or neighbour + const label& own = owner[faceI]; + const label& nei = neighbour[faceI]; + + // Either owner or neighbour may be hole, not both + if (insideMask[own]) + { + // Owner cell is a hole, insert it + newAcceptors.insert(own); + + // Update hole mask + insideMask[own] = false; + } + else if (insideMask[nei]) + { + // Neighbour cell is a hole, insert it + newAcceptors.insert(nei); + + // Update hole mask + insideMask[nei] = false; + } + + // Also update hole face mask for next iteration + insideFaceMask[faceI] = false; + } + } + + // Loop through boundary faces + for + ( + label faceI = mesh.nInternalFaces(); + faceI < mesh.nFaces(); + ++faceI + ) + { + if (propagateFace[faceI]) + { + // Face is marked, select owner if this is the right + // side. Neighbour handled on the other side + const label& own = owner[faceI]; + + if (insideMask[own]) + { + // Face cell is a hole, insert it + newAcceptors.insert(own); + + // Update hole mask + insideMask[own] = false; + } + + // Also update hole face mask for next iteration + insideFaceMask[faceI] = false; + } + } + + // Transfer newAcceptors into acceptors for next iteration or + // for final assembly. Resize newAcceptors accordingly + acceptors.transfer(newAcceptors); + newAcceptors.resize(acceptors.size()); + + } // End for specified number of layers + + // At this point, we have the final set of acceptors and we marked + // all cells that should be holes. Collect them into the list + dynamicLabelList fringeHoles(myRegionCells.size()/10); + + forAll (insideMask, cellI) + { + if (insideMask[cellI]) + { + fringeHoles.append(cellI); + } + } + + // Set acceptors and holes from the data for all regions + acceptorsPtr_ = new labelList(acceptors.sortedToc()); + fringeHolesPtr_ = new labelList(fringeHoles.xfer()); + + if (debug) + { + Info<< "In cuttingPatchFringe::calcAddressing() const" << nl + << "Found " << acceptorsPtr_->size() << " acceptors." << nl + << "Found " << fringeHolesPtr_->size() << " fringe holes." << endl; + } +} + + +void Foam::cuttingPatchFringe::clearAddressing() const +{ + deleteDemandDrivenData(fringeHolesPtr_); + deleteDemandDrivenData(acceptorsPtr_); + deleteDemandDrivenData(finalDonorAcceptorsPtr_); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +// Construct from dictionary +Foam::cuttingPatchFringe::cuttingPatchFringe +( + const fvMesh& mesh, + const oversetRegion& region, + const dictionary& dict +) +: + oversetFringe(mesh, region, dict), + cuttingPatchNames_(dict.lookup("cuttingPatches")), + nLayers_(readLabel(dict.lookup("nLayers"))), + fringeHolesPtr_(nullptr), + acceptorsPtr_(nullptr), + finalDonorAcceptorsPtr_(nullptr) +{ + // Sanity check number of layers: must be greater than 0 + if (nLayers_ < 1) + { + FatalIOErrorIn + ( + "cuttingPatchFringe::" + "cuttingPatchFringe\n" + "(\n" + " const fvMesh& mesh,\n" + " const oversetRegion& region,\n" + " const dictionary& dict\n" + ")", + dict + ) << "Invalid number of layers specified, nLayers = " << nLayers_ + << nl + << "The number should be greater than 0." + << exit(FatalError); + } + + // Preferably, the number of layers should be at least 2 + if (nLayers_ == 1) + { + WarningIn + ( + "cuttingPatchFringe::" + "cuttingPatchFringe\n" + "(\n" + " const fvMesh& mesh,\n" + " const oversetRegion& region,\n" + " const dictionary& dict\n" + ")" + ) << "You have specified nLayers = " << nLayers_ + << nl + << "We strongly advise to use at least 2 layers to avoid" << nl + << "possibility of having acceptors that cannot find decent" << nl + << "donors on the other side." + << endl; + } +} + + +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // + +Foam::cuttingPatchFringe::~cuttingPatchFringe() +{ + clearAddressing(); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +bool Foam::cuttingPatchFringe::updateIteration +( + donorAcceptorList& donorAcceptorRegionData +) const +{ + // If the donorAcceptor list has been allocated, something went wrong with + // the iteration procedure (not-updated flag): this function has been called + // more than once, which should not happen for cuttingPatchFringe + if (finalDonorAcceptorsPtr_) + { + FatalErrorIn + ( + "cuttingPatchFringe::" + "updateIteration(donorAcceptorList&) const" + ) << "finalDonorAcceptorPtr_ already allocated. Something went " + << "wrong with the iteration procedure (flag was not updated)." + << nl + << "This should not happen for cuttingPatchFringe." + << abort(FatalError); + } + + // Allocate the list by reusing the argument list + finalDonorAcceptorsPtr_ = new donorAcceptorList + ( + donorAcceptorRegionData, + true + ); + + // Set the flag to true and return + updateSuitableOverlapFlag(true); + + return foundSuitableOverlap(); +} + + +const Foam::labelList& Foam::cuttingPatchFringe::fringeHoles() const +{ + if (!fringeHolesPtr_) + { + calcAddressing(); + } + + return *fringeHolesPtr_; +} + + +const Foam::labelList& Foam::cuttingPatchFringe::candidateAcceptors() const +{ + if (!acceptorsPtr_) + { + calcAddressing(); + } + + return *acceptorsPtr_; +} + + +Foam::donorAcceptorList& +Foam::cuttingPatchFringe::finalDonorAcceptors() const +{ + if (!finalDonorAcceptorsPtr_) + { + FatalErrorIn("cuttingPatchFringe::finalDonorAcceptors()") + << "finalDonorAcceptorPtr_ not allocated. Make sure you have" + << " called cuttingPatchFringe::updateIteration() before" + << " asking for final set of donor/acceptor pairs." + << abort(FatalError); + } + + if (!foundSuitableOverlap()) + { + FatalErrorIn("cuttingPatchFringe::finalDonorAcceptors()") + << "Attemted to access finalDonorAcceptors but suitable overlap " + << "has not been found. This is not allowed. " + << abort(FatalError); + } + + return *finalDonorAcceptorsPtr_; +} + + +void Foam::cuttingPatchFringe::update() const +{ + Info<< "cuttingPatchFringe::update() const" << endl; + + // Clear out + clearAddressing(); + + // Set flag to false and clear final donor/acceptors only + deleteDemandDrivenData(finalDonorAcceptorsPtr_); + updateSuitableOverlapFlag(false); +} + + +// ************************************************************************* // diff --git a/src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.H b/src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.H new file mode 100644 index 000000000..7feb6cbe2 --- /dev/null +++ b/src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.H @@ -0,0 +1,151 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.1 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +Class + cuttingPatchFringe + +Description + Fringe algorithm based on collecting all the cells cut by its donor region + patch (assuming faceCells algorithm is used). Inside cells are marked and + moved nLayers towards the interior to define the acceptors as robustly as + possible. + +Author + Vuko Vukcevic, Wikki Ltd. All rights reserved. + +SourceFiles + cuttingPatchFringe.C + +\*---------------------------------------------------------------------------*/ + +#ifndef cuttingPatchFringe_H +#define cuttingPatchFringe_H + +#include "oversetFringe.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class cuttingPatchFringe Declaration +\*---------------------------------------------------------------------------*/ + +class cuttingPatchFringe +: + public oversetFringe +{ + // Private data + + //- Names of the cuttingPatches + wordList cuttingPatchNames_; + + //- How many layers to move away from connected region donors to define + // acceptor (and holes) + label nLayers_; + + //- Fringe hole cells + mutable labelList* fringeHolesPtr_; + + //- Acceptor cells + mutable labelList* acceptorsPtr_; + + //- Final donor/acceptor pairs for this region (fringe) + mutable donorAcceptorList* finalDonorAcceptorsPtr_; + + + // Private Member Functions + + //- Calculate hole and acceptor addressing + void calcAddressing() const; + + //- Clear addressing + void clearAddressing() const; + + +public: + + //- Runtime type information + TypeName("cuttingPatch"); + + + // Constructors + + //- Construct from dictionary + cuttingPatchFringe + ( + const fvMesh& mesh, + const oversetRegion& region, + const dictionary& dict + ); + + //- Disallow default bitwise copy construct + cuttingPatchFringe + ( + const cuttingPatchFringe& + ) = delete; + + //- Disallow default bitwise assignment + void operator=(const cuttingPatchFringe&) = delete; + + + // Destructor + + virtual ~cuttingPatchFringe(); + + + // Member Functions + + //- Update iteration. Note: invalidates parameter + virtual bool updateIteration + ( + donorAcceptorList& donorAcceptorRegionData + ) const; + + //- Return list of deactivated (hole) cells + // Fringe hole cells are collected in addition to geometric hole + // cells, which fall outside of all donor regions + virtual const labelList& fringeHoles() const; + + //- Return list of acceptor cells + virtual const labelList& candidateAcceptors() const; + + //- Return list of final donor acceptor pairs. Note: caller may + // invalidate finalDonorAcceptorsPtr_ for optimisation purposes + virtual donorAcceptorList& finalDonorAcceptors() const; + + //- Update the fringe + virtual void update() const; +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/overset/oversetMesh/oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.C b/src/overset/oversetMesh/oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.C new file mode 100644 index 000000000..b075e9999 --- /dev/null +++ b/src/overset/oversetMesh/oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.C @@ -0,0 +1,887 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | foam-extend: Open Source CFD + \\ / O peration | Version: 4.1 + \\ / A nd | Web: http://www.foam-extend.org + \\/ M anipulation | For copyright notice see file Copyright +------------------------------------------------------------------------------- +License + This file is part of foam-extend. + + foam-extend 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. + + foam-extend 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 foam-extend. If not, see . + +\*---------------------------------------------------------------------------*/ + +#include "donorBasedLayeredOverlapFringe.H" +#include "oversetMesh.H" +#include "oversetRegion.H" +#include "faceCellsFringe.H" +#include "oversetRegion.H" +#include "addToRunTimeSelectionTable.H" +#include "syncTools.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(donorBasedLayeredOverlapFringe, 0); + addToRunTimeSelectionTable + ( + oversetFringe, + donorBasedLayeredOverlapFringe, + dictionary + ); +} + + +const Foam::debug::tolerancesSwitch +Foam::donorBasedLayeredOverlapFringe::distTol_ +( + "donorBasedLayeredOverlapDistanceTolerance", + 0.0 +); + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::donorBasedLayeredOverlapFringe::init() const +{ + // Set size of the list containing IDs + connectedRegionIDs_.setSize(connectedRegionNames_.size()); + + // Get list of all overset regions + const PtrList& allRegions = + this->region().overset().regions(); + + // Create list of all region names for easy lookup + wordList allRegionNames(allRegions.size()); + forAll (allRegionNames, arI) + { + allRegionNames[arI] = allRegions[arI].name(); + } + + // Loop through all regions and check whether the overlap has been found + forAll (connectedRegionNames_, crI) + { + // Get name of this connected region + const word& crName = connectedRegionNames_[crI]; + + // Find this region in the list of all regions + const label regionID = findIndex(allRegionNames, crName); + + if (regionID == -1) + { + FatalErrorIn("void donorBasedLayeredOverlapFringe::init() const") + << "Region " << crName << " not found in list of regions." + << "List of overset regions: " << allRegionNames + << abort(FatalError); + } + + // Collect the region index in the list + connectedRegionIDs_[crI] = regionID; + + // Sanity check: if the specified connected donor region has more than 1 + // donor regions, this fringe algorithm is attempted to be used for + // something that's not intended. Issue an error + if (allRegions[regionID].donorRegions().size() != 1) + { + FatalErrorIn("void donorBasedLayeredOverlapFringe::init() const") + << "Region " << crName << " specified as connected region, but" + << " that region has " + << allRegions[regionID].donorRegions().size() + << " donor regions." + << abort(FatalError); + } + + // Sanity check whether the donor region of connected region is actually + // this region + if (allRegions[regionID].donorRegions()[0] != this->region().index()) + { + FatalErrorIn("void donorBasedLayeredOverlapFringe::init() const") + << "The donor region of region " << crName + << " should be only region " << this->region().name() + << abort(FatalError); + } + } + + // Sanity check: number of (optionally) specified centre points must be + // equal to the number of connected regions + if + ( + !regionCentrePoints_.empty() + && (regionCentrePoints_.size() != connectedRegionIDs_.size()) + ) + { + FatalErrorIn("void donorBasedLayeredOverlapFringe::init() const") + << "You have specified " + << regionCentrePoints_.size() + << " regionCentrePoints, while specifying " + << connectedRegionIDs_.size() + << " connectedRegions." + << nl + << "If you'd like to avoid using automatic centre point detection," + << " make sure to specify centre points for all connected regions." + << abort(FatalError); + } +} + + +void Foam::donorBasedLayeredOverlapFringe::calcAddressing() const +{ + // Make sure that either acceptorsPtr is unnalocated or if it is allocated, + // that it is empty + if (acceptorsPtr_ && !acceptorsPtr_->empty()) + { + FatalErrorIn + ( + "void Foam::donorBasedLayeredOverlapFringe::calcAddressing() const" + ) << "Addressing already calculated" + << abort(FatalError); + } + + if (!isInitialized_) + { + // This is the first call, initialize the data and set flag to true + init(); + isInitialized_ = true; + } + + // Get list of all overset regions + const PtrList& allRegions = + this->region().overset().regions(); + + // Loop through all connected regions and check whether the fringe overlap + // has been found for all of them + bool allFringesReady = true; + forAll (connectedRegionIDs_, crI) + { + // Get ID of this region + const label& regionID = connectedRegionIDs_[crI]; + + // Get the overset region + const oversetRegion& region = allRegions[regionID]; + + // Get the fringe algorithm from the region + const oversetFringe& fringe = region.fringe(); + + // If this is not faceCells fringe, issue a Warning. This fringe + // selection algorithm is intended to work only with faceCells fringe on + // the other side. VV, 9/Apr/2019 + if (!isA(fringe)) + { + WarningIn + ( + "void Foam::donorBasedLayeredOverlapFringe::" + "updateIteration(donorAcceptorList&) const" + ) << "donorBasedLayeredOverlap fringe is designed to work" + << " with faceCells fringe as a connected region fringe." + << nl + << "Connected overset region " << region.name() + << " has " << fringe.type() << " fringe type. " + << nl + << "Proceed with care!" + << endl; + } + + // Update flag collecting whether all connected regions found the + // overlap + allFringesReady &= fringe.foundSuitableOverlap(); + } + + // Sets containing all acceptors and all holes for all connected regions + const polyMesh& mesh = this->mesh(); + labelHashSet allAcceptors(0.02*mesh.nCells()); + labelHashSet allFringeHoles(0.02*mesh.nCells()); + + // Communicate state across processors + reduce(allFringesReady, andOp()); + + if (allFringesReady) + { + if (debug) + { + Info<< "All dependent fringes are ready." + << " Starting donor based layered overlap assembly..." << endl; + } + + // Loop through connected regions + forAll (connectedRegionIDs_, crI) + { + // Get ID of this region + const label& regionID = connectedRegionIDs_[crI]; + + // Get fringe of the connected region + const oversetFringe& fringe = allRegions[regionID].fringe(); + + // The fringe should be finalized, which means we may take a const + // reference to its final donor acceptors + const donorAcceptorList& crDonorAcceptorPairs = + fringe.finalDonorAcceptors(); + + // Need to gather/scatter the donor-acceptor pairs across all + // processors because these pairs only represent acceptors found on + // my processor. It would be possible to optimize this a bit using + // the mapDistribute tool, but I don't think it will represent a big + // overhead, especially since it's done once. + List allDonorAcceptorPairs(Pstream::nProcs()); + + // Fill in my part and communicate + allDonorAcceptorPairs[Pstream::myProcNo()] = crDonorAcceptorPairs; + Pstream::gatherList(allDonorAcceptorPairs); + Pstream::scatterList(allDonorAcceptorPairs); + + // Count approximate number of acceptors to guess the size for the + // hash set containing donors + label nAllAcceptors = 0; + forAll (allDonorAcceptorPairs, procI) + { + nAllAcceptors += allDonorAcceptorPairs[procI].size(); + } + + // Hash set containing donors + labelHashSet donors(6*nAllAcceptors); + + // Initialize centre of the donors of this connected region in order + // to search in a given direction + vector centrePoint(vector::zero); + + // Loop through all processors + forAll (allDonorAcceptorPairs, procI) + { + // Get all donor/acceptor pairs found on this processor + const donorAcceptorList& procDonorAcceptorPairs = + allDonorAcceptorPairs[procI]; + + // Loop through all donor/acceptors + forAll (procDonorAcceptorPairs, daI) + { + // Get this donor/acceptor pair + const donorAcceptor& daPair = procDonorAcceptorPairs[daI]; + + // Check whether all donors have been found + if (!daPair.donorFound()) + { + FatalErrorIn + ( + "donorBasedLayeredOverlapFringe::" + "updateIteration(donorAcceptorList&) const" + ) << "Donor not found for donor/acceptor pair " << daI + << nl + << "Donor/acceptor data: " << daPair + << nl + << "In connected region: " + << allRegions[regionID].name() + << abort(FatalError); + } + + // Mark donors on my processor from this connected region. + // Note that the check has been made in constructor to make + // sure that this region is the only donor region for the + // connected region + if (daPair.donorProcNo() == Pstream::myProcNo()) + { + // Get donor index + const label& dI = daPair.donorCell(); + + // Insert donor into the hash set + if (donors.insert(dI)) + { + // Donor has been inserted (not previously found in + // the hash set), add donor point to centre point + // (the centre point will be calculated later on as + // arithmetic mean) + centrePoint += daPair.donorPoint(); + } + + // Loop through extended donor cells + const donorAcceptor::DynamicLabelList& extDonors = + daPair.extendedDonorCells(); + const donorAcceptor::DynamicPointList& extDonorPoints = + daPair.extendedDonorPoints(); + + forAll (extDonors, i) + { + // Get donor index + const label& edI = extDonors[i]; + + // Inser extended donor into the hash set + if (donors.insert(edI)) + { + // Donor has been inserted (not previously found + // in the hash set), add extended donor point as + // well + centrePoint += extDonorPoints[i]; + } + } // End for all extended donors + } // End if this donor is on my processor + } // End for all (master) donor cells + } // End for all processors + + // Use the centre point as specified by the user if it was specified + // (if the regionCentrePoints_ list is not empty). This avoids + // parallel communication as well. + if (!regionCentrePoints_.empty()) + { + // Use specified centre point, discarding the data we calculated + // above + centrePoint = regionCentrePoints_[crI]; + } + else + { + // User did not specify centre points and the centre point holds + // the sum of all the points. Reduce centre point and divide it + // with global number of unique donors + reduce(centrePoint, sumOp()); + centrePoint /= returnReduce(donors.size(), sumOp