exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Additional overset fringe strategy. Vuko Vukcevic

Browse source 
cuttingPatchFringe (based on cutting patches)
This commit is contained in:
Hrvoje Jasak 2019-07-11 13:42:23 +01:00
commit 1416e3c96c
14 changed files with 1963 additions and 28 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
src/overset/oversetMesh/Make/files
View file

@ -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

6
src/overset/oversetMesh/Make/options
View file

@ -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

10
src/overset/oversetMesh/oversetAdjustPhi/oversetAdjustPhi.C
View file

@ -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
(

19
src/overset/oversetMesh/oversetAdjustPhi/regionWiseOversetAdjustPhi.C
View file

@ -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]];

649
src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<triFaceList> allTriFaces(Pstream::nProcs());
List<pointField> 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<bool>(),
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<bool>(),
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);
}
// ************************************************************************* //

151
src/overset/oversetMesh/oversetFringe/cuttingPatchFringe/cuttingPatchFringe.H Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
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
// ************************************************************************* //

887
src/overset/oversetMesh/oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.C Normal file
View file

@ -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 <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#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<oversetRegion>& 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<oversetRegion>& 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<faceCellsFringe>(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<bool>());
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<donorAcceptorList> 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<vector>());
centrePoint /= returnReduce(donors.size(), sumOp<label>());
}
if (debug)
{
Info<< "Centre point for donors for region "
<< allRegions[regionID].name() << " is : " << centrePoint
<< endl;
}
// We now have a collection of all donors for this connected region
// and the centre point to move to. Let's collect the acceptors
labelHashSet acceptors(donors.size()); // Reasonable size estimate
// Get necessary mesh data (from polyMesh/primitiveMesh)
const vectorField& cc = mesh.cellCentres();
const vectorField& fc = mesh.faceCentres();
const cellList& meshCells = mesh.cells();
const unallocLabelList& owner = mesh.faceOwner();
const unallocLabelList& neighbour = mesh.faceNeighbour();
// Get bounding box of this region for additional check when marking
// acceptors
const boundBox& bb = allRegions[regionID].globalBounds();
// Mark cells that are eligible to be acceptors (not donors)
boolList eligibleCells(mesh.nCells(), true);
forAllConstIter (labelHashSet, donors, iter)
{
eligibleCells[iter.key()] = false;
}
// 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 donors and mark faces that are pointing
// towards the centre point and have an eligible neighbour
forAllConstIter (labelHashSet, donors, iter)
{
// Get the cell index and the cell
const label& cellI = iter.key();
const cell& cFaces = meshCells[cellI];
// Get cell centre of this donor and calculate distance to
// centre point
const vector& donorCentre = cc[cellI];
const scalar donorCentreToRegionCentreDist =
mag(donorCentre - centrePoint);
// Loop through all faces of the cell
forAll (cFaces, i)
{
// Get face index (global)
const label& faceI = cFaces[i];
// Get face centre and calculate distance to centre
// point
const vector& faceCentre = fc[faceI];
const scalar faceCentreToRegionCentreDist =
mag(faceCentre - centrePoint);
if
(
// First, distance based criterium
(
faceCentreToRegionCentreDist
- donorCentreToRegionCentreDist
< distTol_
)
// Second, bounding box based criterium
&& (bb.contains(faceCentre))
)
{
// Face is closer to the centre point than cell: we
// are moving in the right direction. Mark the face
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,
orOp<bool>(),
false
);
// Loop through all 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 eligible, not both
if (eligibleCells[own])
{
// Owner cell is not a donor, insert it
acceptors.insert(own);
// Mark as ineligible
eligibleCells[own] = false;
}
else if (eligibleCells[nei])
{
// Neighbour cell is not a donor, insert it
acceptors.insert(nei);
// Mark as ineligible
eligibleCells[nei] = 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 (eligibleCells[own])
{
// Face cell is not a donor, insert it
acceptors.insert(own);
// Mark as ineligible
eligibleCells[own] = false;
}
}
}
// Special treatment for all non-final iterations
if (i < nLayers_ - 1)
{
// This is not the last iteration, transfer acceptors into
// donors
donors.transfer(acceptors);
// Resize acceptors list
acceptors.resize(donors.size());
}
} // End for specified number of layers
// At this point, we have the final set of acceptors and we marked
// all cells that are ineligible (either donor or acceptor). The
// remaining thing to do is to mark the interior holes
// Create a hash set that will contain fringe holes
labelHashSet fringeHoles(10*acceptors.size());
// Collect holes until there are no holes to collect. Note: for the
// first iteration, we will start with acceptors, otherwise, we will
// start from all fringe holes
label nAddedHoles;
bool firstIteration = true;
do
{
// Reset number of newly added holes
nAddedHoles = 0;
// Face markup for propagation
boolList propagateFace(mesh.nFaces(), false);
// Get collection of cells from which to start the search:
// acceptors for the first iteration and fringe holes otherwise
labelHashSet* curSetPtr = nullptr;
if (firstIteration)
{
// Point current set to acceptors and switch off the flag
curSetPtr = &acceptors;
firstIteration = false;
}
else
{
// Point current set to fringe holes
curSetPtr = &fringeHoles;
}
const labelHashSet& curSet = *curSetPtr;
// Loop through all cells in the set and mark faces that are
// pointing towards the centre point and have eligible neighbour
// (not an acceptor or donor).
forAllConstIter (labelHashSet, curSet, iter)
{
// Get the cell index and the cell
const label& cellI = iter.key();
const cell& cFaces = meshCells[cellI];
// Note: there's no need to check for the distance here
// because there's always at least one "buffer" layer
// towards the outer side that consists of donors, which are
// marked as ineligible at the beginning
// Loop through all faces of the cell and mark all of them
// for propagation
forAll (cFaces, i)
{
propagateFace[cFaces[i]] = true;
}
}
// Sync the face list across processor boundaries
syncTools::syncFaceList
(
mesh,
propagateFace,
orOp<bool>(),
false
);
// Loop through all faces and append interior holes
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 eligible, not both
if (eligibleCells[own])
{
// Owner cell is not a hole (or an acceptor in the
// first iteration), insert it
if (fringeHoles.insert(own))
{
// Count number of added holes in this iteration
++nAddedHoles;
}
// Mark as ineligible
eligibleCells[own] = false;
}
else if (eligibleCells[nei])
{
// Neighbour cell is not a hole (or an acceptor in
// the first iteration), insert it
if (fringeHoles.insert(nei))
{
// Count number of added holes in this iteration
++nAddedHoles;
}
// Mark as ineligible
eligibleCells[nei] = 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 (eligibleCells[own])
{
// Face cell is not a hole (or an acceptor in the
// first iteration), inser it
if (fringeHoles.insert(own))
{
++nAddedHoles;
}
// Mark as ineligible
eligibleCells[own] = false;
}
}
}
// Need to reduce nAddedHoles in order to have synced loops in
// parallel
reduce(nAddedHoles, sumOp<label>());
// We moved one layer "inside" the fringe. Keep going until
// there are no more holes to add
} while (nAddedHoles != 0);
// Finally, we have collected all the fringe holes and acceptors
// from this connected region. Append them to the global sets
allAcceptors += acceptors;
allFringeHoles += fringeHoles;
} // End for all connected regions
// Set acceptors and holes from the data for all regions
acceptorsPtr_ = new labelList(allAcceptors.sortedToc());
fringeHolesPtr_ = new labelList(allFringeHoles.sortedToc());
}
else
{
// Connected fringes are not ready, allocate empty lists for acceptors
// and holes, which will be deleted when asked for again from the
// iterative procedure (see candidateAcceptors() and fringeHoles())
acceptorsPtr_ = new labelList;
fringeHolesPtr_ = new labelList;
}
if (debug)
{
Info<< "In donorBasedLayeredOverlapFringe::calcAddressing() const" << nl
<< "Found " << acceptorsPtr_->size() << " acceptors." << nl
<< "Found " << fringeHolesPtr_->size() << " fringe holes." << endl;
}
}
void Foam::donorBasedLayeredOverlapFringe::clearAddressing() const
{
deleteDemandDrivenData(fringeHolesPtr_);
deleteDemandDrivenData(acceptorsPtr_);
deleteDemandDrivenData(finalDonorAcceptorsPtr_);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct from dictionary
Foam::donorBasedLayeredOverlapFringe::donorBasedLayeredOverlapFringe
(
const fvMesh& mesh,
const oversetRegion& region,
const dictionary& dict
)
:
oversetFringe(mesh, region, dict),
connectedRegionNames_(dict.lookup("connectedRegions")),
connectedRegionIDs_(),
regionCentrePoints_
(
dict.lookupOrDefault<List<point> >
(
"regionCentrePoints",
List<point>(0)
)
),
nLayers_(readLabel(dict.lookup("nLayers"))),
fringeHolesPtr_(nullptr),
acceptorsPtr_(nullptr),
finalDonorAcceptorsPtr_(nullptr),
isInitialized_(false)
{
// Sanity check number of layers: must be greater than 0
if (nLayers_ < 1)
{
FatalIOErrorIn
(
"donorBasedLayeredOverlapFringe::"
"donorBasedLayeredOverlapFringe\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."
<< abort(FatalError);
}
// Recommendation: use at least two layers in order to avoid going defining
// acceptors on the wrong side and filling in the whole region with holes
if (nLayers_ == 1)
{
WarningIn
(
"donorBasedLayeredOverlapFringe::"
"donorBasedLayeredOverlapFringe\n"
"(\n"
" const fvMesh& mesh,\n"
" const oversetRegion& region,\n"
" const dictionary& dict\n"
")"
) << "It is recommended to use at least 2 layers (nLayers = 2) in"
<< " order to avoid defining acceptor cells on the wrong side"
<< " of the region."
<< nl
<< "Be sure to check the fringe layer for this region."
<< endl;
}
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::donorBasedLayeredOverlapFringe::~donorBasedLayeredOverlapFringe()
{
clearAddressing();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
bool Foam::donorBasedLayeredOverlapFringe::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
// donorBasedLayeredOverlapFringe
if (finalDonorAcceptorsPtr_)
{
FatalErrorIn
(
"donorBasedLayeredOverlapFringe::"
"updateIteration(donorAcceptorList&) const"
) << "finalDonorAcceptorPtr_ already allocated. Something went "
<< "wrong with the iteration procedure (flag was not updated)."
<< nl
<< "This should not happen for donorBasedLayeredOverlapFringe."
<< abort(FatalError);
}
if
(
fringeHolesPtr_ && acceptorsPtr_
&& returnReduce(!acceptorsPtr_->empty(), orOp<bool>())
)
{
// Note: we first check whether fringeHoles and acceptors pointers are
// allocated. If they are, there must be at least one processor with
// more than 0 acceptors in order for this fringe to be valid.
// Allocate the list by reusing the argument list
finalDonorAcceptorsPtr_ = new donorAcceptorList
(
donorAcceptorRegionData,
true
);
// Set the flag to true (for all processors due to reduction in the if
// statement)
updateSuitableOverlapFlag(true);
}
else
{
// Delete fringeHolesPtr and acceptorsPtr to trigger calculation of
// addressing, this time with other fringes up-to-date
deleteDemandDrivenData(fringeHolesPtr_);
deleteDemandDrivenData(acceptorsPtr_);
}
return foundSuitableOverlap();
}
const Foam::labelList& Foam::donorBasedLayeredOverlapFringe::fringeHoles() const
{
// Note: updateIteration deletes fringeHolesPtr if the suitable overlap is
// not found, thus preparing it for the next iteration when the other
// fringes should be ready
if (!fringeHolesPtr_)
{
calcAddressing();
}
return *fringeHolesPtr_;
}
const Foam::labelList&
Foam::donorBasedLayeredOverlapFringe::candidateAcceptors() const
{
// Note: updateIteration deletes acceptorsPtr if the suitable overlap is
// not found, thus preparing it for the next iteration when the other
// fringes should be ready
if (!acceptorsPtr_)
{
calcAddressing();
}
return *acceptorsPtr_;
}
Foam::donorAcceptorList&
Foam::donorBasedLayeredOverlapFringe::finalDonorAcceptors() const
{
if (!finalDonorAcceptorsPtr_)
{
FatalErrorIn("donorBasedLayeredOverlapFringe::finalDonorAcceptors()")
<< "finalDonorAcceptorPtr_ not allocated. Make sure you have"
<< " called donorBasedLayeredOverlapFringe::updateIteration() before"
<< " asking for final set of donor/acceptor pairs."
<< abort(FatalError);
}
if (!foundSuitableOverlap())
{
FatalErrorIn("donorBasedLayeredOverlapFringe::finalDonorAcceptors()")
<< "Attemted to access finalDonorAcceptors but suitable overlap "
<< "has not been found. This is not allowed. "
<< abort(FatalError);
}
return *finalDonorAcceptorsPtr_;
}
void Foam::donorBasedLayeredOverlapFringe::update() const
{
Info<< "donorBasedLayeredOverlapFringe::update() const" << endl;
// Clear out
clearAddressing();
// Set flag to false and clear final donor/acceptors only
deleteDemandDrivenData(finalDonorAcceptorsPtr_);
updateSuitableOverlapFlag(false);
}
// ************************************************************************* //

190
src/overset/oversetMesh/oversetFringe/donorBasedLayeredOverlapFringe/donorBasedLayeredOverlapFringe.H Normal file
View file

@ -0,0 +1,190 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / 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 <http://www.gnu.org/licenses/>.
Class
donorBasedLayeredOverlapFringe
Description
Fringe algorithm based on donors from different regions. The algorithm waits
until the final donor/acceptor assembly has been performed for all regions
whose donor region is this region. Then, all donors are collected and the
acceptors are cells neighbouring the donors nLayers towards the interior.
Interior is defined either by user-specified points for each region or as a
centre of volume of the donor cells in particular region.
This fringe algorithm is intended to be used along with the faceCells fringe
on the other side, where the cell sizes are not significantly different from
each other (e.g. 10:1 cell ratio, where the finer cells are found on the
background mesh will probably be problematic to correctly set-up).
Note: based on distance tolerance (see distTol_ member), specified number of
layers and the bounding box of the connected donor region, it is possible
that acceptors end also on the wrong side (from the region centre as opposed
to towards the region centre). I do all I can to prevent this by looking at
the bounding box and hinting that at least 2 layers should be used. Be sure
to check the overset assembly for this region using calcOverset utility.
Author
Vuko Vukcevic, Wikki Ltd. All rights reserved.
SourceFiles
donorBasedLayeredOverlapFringe.C
\*---------------------------------------------------------------------------*/
#ifndef donorBasedLayeredOverlapFringe_H
#define donorBasedLayeredOverlapFringe_H
#include "oversetFringe.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class donorBasedLayeredOverlapFringe Declaration
\*---------------------------------------------------------------------------*/
class donorBasedLayeredOverlapFringe
:
public oversetFringe
{
// Private data
//- Names of connected regions. Looked up on construction
wordList connectedRegionNames_;
//- Regions IDs from which the donors will be collected as a starting
// point. Note: initialized in init private member function because we
// cannot initialize it in constructor. This is because certain overset
// regions (and their fringes) may not be initialized at this point.
mutable labelList connectedRegionIDs_;
//- Optional list of points representing a rough estimate of the centre
// for each underlying connected region. If these are not provided, the
// centre is calculated as the centre of all donors for a given
// connected region
List<point> regionCentrePoints_;
//- 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_;
//- Initialization helper
mutable bool isInitialized_;
// Private static data
//- Distance tolerance to determine propagation direction. Note:
// absolute value, default = 0. Might be useful is some strange cases
static const debug::tolerancesSwitch distTol_;
// Private Member Functions
//- Initialization
void init() const;
//- Calculate hole and acceptor addressing
void calcAddressing() const;
//- Clear addressing
void clearAddressing() const;
public:
//- Runtime type information
TypeName("donorBasedLayeredOverlap");
// Constructors
//- Construct from dictionary
donorBasedLayeredOverlapFringe
(
const fvMesh& mesh,
const oversetRegion& region,
const dictionary& dict
);
//- Disallow default bitwise copy construct
donorBasedLayeredOverlapFringe
(
const donorBasedLayeredOverlapFringe&
) = delete;
//- Disallow default bitwise assignment
void operator=(const donorBasedLayeredOverlapFringe&) = delete;
// Destructor
virtual ~donorBasedLayeredOverlapFringe();
// 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
// ************************************************************************* //

2
src/overset/oversetMesh/oversetFringe/faceCellsFringe/faceCellsFringe.C
View file

@ -72,7 +72,7 @@ void Foam::faceCellsFringe::calcAddressing() const
(
"void faceCellsFringe::calcAddressing() const"
) << "Fringe patch " << patchNames_[nameI]
<< " cannot be found"
<< " cannot be found."
<< abort(FatalError);
}

6
src/overset/oversetMesh/oversetFringe/faceCellsFringe/faceCellsFringe.H
View file

@ -115,6 +115,12 @@ public:
// Member Functions
//- Access to patch names
const wordList& patchNames() const
{
return patchNames_;
}
//- Update iteration. Note: invalidates parameter
virtual bool updateIteration
(

7
src/overset/oversetMesh/oversetMesh/oversetMesh.H
View file

@ -121,8 +121,9 @@ private:
//- Return overset type indicator field
mutable volScalarField* oversetTypesPtr_;
//- Region ID: region index for each cell
mutable labelList* regionIDPtr_;
//- Region ID: region index for each cell as a volScalarField for
// visualization. VV, 15/Apr/2019
mutable volScalarField* regionIDPtr_;
// Overset discretisation support
@ -286,7 +287,7 @@ public:
const volScalarField& oversetTypes() const;
//- Return region indicator
const labelList& regionID() const;
const volScalarField& regionID() const;
// Overset discretisation support

45
src/overset/oversetMesh/oversetMesh/oversetMeshAddressing.C
View file

@ -27,6 +27,7 @@ License
#include "surfaceFields.H"
#include "volFields.H"
#include "polyPatchID.H"
#include "oversetFvPatchFields.H"
#include "demandDrivenData.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -190,24 +191,52 @@ void Foam::oversetMesh::calcDomainMarkup() const
}
}
// Region ID
regionIDPtr_ = new labelList(mesh().nCells(), -1);
labelList& rID = *regionIDPtr_;
// Mark regions
// Region ID, initialized with -1 for sanity check later on
regionIDPtr_ = new volScalarField
(
IOobject
(
"regionIndex",
mesh().time().timeName(),
mesh(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh(),
dimensionedScalar("minusOne", dimless, -1.0),
"zeroGradient"
);
volScalarField& regionID = *regionIDPtr_;
scalarField& regionIDIn = regionID.internalField();
// Mark regions (internal field)
forAll (regions_, regionI)
{
const labelList& curCells = regions_[regionI].zone();
forAll (curCells, curCellI)
{
rID[curCells[curCellI]] = regionI;
regionIDIn[curCells[curCellI]] = regionI;
}
}
// Update boundary values, making sure that we skip the overset patch
volScalarField::GeometricBoundaryField& regionIDb =
regionID.boundaryField();
forAll (regionIDb, patchI)
{
// Get the patch field
fvPatchScalarField& ripf = regionIDb[patchI];
if (!isA<oversetFvPatchScalarField>(ripf))
{
ripf = ripf.patchInternalField();
}
}
// Check regions
if (min(rID) < 0)
if (min(regionID).value() < 0)
{
FatalErrorIn("void oversetMesh::calcDomainMarkup() const")
<< "Found cells without region ID. Please check overset setup"
@ -1166,7 +1195,7 @@ const Foam::volScalarField& Foam::oversetMesh::oversetTypes() const
}
const Foam::labelList& Foam::oversetMesh::regionID() const
const Foam::volScalarField& Foam::oversetMesh::regionID() const
{
if (!regionIDPtr_)
{

14
src/overset/oversetMesh/oversetRegion/oversetRegion.C
View file

@ -1730,6 +1730,20 @@ Foam::oversetRegion::~oversetRegion()
// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
const Foam::oversetFringe& Foam::oversetRegion::fringe() const
{
if (fringePtr_.empty())
{
FatalErrorIn("const oversetFringe& oversetRegion::fringe() const")
<< "Fringe pointer not allocated. It should have been initialized"
<< " properly at construction. Something went wrong..."
<< abort(FatalError);
}
return fringePtr_();
}
const Foam::labelList& Foam::oversetRegion::donorRegions() const
{
if (!donorRegionsPtr_)

3
src/overset/oversetMesh/oversetRegion/oversetRegion.H
View file

@ -327,6 +327,9 @@ public:
return zone();
}
//- Return overset fringe algorithm for this region
const oversetFringe& fringe() const;
//- Return list of donor region indices
const labelList& donorRegions() const;