zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'feature-finiteArea-fixes' into 'develop'

Browse Source 
parallel construct finiteArea with arbitrary connections

See merge request Development/openfoam!490
This commit is contained in:
Andrew Heather
2021-11-05 18:08:12 +00:00
parent 38bf30167d 7fc943c178
commit 6102b76377
31 changed files with 3217 additions and 1147 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
applications/utilities/finiteArea/checkFaMesh/Make/options
View File

@ -1,9 +1,11 @@
EXE_INC = \
-I$(LIB_SRC)/finiteArea/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude
EXE_LIBS = \
-lfiniteArea \
-lfiniteVolume \
-lfileFormats \
-lsurfMesh \
-lmeshTools

54
applications/utilities/finiteArea/checkFaMesh/checkFaMesh.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -29,14 +30,22 @@ Application
Description
Check a finiteArea mesh
Author
Original Authors
Zeljko Tukovic, FAMENA
Hrvoje Jasak, Wikki Ltd.
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "faCFD.H"
#include "Time.H"
#include "argList.H"
#include "faMesh.H"
#include "polyMesh.H"
#include "areaFaMesh.H"
#include "edgeFaMesh.H"
#include "areaFields.H"
#include "edgeFields.H"
#include "processorFaPatch.H"
#include "foamVtkUIndPatchWriter.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -51,37 +60,32 @@ int main(int argc, char *argv[])
"Check a finiteArea mesh"
);
#include "addRegionOption.H"
argList::addBoolOption
(
"write-vtk",
"Write mesh as a vtp (vtk) file for display or debugging"
);
#include "addRegionOption.H"
#include "setRootCase.H"
#include "createTime.H"
#include "createNamedMesh.H"
#include "createFaMesh.H"
#include "createNamedPolyMesh.H"
// Create
faMesh aMesh(mesh);
Info<< "Time = " << runTime.timeName() << nl << endl;
// General mesh statistics
Info<< "Number of points: " << aMesh.nPoints() << nl
<< "Number of internal edges: " << aMesh.nInternalEdges() << nl
<< "Number of edges: " << aMesh.nEdges() << nl
<< "Number of faces: " << aMesh.nFaces() << nl
<< endl;
#include "printMeshSummary.H"
// Check geometry
Info<< "Face area: min = " << min(aMesh.S().field())
<< " max = " << max(aMesh.S().field()) << nl
<< "Internal edge length: min = "
<< min(aMesh.magLe().internalField()) << nl
<< " max = " << max(aMesh.magLe().internalField()) << nl
<< "Edge length: min = "
<< min(aMesh.magLe()).value() << nl
<< " max = " << max(aMesh.magLe()).value() << nl
<< "Face area normals: min = " << min(aMesh.faceAreaNormals().field())
<< " max = " << max(aMesh.faceAreaNormals().field()) << nl
<< endl;
if (args.found("write-vtk"))
{
#include "faMeshWriteVTK.H"
}
Info<< "\nEnd\n" << endl;
Info << "\nEnd" << endl;
return 0;
}

47
applications/utilities/finiteArea/checkFaMesh/faMeshWriteVTK.H Normal file
View File

@ -0,0 +1,47 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
Description
VTK output of faMesh with some geometric or debug fields
\*---------------------------------------------------------------------------*/
{
// finiteArea
vtk::uindirectPatchWriter writer
(
aMesh.patch(),
// vtk::formatType::INLINE_ASCII,
fileName
(
aMesh.mesh().time().globalPath() / "finiteArea"
)
);
writer.writeGeometry();
// CellData
writer.beginCellData(3);
writer.writeProcIDs();
writer.write("area", aMesh.S().field());
writer.write("normal", aMesh.faceAreaNormals());
// PointData
writer.beginPointData(1);
writer.write("normal", aMesh.pointAreaNormals());
Info<< nl
<< "Wrote faMesh in vtk format: " << writer.output().name() << nl;
}
// ************************************************************************* //

134
applications/utilities/finiteArea/checkFaMesh/printMeshSummary.H Normal file
View File

@ -0,0 +1,134 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
Description
Summary of faMesh information
\*---------------------------------------------------------------------------*/
{
const faBoundaryMesh& patches = aMesh.boundary();
const label nNonProcessor = patches.nNonProcessor();
const label nPatches = patches.size();
Info<< "----------------" << nl
<< "Mesh Information" << nl
<< "----------------" << nl
<< " " << "boundingBox: " << boundBox(aMesh.points()) << nl;
Info<< " Number of points: "
<< returnReduce(aMesh.nPoints(), sumOp<label>()) << nl
<< " Number of faces: "
<< returnReduce(aMesh.nFaces(), sumOp<label>()) << nl;
Info<< " Number of edges: "
<< returnReduce(aMesh.nEdges(), sumOp<label>()) << nl
<< " Number of internal edges: "
<< returnReduce(aMesh.nInternalEdges(), sumOp<label>()) << nl;
label nProcEdges = 0;
if (Pstream::parRun())
{
for (const faPatch& fap : patches)
{
const auto* cpp = isA<processorFaPatch>(fap);
if (cpp)
{
nProcEdges += fap.nEdges();
}
}
}
const label nBoundEdges = aMesh.nBoundaryEdges() - nProcEdges;
Info<< " Number of boundary edges: "
<< returnReduce(nBoundEdges - nProcEdges, sumOp<label>()) << nl;
if (Pstream::parRun())
{
Info<< " Number of processor edges: "
<< returnReduce(nProcEdges, sumOp<label>()) << nl;
}
Info<< "----------------" << nl
<< "Patches" << nl
<< "----------------" << nl;
for (label patchi = 0; patchi < nNonProcessor; ++patchi)
{
const faPatch& p = patches[patchi];
Info<< " " << "patch " << p.index()
<< " (size: " << returnReduce(p.size(), sumOp<label>())
<< ") name: " << p.name()
<< nl;
}
// Geometry information
Info<< nl;
{
scalarMinMax limit(gMinMax(aMesh.S().field()));
Info<< "Face area:" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
}
{
scalarMinMax limit(minMax(aMesh.magLe().primitiveField()));
// Include processor boundaries into 'internal' edges
if (Pstream::parRun())
{
for (label patchi = nNonProcessor; patchi < nPatches; ++patchi)
{
limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
}
reduce(limit, minMaxOp<scalar>());
}
Info<< "Edge length (internal):" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
// Include (non-processor) boundaries
for (label patchi = 0; patchi < nNonProcessor; ++patchi)
{
limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
}
if (Pstream::parRun())
{
reduce(limit, minMaxOp<scalar>());
}
Info<< "Edge length:" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
}
// Not particularly meaningful
#if 0
{
MinMax<vector> limit(gMinMax(aMesh.faceAreaNormals().field()));
Info<< "Face area normals:" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
}
#endif
}
// ************************************************************************* //

2
applications/utilities/finiteArea/makeFaMesh/Make/options
View File

@ -1,6 +1,7 @@
EXE_INC = \
-I$(LIB_SRC)/finiteArea/lnInclude \
-I$(LIB_SRC)/fileFormats/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/parallel/decompose/faDecompose/lnInclude \
-I$(LIB_SRC)/parallel/reconstruct/faReconstruct/lnInclude
@ -8,6 +9,7 @@ EXE_INC = \
EXE_LIBS = \
-lfiniteArea \
-lfileFormats \
-lsurfMesh \
-lmeshTools \
-lfaDecompose \
-lfaReconstruct

49
applications/utilities/finiteArea/makeFaMesh/decomposeFaFields.H
View File

@ -17,9 +17,11 @@ Description
if (Pstream::parRun())
{
faMeshReconstructor reconstructor(areaMesh);
faMeshReconstructor reconstructor(aMesh);
reconstructor.writeAddressing();
Info<< "Wrote proc-addressing" << nl << endl;
// Handle area fields
// ------------------
@ -72,22 +74,53 @@ if (Pstream::parRun())
if (nAreaFields)
{
Info<< "Decomposing " << nAreaFields << " area fields" << nl << endl;
Info<< "Decomposing " << nAreaFields << " area fields" << nl;
faFieldDecomposer fieldDecomposer
(
fullMesh,
areaMesh,
aMesh,
reconstructor.edgeProcAddressing(),
reconstructor.faceProcAddressing(),
reconstructor.boundaryProcAddressing()
);
fieldDecomposer.decomposeFields(areaScalarFields);
fieldDecomposer.decomposeFields(areaVectorFields);
fieldDecomposer.decomposeFields(areaSphTensorFields);
fieldDecomposer.decomposeFields(areaSymmTensorFields);
fieldDecomposer.decomposeFields(areaTensorFields);
if (areaScalarFields.size())
{
Info<< " scalars: "
<< flatOutput(PtrListOps::names(areaScalarFields)) << nl;
fieldDecomposer.decomposeFields(areaScalarFields);
}
if (areaVectorFields.size())
{
Info<< " vectors: "
<< flatOutput(PtrListOps::names(areaVectorFields)) << nl;
fieldDecomposer.decomposeFields(areaVectorFields);
}
if (areaSphTensorFields.size())
{
Info<< " sphTensors: "
<< flatOutput(PtrListOps::names(areaSphTensorFields)) << nl;
fieldDecomposer.decomposeFields(areaSphTensorFields);
}
if (areaSymmTensorFields.size())
{
Info<< " symmTensors: "
<< flatOutput(PtrListOps::names(areaSymmTensorFields)) << nl;
fieldDecomposer.decomposeFields(areaSymmTensorFields);
}
if (areaTensorFields.size())
{
Info<< " tensors: "
<< flatOutput(PtrListOps::names(areaTensorFields)) << nl;
fieldDecomposer.decomposeFields(areaTensorFields);
}
Info<< endl;
}
}

11
applications/utilities/finiteArea/makeFaMesh/faMeshWriteEdgesOBJ.H
View File

@ -16,15 +16,18 @@ Description
\*---------------------------------------------------------------------------*/
{
Info<< "Writing edges in obj format" << endl;
Info<< nl
<< "Writing edges in obj format" << endl;
word outputName("faMesh-edges.obj");
word outputName("finiteArea-edges.obj");
if (Pstream::parRun())
{
outputName = word
(
"faMesh-edges-" + Foam::name(Pstream::myProcNo()) + ".obj"
"finiteArea-edges-proc"
+ Foam::name(Pstream::myProcNo())
+ ".obj"
);
}
@ -36,7 +39,7 @@ Description
false
);
os.write(areaMesh.patch().edges(), areaMesh.patch().localPoints());
os.write(aMesh.patch().edges(), aMesh.patch().localPoints());
}

47
applications/utilities/finiteArea/makeFaMesh/faMeshWriteVTK.H Normal file
View File

@ -0,0 +1,47 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
Description
VTK output of faMesh with some geometric or debug fields
\*---------------------------------------------------------------------------*/
{
// finiteArea
vtk::uindirectPatchWriter writer
(
aMesh.patch(),
// vtk::formatType::INLINE_ASCII,
fileName
(
aMesh.mesh().time().globalPath() / "finiteArea"
)
);
writer.writeGeometry();
// CellData
writer.beginCellData(3);
writer.writeProcIDs();
writer.write("area", aMesh.S().field());
writer.write("normal", aMesh.faceAreaNormals());
// PointData
writer.beginPointData(1);
writer.write("normal", aMesh.pointAreaNormals());
Info<< nl
<< "Wrote faMesh in vtk format: " << writer.output().name() << nl;
}
// ************************************************************************* //

55
applications/utilities/finiteArea/makeFaMesh/makeFaMesh.C
View File

@ -32,7 +32,7 @@ Description
When called in parallel, it will also try to act like decomposePar,
create procAddressing and decompose serial finite-area fields.
Author
Original Authors
Zeljko Tukovic, FAMENA
Hrvoje Jasak, Wikki Ltd.
@ -48,6 +48,8 @@ Author
#include "areaFields.H"
#include "faFieldDecomposer.H"
#include "faMeshReconstructor.H"
#include "PtrListOps.H"
#include "foamVtkUIndPatchWriter.H"
#include "OBJstream.H"
using namespace Foam;
@ -69,11 +71,21 @@ int main(int argc, char *argv[])
);
argList::addOption("dict", "file", "Alternative faMeshDefinition");
argList::addBoolOption
(
"dry-run",
"Create but do not write"
);
argList::addBoolOption
(
"write-vtk",
"Write mesh as a vtp (vtk) file for display or debugging"
);
argList::addBoolOption
(
"write-edges-obj",
"Write mesh edges as obj files and exit",
false // could make an advanced option
"Write mesh edges as obj files (one per processor)",
true // advanced option (debugging only)
);
#include "addRegionOption.H"
@ -81,6 +93,8 @@ int main(int argc, char *argv[])
#include "createTime.H"
#include "createNamedPolyMesh.H"
const bool dryrun = args.found("dry-run");
// Reading faMeshDefinition dictionary
#include "findMeshDefinitionDict.H"
@ -91,33 +105,40 @@ int main(int argc, char *argv[])
meshDefDict.add("emptyPatch", patchName, true);
}
// Create
faMesh areaMesh(mesh, meshDefDict);
bool quickExit = false;
// Create
faMesh aMesh(mesh, meshDefDict);
// Mesh information
#include "printMeshSummary.H"
if (args.found("write-edges-obj"))
{
quickExit = true;
#include "faMeshWriteEdgesOBJ.H"
}
if (quickExit)
if (args.found("write-vtk"))
{
Info<< "\nEnd\n" << endl;
return 0;
#include "faMeshWriteVTK.H"
}
// Set the precision of the points data to 10
IOstream::defaultPrecision(10);
if (dryrun)
{
Info<< "\ndry-run: not writing mesh or decomposing fields\n" << nl;
}
else
{
// Set the precision of the points data to 10
IOstream::defaultPrecision(10);
Info<< nl << "Write finite area mesh." << nl;
areaMesh.write();
Info<< endl;
Info<< nl << "Write finite area mesh." << nl;
aMesh.write();
#include "decomposeFaFields.H"
Info<< endl;
#include "decomposeFaFields.H"
}
Info << "\nEnd\n" << endl;
Info<< "\nEnd\n" << endl;
return 0;
}

99
applications/utilities/finiteArea/makeFaMesh/printMeshSummary.H Normal file
View File

@ -0,0 +1,99 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
Description
Summary of faMesh information
\*---------------------------------------------------------------------------*/
{
const faBoundaryMesh& patches = aMesh.boundary();
const label nNonProcessor = patches.nNonProcessor();
const label nPatches = patches.size();
Info<< "----------------" << nl
<< "Mesh Information" << nl
<< "----------------" << nl
<< " " << "boundingBox: " << boundBox(aMesh.points()) << nl
<< " " << "nFaces: " << returnReduce(aMesh.nFaces(), sumOp<label>())
<< nl;
Info<< "----------------" << nl
<< "Patches" << nl
<< "----------------" << nl;
for (label patchi = 0; patchi < nNonProcessor; ++patchi)
{
const faPatch& p = patches[patchi];
Info<< " " << "patch " << p.index()
<< " (size: " << returnReduce(p.size(), sumOp<label>())
<< ") name: " << p.name()
<< nl;
}
// Geometry information
Info<< nl;
{
scalarMinMax limit(gMinMax(aMesh.S().field()));
Info<< "Face area:" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
}
{
scalarMinMax limit(minMax(aMesh.magLe().primitiveField()));
// Include processor boundaries into 'internal' edges
if (Pstream::parRun())
{
for (label patchi = nNonProcessor; patchi < nPatches; ++patchi)
{
limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
}
reduce(limit, minMaxOp<scalar>());
}
Info<< "Edge length (internal):" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
// Include (non-processor) boundaries
for (label patchi = 0; patchi < nNonProcessor; ++patchi)
{
limit.add(minMax(aMesh.magLe().boundaryField()[patchi]));
}
if (Pstream::parRun())
{
reduce(limit, minMaxOp<scalar>());
}
Info<< "Edge length:" << nl
<< " min = " << limit.min()
<< " max = " << limit.max() << nl;
}
// Not particularly meaningful
#if 0
{
MinMax<vector> limit(gMinMax(aMesh.faceAreaNormals().field()));
Info<< "Face area normals:" << nl
<< " min = " << limit.min() << " max = " << limit.max() << nl;
}
#endif
}
// ************************************************************************* //

2
src/finiteArea/Make/files
View File

@ -2,7 +2,9 @@ faMesh/faGlobalMeshData/faGlobalMeshData.C
faMesh/faMesh.C
faMesh/faMeshDemandDrivenData.C
faMesh/faMeshPatches.C
faMesh/faMeshTopology.C
faMesh/faMeshUpdate.C
faMesh/faMeshBoundaryHalo.C
faMesh/faBoundaryMesh/faBoundaryMesh.C
faPatches = faMesh/faPatches

9
src/finiteArea/faMesh/faBoundaryMesh/faBoundaryMesh.C
View File

@ -75,8 +75,15 @@ Foam::faBoundaryMesh::faBoundaryMesh
regIOobject(io),
mesh_(mesh)
{
if (readOpt() == IOobject::MUST_READ)
if
(
readOpt() == IOobject::MUST_READ
|| readOpt() == IOobject::MUST_READ_IF_MODIFIED
)
{
// Warn for MUST_READ_IF_MODIFIED
warnNoRereading<faBoundaryMesh>();
faPatchList& patches = *this;
// Read faPatch list

97
src/finiteArea/faMesh/faMesh.C
View File

@ -27,6 +27,7 @@ License
\*---------------------------------------------------------------------------*/
#include "faMesh.H"
#include "faMeshBoundaryHalo.H"
#include "faGlobalMeshData.H"
#include "Time.H"
#include "polyMesh.H"
@ -51,6 +52,8 @@ const Foam::word Foam::faMesh::prefix("finite-area");
Foam::word Foam::faMesh::meshSubDir = "faMesh";
int Foam::faMesh::origPointAreaMethod_ = 0; // Tuning
const int Foam::faMesh::quadricsFit_ = 0; // Tuning
@ -109,17 +112,50 @@ static labelList selectPatchFaces
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::faMesh::checkBoundaryEdgeLabelRange
(
const labelUList& edgeLabels
) const
{
label nErrors = 0;
for (const label edgei : edgeLabels)
{
if (edgei < nInternalEdges_ || edgei >= nEdges_)
{
if (!nErrors++)
{
FatalErrorInFunction
<< "Boundary edge label out of range "
<< nInternalEdges_ << ".." << (nEdges_-1) << nl
<< " ";
}
FatalError<< ' ' << edgei;
}
}
if (nErrors)
{
FatalError << nl << exit(FatalError);
}
}
void Foam::faMesh::initPatch() const
{
if (patchPtr_)
{
delete patchPtr_;
}
patchPtr_ = new uindirectPrimitivePatch
patchPtr_.reset
(
UIndirectList<face>(mesh().faces(), faceLabels_),
mesh().points()
new uindirectPrimitivePatch
(
UIndirectList<face>(mesh().faces(), faceLabels_),
mesh().points()
)
);
bndConnectPtr_.reset(nullptr);
haloMapPtr_.reset(nullptr);
haloFaceCentresPtr_.reset(nullptr);
haloFaceNormalsPtr_.reset(nullptr);
}
@ -168,12 +204,24 @@ void Foam::faMesh::setPrimitiveMeshData()
}
void Foam::faMesh::clearHalo() const
{
DebugInFunction << "Clearing halo information" << endl;
haloMapPtr_.reset(nullptr);
haloFaceCentresPtr_.reset(nullptr);
haloFaceNormalsPtr_.reset(nullptr);
}
void Foam::faMesh::clearGeomNotAreas() const
{
DebugInFunction << "Clearing geometry" << endl;
clearHalo();
patchPtr_.reset(nullptr);
bndConnectPtr_.reset(nullptr);
deleteDemandDrivenData(SPtr_);
deleteDemandDrivenData(patchPtr_);
deleteDemandDrivenData(patchStartsPtr_);
deleteDemandDrivenData(LePtr_);
deleteDemandDrivenData(magLePtr_);
@ -181,7 +229,7 @@ void Foam::faMesh::clearGeomNotAreas() const
deleteDemandDrivenData(edgeCentresPtr_);
deleteDemandDrivenData(faceAreaNormalsPtr_);
deleteDemandDrivenData(edgeAreaNormalsPtr_);
deleteDemandDrivenData(pointAreaNormalsPtr_);
pointAreaNormalsPtr_.reset(nullptr);
deleteDemandDrivenData(faceCurvaturesPtr_);
deleteDemandDrivenData(edgeTransformTensorsPtr_);
}
@ -256,6 +304,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
),
comm_(Pstream::worldComm),
patchPtr_(nullptr),
bndConnectPtr_(nullptr),
lduPtr_(nullptr),
curTimeIndex_(time().timeIndex()),
SPtr_(nullptr),
@ -272,7 +321,11 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
faceCurvaturesPtr_(nullptr),
edgeTransformTensorsPtr_(nullptr),
correctPatchPointNormalsPtr_(nullptr),
globalMeshDataPtr_(nullptr)
globalMeshDataPtr_(nullptr),
haloMapPtr_(nullptr),
haloFaceCentresPtr_(nullptr),
haloFaceNormalsPtr_(nullptr)
{
DebugInFunction << "Creating from IOobject" << endl;
@ -290,7 +343,7 @@ Foam::faMesh::faMesh(const polyMesh& pMesh)
// Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry();
if (isFile(pMesh.time().timePath()/mesh().dbDir()/"S0"))
if (fileHandler().isFile(pMesh.time().timePath()/"S0"))
{
S0Ptr_ = new DimensionedField<scalar, areaMesh>
(
@ -349,6 +402,7 @@ Foam::faMesh::faMesh
),
comm_(Pstream::worldComm),
patchPtr_(nullptr),
bndConnectPtr_(nullptr),
lduPtr_(nullptr),
curTimeIndex_(time().timeIndex()),
SPtr_(nullptr),
@ -365,7 +419,11 @@ Foam::faMesh::faMesh
faceCurvaturesPtr_(nullptr),
edgeTransformTensorsPtr_(nullptr),
correctPatchPointNormalsPtr_(nullptr),
globalMeshDataPtr_(nullptr)
globalMeshDataPtr_(nullptr),
haloMapPtr_(nullptr),
haloFaceCentresPtr_(nullptr),
haloFaceNormalsPtr_(nullptr)
{}
@ -450,7 +508,7 @@ Foam::faMesh::faMesh
// Calculate the geometry for the patches (transformation tensors etc.)
boundary_.calcGeometry();
if (isFile(mesh().time().timePath()/"S0"))
if (fileHandler().isFile(pMesh.time().timePath()/"S0"))
{
S0Ptr_ = new DimensionedField<scalar, areaMesh>
(
@ -647,11 +705,20 @@ const Foam::vectorField& Foam::faMesh::pointAreaNormals() const
{
if (!pointAreaNormalsPtr_)
{
calcPointAreaNormals();
pointAreaNormalsPtr_.reset(new vectorField(nPoints()));
if (origPointAreaMethod_)
{
calcPointAreaNormals_orig(*pointAreaNormalsPtr_);
}
else
{
calcPointAreaNormals(*pointAreaNormalsPtr_);
}
if (quadricsFit_ > 0)
{
calcPointAreaNormalsByQuadricsFit();
calcPointAreaNormalsByQuadricsFit(*pointAreaNormalsPtr_);
}
}

221
src/finiteArea/faMesh/faMesh.H
View File

@ -33,6 +33,9 @@ Description
SourceFiles
faMesh.C
faMeshDemandDrivenData.C
faMeshPatches.C
faMeshTopology.C
faMeshUpdate.C
Author
Zeljko Tukovic, FMENA
@ -68,10 +71,10 @@ namespace Foam
{
// Forward Declarations
class faMeshBoundaryHalo;
class faMeshLduAddressing;
class faMeshMapper;
class faPatchData;
template<class T> class LabelledItem;
/*---------------------------------------------------------------------------*\
Class faMesh Declaration
@ -86,6 +89,107 @@ class faMesh
public faSolution,
public data
{
// Private (internal) classes/structures
//- A (proc, patchi, patchEdgei) tuple used internally for managing
//- patch/patch bookkeeping during construction.
// Finite-area patches are stored with negated indices, which makes
// them readily identifiable and always sort before normal patches.
// Note
struct patchTuple
:
public FixedList<label, 4>
{
// Constructors
// Inherit constructors
using FixedList<label, 4>::FixedList;
//- Default construct as 'invalid'
patchTuple()
{
clear();
}
// Static Member Functions
// Globally consistent ordering
//
// 1. sort left/right as lower/higher processor connection
// 2. sort by proc/patch/patch index
static void sort(UList<Pair<patchTuple>>& list)
{
for (auto& tuples : list)
{
tuples.sort();
}
Foam::stableSort(list);
}
// Member Functions
//- Reset to 'invalid'
void clear()
{
procNo(-1);
patchi(labelMax);
patchEdgei(-1);
meshFacei(-1);
}
//- Valid if proc and edge are non-negative
bool valid() const noexcept
{
return (procNo() >= 0 && patchEdgei() >= 0);
}
// Processor is the first sort index
label procNo() const { return (*this)[0]; }
void procNo(label val) { (*this)[0] = val; }
// PatchId (-ve for finiteArea patches) is the second sort index
label patchi() const { return (*this)[1]; }
void patchi(label val) { (*this)[1] = val; }
// The patch edge index (on the finiteArea patch)
// is the third sort index
label patchEdgei() const { return (*this)[2]; }
void patchEdgei(label val) { (*this)[2] = val; }
// The processor-local mesh face is the fourth sort index
label meshFacei() const { return (*this)[3]; }
void meshFacei(label val) { (*this)[3] = val; }
//- Return the real patchId
label realPatchi() const
{
const label id = patchi();
return (id < 0 ? -(id + 1) : id);
}
//- Set patchId as finiteArea
void faPatchi(label val)
{
patchi(-(val + 1));
}
//- Considered to be finiteArea if patchi < 0
bool is_finiteArea() const noexcept
{
return (patchi() < 0);
}
//- Considered to be processor local
bool is_localProc() const noexcept
{
return (procNo() == Pstream::myProcNo());
}
};
// Private Data
//- Face labels
@ -131,7 +235,10 @@ class faMesh
// Demand-driven data
//- Primitive patch
mutable uindirectPrimitivePatch* patchPtr_;
mutable std::unique_ptr<uindirectPrimitivePatch> patchPtr_;
//- List of proc/mesh-face for boundary edge neighbours
mutable std::unique_ptr<List<labelPair>> bndConnectPtr_;
//- Ldu addressing data
mutable faMeshLduAddressing* lduPtr_;
@ -172,8 +279,8 @@ class faMesh
//- Edge area normals
mutable edgeVectorField* edgeAreaNormalsPtr_;
//- Edge area normals
mutable vectorField* pointAreaNormalsPtr_;
//- Point area normals
mutable std::unique_ptr<vectorField> pointAreaNormalsPtr_;
//- Face curvatures
mutable areaScalarField* faceCurvaturesPtr_;
@ -185,14 +292,26 @@ class faMesh
mutable boolList* correctPatchPointNormalsPtr_;
// Other mesh-related data
// Other mesh-related data
//- Parallel info
mutable autoPtr<faGlobalMeshData> globalMeshDataPtr_;
//- Mapping/swapping for boundary edge halo neighbours
mutable std::unique_ptr<faMeshBoundaryHalo> haloMapPtr_;
//- Face centres for boundary edge halo neighbours
mutable std::unique_ptr<pointField> haloFaceCentresPtr_;
//- Face normals for boundary edge halo neighbours
mutable std::unique_ptr<vectorField> haloFaceNormalsPtr_;
// Static Private Data
//- Use the original method for point normals
static int origPointAreaMethod_;
//- Use quadrics fit
static const int quadricsFit_;
@ -211,6 +330,20 @@ class faMesh
//- Set primitive mesh data
void setPrimitiveMeshData();
//- Get list of (proc/patchi/patchEdgei/meshFacei) tuple pairs in an
//- globally consistent ordering
List<Pair<patchTuple>> getBoundaryEdgeConnections() const;
//- Determine the boundary edge neighbour connections
void calcBoundaryConnections() const;
//- Define boundary edge neighbours (proc/face) based on
//- gathered topology information
void setBoundaryConnections
(
const List<Pair<patchTuple>>& bndEdgeConnections
) const;
// Private member functions to calculate demand driven data
@ -242,10 +375,13 @@ class faMesh
void calcEdgeAreaNormals() const;
//- Calculate point area normals
void calcPointAreaNormals() const;
void calcPointAreaNormals_orig(vectorField& result) const;
//- Calculate point area normals
void calcPointAreaNormals(vectorField& result) const;
//- Calculate point area normals by quadrics fit
void calcPointAreaNormalsByQuadricsFit() const;
void calcPointAreaNormalsByQuadricsFit(vectorField& result) const;
//- Calculate face curvatures
void calcFaceCurvatures() const;
@ -256,6 +392,9 @@ class faMesh
//- Clear geometry but not the face areas
void clearGeomNotAreas() const;
//- Clear boundary halo information
void clearHalo() const;
//- Clear geometry
void clearGeom() const;
@ -266,10 +405,13 @@ class faMesh
void clearOut() const;
// Helpers
// Halo handling
//- Get the polyPatch pairs for the boundary edges (natural order)
List<LabelledItem<edge>> getBoundaryEdgePatchPairs() const;
//- Calculate halo centres/normals
void calcHaloFaceGeometry() const;
// Helpers
//- Create a single patch
PtrList<faPatch> createOnePatch
@ -286,13 +428,30 @@ class faMesh
const dictionary* defaultPatchDefinition = nullptr
) const;
//- Reorder processor edges using order of the
//- neighbour processorPolyPatch
void reorderProcEdges
//- Fatal error if edge labels are out of range
void checkBoundaryEdgeLabelRange(const labelUList& edgeLabels) const;
//- Extract list from contiguous (unordered) boundary data
//- to the locally sorted order.
template<class T>
List<T> boundarySubset
(
faPatchData& patchDef,
const List<LabelledItem<edge>>& bndEdgePatchPairs
) const;
const UList<T>& bndField,
const labelUList& edgeLabels
) const
{
#ifdef FULLDEBUG
checkBoundaryEdgeLabelRange(edgeLabels);
#endif
// Like UIndirectList but with an offset
List<T> result(edgeLabels.size());
forAll(edgeLabels, i)
{
result[i] = bndField[edgeLabels[i] - nInternalEdges_];
}
return result;
}
public:
@ -451,7 +610,7 @@ public:
//- Return constant reference to boundary mesh
inline const faBoundaryMesh& boundary() const noexcept;
//- Return faMesh face labels
//- Return the underlying polyMesh face labels
inline const labelList& faceLabels() const noexcept;
@ -486,6 +645,33 @@ public:
return (edgeIndex < nInternalEdges_);
}
//- List of proc/face for the boundary edge neighbours
//- using primitive patch edge numbering.
inline const List<labelPair>& boundaryConnections() const;
//- Boundary edge neighbour processors
//- (does not include own proc)
labelList boundaryProcs() const;
//- List of proc/size for the boundary edge neighbour processors
//- (does not include own proc)
List<labelPair> boundaryProcSizes() const;
//- Mapping/swapping for boundary halo neighbours
const faMeshBoundaryHalo& boundaryHaloMap() const;
//- Face centres of boundary halo neighbours
const pointField& haloFaceCentres() const;
//- Face normals of boundary halo neighbours
const vectorField& haloFaceNormals() const;
//- Face centres of boundary halo neighbours for specified patch
tmp<pointField> haloFaceCentres(const label patchi) const;
//- Face normals of boundary halo neighbours for specified patch
tmp<vectorField> haloFaceNormals(const label patchi) const;
// Mesh motion and morphing
@ -590,6 +776,7 @@ public:
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "faMeshI.H"

194
src/finiteArea/faMesh/faMeshBoundaryHalo.C Normal file
View File

@ -0,0 +1,194 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "faMeshBoundaryHalo.H"
#include "faMesh.H"
#include "globalIndex.H"
#include "Pstream.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(faMeshBoundaryHalo, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::faMeshBoundaryHalo::faMeshBoundaryHalo(const label comm)
:
mapDistributeBase(comm),
inputMeshFaces_(),
boundaryToCompact_()
{}
Foam::faMeshBoundaryHalo::faMeshBoundaryHalo(const faMesh& areaMesh)
:
mapDistributeBase(),
inputMeshFaces_(),
boundaryToCompact_()
{
reset(areaMesh);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::faMeshBoundaryHalo::clear()
{
static_cast<mapDistributeBase&>(*this) = mapDistributeBase();
inputMeshFaces_.clear();
boundaryToCompact_.clear();
}
Foam::label Foam::faMeshBoundaryHalo::haloSize() const
{
if (Pstream::parRun())
{
return boundaryToCompact_.size();
}
else
{
return inputMeshFaces_.size();
}
}
void Foam::faMeshBoundaryHalo::reset(const faMesh& areaMesh)
{
inputMeshFaces_.clear();
boundaryToCompact_.clear();
const auto& procConnections = areaMesh.boundaryConnections();
if (!Pstream::parRun())
{
// Serial - extract halo numbers directly
inputMeshFaces_.resize(procConnections.size());
forAll(procConnections, connecti)
{
// Connected neighbour, non-parallel = must be local
const auto& tuple = procConnections[connecti];
// const label nbrProci = tuple.first();
const label nbrFacei = tuple.second();
inputMeshFaces_[connecti] = nbrFacei;
}
return;
}
const label nProcs = Pstream::nProcs(comm_);
const label myRank = Pstream::myProcNo(comm_);
const globalIndex globalFaceNum(areaMesh.mesh().nFaces());
// Boundary inside faces in polyMesh face ids
const labelList insideFaces
(
UIndirectList<label>
(
areaMesh.faceLabels(),
areaMesh.patch().boundaryFaces()
)
);
// Slightly circuitous, but allows maximum reuse of mapDistributeBase
// 1. Construct a connectivity map using global face numbers
labelListList connectivity(areaMesh.nBoundaryEdges());
List<Map<label>> compactMap(nProcs, Map<label>(0));
// All local mesh faces used
labelHashSet localUsed(insideFaces);
forAll(connectivity, connecti)
{
labelList& edgeFaces = connectivity[connecti];
edgeFaces.resize(2);
// Owner is the boundary inside face (our side)
// Neighbour is the boundary outside face
// Connected neighbour
const auto& tuple = procConnections[connecti];
const label nbrProci = tuple.first();
const label nbrFacei = tuple.second();
if (myRank == nbrProci)
{
// Processor-local connectivity
localUsed.insert(nbrFacei);
}
// Global addressing for the connectivity
edgeFaces[0] = globalFaceNum.toGlobal(insideFaces[connecti]);
edgeFaces[1] = globalFaceNum.toGlobal(nbrProci, nbrFacei);
}
// Create and replace mapping
static_cast<mapDistributeBase&>(*this) = mapDistributeBase
(
globalFaceNum,
connectivity,
compactMap,
Pstream::msgType(),
comm_
);
// List of local mesh faces referenced.
// Includes inside and locally connected outside faces
inputMeshFaces_ = localUsed.sortedToc();
boundaryToCompact_.clear();
boundaryToCompact_.resize(connectivity.size());
// After creating the map, connectivity is localized *and*
// uses compact numbering!
// Extract the neighbour connection (compact numbering)
forAll(connectivity, connecti)
{
const labelList& edgeFaces = connectivity[connecti];
// const label face0 = edgeFaces[0];
const label face1 = edgeFaces[1];
boundaryToCompact_[connecti] = face1;
}
}
// ************************************************************************* //

154
src/finiteArea/faMesh/faMeshBoundaryHalo.H Normal file
View File

@ -0,0 +1,154 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::faMeshBoundaryHalo
Description
Class for obtaining halo face data for the boundary edges.
The ordering follows that natural edge ordering of the underlying
primitive patch.
Note
The halo faces can be located on-processor or off-processor.
SourceFiles
faMeshBoundaryHalo.C
faMeshBoundaryHaloTemplates.C
\*---------------------------------------------------------------------------*/
#ifndef faMeshBoundaryHalo_H
#define faMeshBoundaryHalo_H
#include "mapDistributeBase.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward Declarations
class faMesh;
/*---------------------------------------------------------------------------*\
Class faMeshBoundaryHalo Declaration
\*---------------------------------------------------------------------------*/
class faMeshBoundaryHalo
:
public mapDistributeBase
{
// Private Data
//- List of local input mesh faces required
labelList inputMeshFaces_;
//- Internal mapping from boundary index to compact
labelList boundaryToCompact_;
public:
// Declare name of the class and its debug switch
ClassName("faMeshBoundaryHalo");
// Constructors
//- Default construct
explicit faMeshBoundaryHalo(const label comm = UPstream::worldComm);
//- Construct from mesh
explicit faMeshBoundaryHalo(const faMesh& mesh);
// Member Functions
//- Clear out all parameters
void clear();
//- Redefine map and connectivity for a mesh
void reset(const faMesh& mesh);
//- The local data size (output)
label haloSize() const;
//- List of local input mesh faces required.
// \note will not correspond exactly to the boundary inside faces.
// Duplicates have been removed and it also contains the
// processor-local neighbour faces, which would otherwise not be
// handled by the distribute method.
const labelList& inputMeshFaces() const noexcept
{
return inputMeshFaces_;
}
// Other
//- Distribute sparse data.
// On output it is adjusted.
template<class Type>
void distributeSparse
(
List<Type>& fld,
const labelUList& sparseInputLocations,
const labelUList& compactOutputMapping
) const;
//- Distribute sparse data.
// On output it is adjusted.
template<class Type>
void distributeSparse
(
List<Type>& fld,
const labelUList& sparseInputLocations
) const;
//- Distribute sparse data.
// The input field one enty per sparse id (inputMeshFaces).
// On output it will have for the input sparse
// The input field contains location.
template<class Type>
void distributeSparse(List<Type>& fld) const;
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
#include "faMeshBoundaryHaloTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

107
src/finiteArea/faMesh/faMeshBoundaryHaloTemplates.C Normal file
View File

@ -0,0 +1,107 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "faMeshBoundaryHalo.H"
#include "UIndirectList.H"
#include "Pstream.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template<class Type>
void Foam::faMeshBoundaryHalo::distributeSparse
(
List<Type>& fld,
const labelUList& sparseInputLocations,
const labelUList& compactOutputMapping
) const
{
if (!Pstream::parRun())
{
return; // No-op in serial
}
// Construct data in compact addressing
List<Type> compactFld(constructSize_, Zero);
if (sparseInputLocations.empty())
{
// Copy in as dense field
forAll(fld, i)
{
compactFld[i] = fld[i];
}
}
else
{
if (fld.size() != sparseInputLocations.size())
{
FatalErrorInFunction
<< "Input field size (" << fld.size()
<< " != sparse ids size ("
<< sparseInputLocations.size() << ")\n"
<< exit(FatalError);
}
// Copy in sparse locations
forAll(sparseInputLocations, i)
{
const label idx = sparseInputLocations[i];
if (idx != -1)
{
compactFld[idx] = fld[i];
}
}
}
// Pull all data
mapDistributeBase::distribute<Type>(compactFld);
// Rewrite to output
fld = UIndirectList<Type>(compactFld, compactOutputMapping);
}
template<class Type>
void Foam::faMeshBoundaryHalo::distributeSparse
(
List<Type>& fld,
const labelUList& sparseInputLocations
) const
{
this->distributeSparse(fld, sparseInputLocations, boundaryToCompact_);
}
template<class Type>
void Foam::faMeshBoundaryHalo::distributeSparse(List<Type>& fld) const
{
this->distributeSparse(fld, inputMeshFaces_, boundaryToCompact_);
}
// ************************************************************************* //

540
src/finiteArea/faMesh/faMeshDemandDrivenData.C
View File

@ -40,9 +40,73 @@ License
#include "processorFaPatchFields.H"
#include "emptyFaPatchFields.H"
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
// Define an area-weighted normal for three points (defining a triangle)
// (p0, p1, p2) are the base, first, second points respectively
//
// From the original Tukovic code:
//
// vector n = (d1^d2)/mag(d1^d2);
// scalar sinAlpha = mag(d1^d2)/(mag(d1)*mag(d2));
// scalar w = sinAlpha/(mag(d1)*mag(d2));
//
// ie, normal weighted by area, sine angle and inverse distance squared.
// - area : larger weight for larger areas
// - sin : lower weight for narrow angles (eg, shards)
// - inv distance squared : lower weights for distant points
//
// The above refactored, with 0.5 for area:
//
// (d1 ^ d2) / (2 * magSqr(d1) * magSqr(d2))
static inline vector areaInvDistSqrWeightedNormal
(
const vector& a,
const vector& b
)
{
const scalar s(2*magSqr(a)*magSqr(b));
return s < VSMALL ? Zero : (a ^ b) / s;
}
// The area normal for the face dual (around base-point)
// described by the right/left edge points and the centre point
//
// The adjustment for 1/2 edge point (the dual point) is done internally
static inline vector areaInvDistSqrWeightedNormalDualEdge
(
const point& basePoint,
const point& rightPoint,
const point& leftPoint,
const point& centrePoint
)
{
const vector mid(centrePoint - basePoint);
return
(
areaInvDistSqrWeightedNormal
(
0.5*(rightPoint - basePoint), // vector to right 1/2 edge
mid
)
+ areaInvDistSqrWeightedNormal
(
mid,
0.5*(leftPoint - basePoint) // vector to left 1/2 edge
)
);
}
} // End namespace Foam
namespace Foam
{
@ -61,6 +125,7 @@ static Foam::bitSet markupBoundaryPoints(const uindirectPrimitivePatch& p)
return markPoints;
}
} // End namespace Foam
@ -496,108 +561,40 @@ void Foam::faMesh::calcEdgeAreaNormals() const
// Point area normals
const vectorField& pointNormals = pointAreaNormals();
// // Primitive patch edge normals
// const labelListList& patchPointEdges = patch().pointEdges();
// vectorField patchEdgeNormals(nEdges(), Zero);
// forAll(pointNormals, pointI)
// {
// const labelList& curPointEdges = patchPointEdges[pointI];
// forAll(curPointEdges, edgeI)
// {
// label curEdge = curPointEdges[edgeI];
// patchEdgeNormals[curEdge] += 0.5*pointNormals[pointI];
// }
// }
// patchEdgeNormals /= mag(patchEdgeNormals);
// // Edge area normals
// label nIntEdges = patch().nInternalEdges();
// for (label edgeI = 0; edgeI < nIntEdges; ++edgeI)
// {
// edgeAreaNormals.ref()[edgeI] =
// patchEdgeNormals[edgeI];
// }
// forAll(boundary(), patchI)
// {
// const labelList& edgeLabels = boundary()[patchI];
// forAll(edgeAreaNormals.boundaryFieldRef()[patchI], edgeI)
// {
// edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] =
// patchEdgeNormals[edgeLabels[edgeI]];
// }
// }
forAll(edgeAreaNormals.internalField(), edgeI)
forAll(edgeAreaNormals.internalField(), edgei)
{
const vector e = edges()[edgeI].unitVec(points());
const edge& e = edges()[edgei];
const vector edgeVec = e.unitVec(points());
// scalar wStart =
// 1.0 - sqr(mag(e^pointNormals[edges()[edgeI].end()]));
vector& n = edgeAreaNormals.ref()[edgei];
// scalar wEnd =
// 1.0 - sqr(mag(e^pointNormals[edges()[edgeI].start()]));
n = (pointNormals[e.first()] + pointNormals[e.second()]);
// wStart = 1.0;
// wEnd = 1.0;
n -= edgeVec*(edgeVec & n);
// edgeAreaNormals.ref()[edgeI] =
// wStart*pointNormals[edges()[edgeI].start()]
// + wEnd*pointNormals[edges()[edgeI].end()];
// vector eC =
// 0.5
// *(
// points()[edges()[edgeI].start()]
// + points()[edges()[edgeI].end()]
// );
// vector eCp = 0.5*
// (
// points()[edges()[edgeI].start()]
// + pointNormals[edges()[edgeI].start()]
// points()[edges()[edgeI].end()] +
// );
edgeAreaNormals.ref()[edgeI] =
pointNormals[edges()[edgeI].start()]
+ pointNormals[edges()[edgeI].end()];
edgeAreaNormals.ref()[edgeI] -=
e*(e&edgeAreaNormals.internalField()[edgeI]);
n.normalise();
}
edgeAreaNormals.ref() /= mag(edgeAreaNormals.internalField());
forAll(boundary(), patchI)
forAll(boundary(), patchi)
{
const edgeList::subList patchEdges =
boundary()[patchI].patchSlice(edges());
boundary()[patchi].patchSlice(edges());
forAll(patchEdges, edgeI)
vectorField& edgeNorms = edgeAreaNormals.boundaryFieldRef()[patchi];
forAll(patchEdges, edgei)
{
edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] =
pointNormals[patchEdges[edgeI].start()]
+ pointNormals[patchEdges[edgeI].end()];
const edge& e = patchEdges[edgei];
const vector edgeVec = e.unitVec(points());
const vector e = patchEdges[edgeI].unitVec(points());
vector& n = edgeNorms[edgei];
edgeAreaNormals.boundaryFieldRef()[patchI][edgeI] -=
e*(e&edgeAreaNormals.boundaryField()[patchI][edgeI]);
n = (pointNormals[e.first()] + pointNormals[e.second()]);
n -= edgeVec*(edgeVec & n);
n.normalise();
}
edgeAreaNormals.boundaryFieldRef()[patchI] /=
mag(edgeAreaNormals.boundaryField()[patchI]);
}
}
@ -880,19 +877,53 @@ Foam::labelList Foam::faMesh::boundaryPoints() const
}
void Foam::faMesh::calcPointAreaNormals() const
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Point normal calculations
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Original method (general)
// -------------------------
// - For each point, obtain the list of connected patch faces
// (from point-to-face addressing).
//
// - Create a primitive patch for those faces and use that to obtain the
// outer edge loop(s). This is effectively an agglomeration of the patch
// faces connected to a point.
//
// - Perform a pair-wise walk of the edge loop to obtain triangles from
// the originating point outwards (fan-like triangulation).
// Calculate an area-weighted value for each triangle.
//
// NOTE: not sure why internal and boundary point agglomeration was
// handled separately.
//
// Problems:
// - possibly susceptible to edge-loop errors (issue #2233) that cause
// the agglomeration logic to include the current point twice?
// - use of outer edge loop makes it more sensitive to face warpage.
// - relatively expensive with point-face connectivity,
// creation/destruction of a primitive-patch around each point.
//
// Original method (boundary correction)
// -------------------------------------
//
// - correct wedge directly, use processor patch information to exchange
// the current summed values
//
// - explicit correction of other boundaries.
// Polls the patch for the ngbPolyPatchPointNormals(), which internally
// calls ngbPolyPatchFaces and can return -1 for unmatched edges.
// This occurs when the outside perimeter of the faPatch aligns with
// a polyMesh processor. The neighbour face is off-processor and cannot
// be found. Accessing the mesh face at -1 == SEGFAULT.
void Foam::faMesh::calcPointAreaNormals_orig(vectorField& result) const
{
if (pointAreaNormalsPtr_)
{
FatalErrorInFunction
<< "pointAreaNormalsPtr_ already allocated"
<< abort(FatalError);
}
DebugInFunction
<< "Calculating pointAreaNormals : original method" << endl;
pointAreaNormalsPtr_ = new vectorField(nPoints(), Zero);
vectorField& result = *pointAreaNormalsPtr_;
result.resize(nPoints());
result = Zero;
labelList intPoints(internalPoints());
labelList bndPoints(boundaryPoints());
@ -1004,14 +1035,14 @@ void Foam::faMesh::calcPointAreaNormals() const
const labelList& patchPoints = wedgePatch.pointLabels();
vector N =
const vector N
(
transform
(
wedgePatch.edgeT(),
wedgePatch.centreNormal()
);
N /= mag(N);
).normalise()
);
for (const label pti : patchPoints)
{
@ -1157,14 +1188,302 @@ void Foam::faMesh::calcPointAreaNormals() const
}
}
result /= mag(result);
for (vector& n : result)
{
n.normalise();
}
}
void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const
{
vectorField& result = *pointAreaNormalsPtr_;
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Point normal calculations
// ~~~~~~~~~~~~~~~~~~~~~~~~~
// Revised method (general)
// ------------------------
//
// - For each patch face obtain a centre point (mathematical avg)
// and use that to define the face dual as a pair of triangles:
// - tri1: base-point / mid-side of right edge / face centre
// - tri2: base-point / face centre / mid-side of left edge
//
// - Walk all face points, inserting directly into the corresponding
// locations. No distinction between internal or boundary points (yet).
//
// Revised method (boundary correction)
// ------------------------------------
//
// - correct wedge directly, use processor patch information to exchange
// the current summed values. [No change from original].
//
// - explicit correction of other boundaries.
// Use the new boundary halo information for the face normals.
// Calculate the equivalent face-point normals locally and apply
// correction as before.
void Foam::faMesh::calcPointAreaNormals(vectorField& result) const
{
DebugInFunction
<< "Calculating pointAreaNormals : face dual method" << endl;
result.resize(nPoints());
result = Zero;
const pointField& points = patch().localPoints();
const faceList& faces = patch().localFaces();
// Loop over all faces
for (const face& f : faces)
{
const label nVerts(f.size());
point centrePoint(Zero);
for (label i = 0; i < nVerts; ++i)
{
centrePoint += points[f[i]];
}
centrePoint /= nVerts;
for (label i = 0; i < nVerts; ++i)
{
const label pt0 = f.thisLabel(i); // base
result[pt0] +=
areaInvDistSqrWeightedNormalDualEdge
(
points[pt0], // base
points[f.nextLabel(i)], // right
points[f.prevLabel(i)], // left
centrePoint
);
}
}
// Handle the boundary edges
bitSet nbrBoundaryAdjust(boundary().size(), true);
forAll(boundary(), patchi)
{
const faPatch& fap = boundary()[patchi];
if (isA<wedgeFaPatch>(fap))
{
// Correct wedge points
const auto& wedgePatch = refCast<const wedgeFaPatch>(fap);
const labelList& patchPoints = wedgePatch.pointLabels();
const vector N
(
transform
(
wedgePatch.edgeT(),
wedgePatch.centreNormal()
).normalise()
);
for (const label pti : patchPoints)
{
result[pti] -= N*(N&result[pti]);
}
// Axis point correction
if (wedgePatch.axisPoint() > -1)
{
result[wedgePatch.axisPoint()] =
wedgePatch.axis()
*(
wedgePatch.axis()
&result[wedgePatch.axisPoint()]
);
}
// Handled
nbrBoundaryAdjust.unset(patchi);
}
else if (Pstream::parRun() && isA<processorFaPatch>(fap))
{
// Correct processor patch points
const auto& procPatch = refCast<const processorFaPatch>(fap);
const labelList& patchPoints = procPatch.pointLabels();
const labelList& nbrPatchPoints = procPatch.neighbPoints();
const labelList& nonGlobalPatchPoints =
procPatch.nonGlobalPatchPoints();
// Send my values
vectorField patchPointNormals
(
UIndirectList<vector>(result, patchPoints)
);
{
OPstream::write
(
Pstream::commsTypes::blocking,
procPatch.neighbProcNo(),
reinterpret_cast<const char*>(patchPointNormals.cdata()),
patchPointNormals.size_bytes()
);
}
// Receive neighbour values
patchPointNormals.resize(nbrPatchPoints.size());
{
IPstream::read
(
Pstream::commsTypes::blocking,
procPatch.neighbProcNo(),
reinterpret_cast<char*>(patchPointNormals.data()),
patchPointNormals.size_bytes()
);
}
for (const label pti : nonGlobalPatchPoints)
{
result[patchPoints[pti]] +=
patchPointNormals[nbrPatchPoints[pti]];
}
// Handled
nbrBoundaryAdjust.unset(patchi);
}
else if (fap.coupled())
{
// Coupled - no further action for neighbour side
nbrBoundaryAdjust.unset(patchi);
}
// TBD:
/// else if (isA<emptyFaPatch>(fap))
/// {
/// // Ignore this boundary
/// nbrBoundaryAdjust.unset(patchi);
/// }
else if (!correctPatchPointNormals(patchi))
{
// No corrections
nbrBoundaryAdjust.unset(patchi);
}
}
// Correct global points
if (globalData().nGlobalPoints())
{
const labelList& spLabels = globalData().sharedPointLabels();
const labelList& addr = globalData().sharedPointAddr();
vectorField spNormals
(
UIndirectList<vector>(result, spLabels)
);
vectorField gpNormals
(
globalData().nGlobalPoints(),
Zero
);
forAll(addr, i)
{
gpNormals[addr[i]] += spNormals[i];
}
combineReduce(gpNormals, plusEqOp<vectorField>());
// Extract local data
forAll(addr, i)
{
spNormals[i] = gpNormals[addr[i]];
}
forAll(spNormals, pointI)
{
result[spLabels[pointI]] = spNormals[pointI];
}
}
if (returnReduce(nbrBoundaryAdjust.any(), orOp<bool>()))
{
if (debug)
{
PoutInFunction
<< "Apply " << nbrBoundaryAdjust.count()
<< " boundary neighbour corrections" << nl;
}
// Apply boundary points correction
// Collect face normals as point normals
const auto& haloNormals = this->haloFaceNormals();
Map<vector> fpNormals(4*nBoundaryEdges());
for (const label patchi : nbrBoundaryAdjust)
{
const faPatch& fap = boundary()[patchi];
const labelList& edgeLabels = fap.edgeLabels();
if (fap.ngbPolyPatchIndex() < 0)
{
FatalErrorInFunction
<< "Neighbour polyPatch index is not defined "
<< "for faPatch " << fap.name()
<< abort(FatalError);
}
for (const label edgei : edgeLabels)
{
const edge& e = patch().edges()[edgei];
// Halo face unitNormal at boundary edge (starts as 0)
const vector& fnorm = haloNormals[edgei - nInternalEdges_];
fpNormals(e.first()) += fnorm;
fpNormals(e.second()) += fnorm;
}
}
// Apply the correction
// Note from Zeljko Tukovic:
//
// This posibility is used for free-surface tracking
// calculations to enforce 90 deg contact angle between
// finite-area mesh and neighbouring polyPatch. It is very
// important for curvature calculation to have correct normal
// at contact line points.
forAllConstIters(fpNormals, iter)
{
const label pointi = iter.key();
vector fpnorm = iter.val();
fpnorm.normalise();
result[pointi] -= fpnorm*(fpnorm & result[pointi]);
}
}
for (vector& n : result)
{
n.normalise();
}
}
void Foam::faMesh::calcPointAreaNormalsByQuadricsFit(vectorField& result) const
{
const labelList intPoints(internalPoints());
const labelList bndPoints(boundaryPoints());
@ -1724,7 +2043,10 @@ void Foam::faMesh::calcPointAreaNormalsByQuadricsFit() const
}
}
result /= mag(result);
for (vector& n : result)
{
n.normalise();
}
}

11
src/finiteArea/faMesh/faMeshI.H
View File

@ -139,4 +139,15 @@ inline Foam::uindirectPrimitivePatch& Foam::faMesh::patch()
}
inline const Foam::List<Foam::labelPair>&
Foam::faMesh::boundaryConnections() const
{
if (!bndConnectPtr_)
{
calcBoundaryConnections();
}
return *bndConnectPtr_;
}
// ************************************************************************* //

1017
src/finiteArea/faMesh/faMeshPatches.C
View File

File diff suppressed because it is too large Load Diff

930
src/finiteArea/faMesh/faMeshTopology.C Normal file
View File

@ -0,0 +1,930 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software: you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "faMesh.H"
#include "faMeshBoundaryHalo.H"
#include "globalMeshData.H"
#include "indirectPrimitivePatch.H"
#include "edgeHashes.H"
#include "foamVtkLineWriter.H"
#include "foamVtkIndPatchWriter.H"
#include "foamVtkUIndPatchWriter.H"
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
// Print out edges as point pairs
template<class PatchType>
static void printPatchEdges
(
Ostream& os,
const PatchType& p,
const labelList& edgeIds,
label maxOutput = 10
)
{
label nOutput = 0;
for (const label patchEdgei : edgeIds)
{
const edge e(p.meshEdge(patchEdgei));
os << " "
<< p.points()[e.first()] << ' '
<< p.points()[e.second()] << nl;
++nOutput;
if (maxOutput > 0 && nOutput >= maxOutput)
{
os << " ... suppressing further output" << nl;
break;
}
}
}
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::List<Foam::Pair<Foam::faMesh::patchTuple>>
Foam::faMesh::getBoundaryEdgeConnections() const
{
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
const label nNonProcessor = pbm.nNonProcessor();
const label nInternalEdges = patch().nInternalEdges();
const label nBoundaryEdges = patch().nBoundaryEdges();
// The output result:
List<Pair<patchTuple>> bndEdgeConnections(nBoundaryEdges);
// Map edges (mesh numbering) back to a boundary index
EdgeMap<label> edgeToBoundaryIndex(2*nBoundaryEdges);
label nBadEdges(0);
labelHashSet badEdges(2*nBoundaryEdges);
{
// Local collection structure for accounting of patch pairs.
// Based on 'edge' which has some hash-like insertion properties
// that are useful.
struct patchPairingType : public Foam::edge
{
label patchEdgei_ = -1;
label meshFacei_ = -1;
void clear()
{
Foam::edge::clear();
patchEdgei_ = -1;
meshFacei_ = -1;
}
};
List<patchPairingType> patchPairings(nBoundaryEdges);
DebugInFunction
<< "Determining required boundary edge connections, "
<< "resolving locally attached boundary edges." << endl;
// Pass 1:
// - setup lookup (edge -> bnd index)
// - add owner patch for each boundary edge
for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei)
{
const label patchEdgei = (bndEdgei + nInternalEdges);
edgeToBoundaryIndex.insert
(
patch().meshEdge(patchEdgei),
bndEdgei
);
// The attached patch face. Should only be one!
const labelList& edgeFaces = patch().edgeFaces()[patchEdgei];
if (edgeFaces.size() != 1)
{
badEdges.insert(patchEdgei);
continue;
}
const label patchFacei = edgeFaces[0];
const label meshFacei = faceLabels_[patchFacei];
const label bndFacei = (meshFacei - mesh().nInternalFaces());
/// const label patchId = pbm.whichPatch(meshFacei);
const label patchId = pbm.patchID()[bndFacei];
// Primary bookkeeping
{
auto& tuple = bndEdgeConnections[bndEdgei].first();
tuple.procNo(Pstream::myProcNo());
tuple.faPatchi(patchId); // Tag as finiteArea patch
tuple.patchEdgei(patchEdgei);
tuple.meshFacei(meshFacei);
}
// Temporary local bookkeeping
{
auto& pairing = patchPairings[bndEdgei];
pairing.clear(); // invalidate
pairing.insert(patchId); // 'hash' into first location
pairing.patchEdgei_ = patchEdgei;
pairing.meshFacei_ = meshFacei;
}
}
if ((nBadEdges = returnReduce(badEdges.size(), sumOp<label>())) > 0)
{
edgeList dumpEdges(patch().edges(), badEdges.sortedToc());
vtk::lineWriter writer
(
patch().localPoints(),
dumpEdges,
fileName
(
mesh().time().globalPath()
/ ("faMesh-construct.nonManifoldEdges")
)
);
writer.writeGeometry();
// CellData
writer.beginCellData();
writer.writeProcIDs();
FatalErrorInFunction
<< "Boundary edges not singly connected: "
<< nBadEdges << '/' << nBoundaryEdges << nl;
printPatchEdges
(
FatalError,
patch(),
badEdges.sortedToc()
);
InfoInFunction
<< "(debug) wrote " << writer.output().name() << nl;
FatalError << abort(FatalError);
}
badEdges.clear();
// Pass 2:
// Add in first connecting neighbour patch for the boundary edges.
// Need to examine all possibly connecting (non-processor) neighbours,
// but only need to check their boundary edges.
label nMissing = patchPairings.size();
for (label patchi = 0; patchi < nNonProcessor; ++patchi)
{
if (!nMissing) break; // Early exit
const polyPatch& pp = pbm[patchi];
// Check boundary edges
for
(
label patchEdgei = pp.nInternalEdges();
patchEdgei < pp.nEdges();
++patchEdgei
)
{
const label bndEdgei =
edgeToBoundaryIndex.lookup(pp.meshEdge(patchEdgei), -1);
if (bndEdgei != -1)
{
// Has a matching owner boundary edge
auto& pairing = patchPairings[bndEdgei];
// Add neighbour (patchId, patchEdgei, meshFacei)
// 'hash' into the second location
// which does not insert the same value twice
if (pairing.insert(patchi))
{
// The attached patch face. Should only be one!
const labelList& edgeFaces = pp.edgeFaces()[patchEdgei];
if (edgeFaces.size() != 1)
{
pairing.erase(patchi);
badEdges.insert(badEdges.size());
continue;
}
const label patchFacei = edgeFaces[0];
const label meshFacei = patchFacei + pp.start();
// The neighbour information
pairing.patchEdgei_ = patchEdgei;
pairing.meshFacei_ = meshFacei;
--nMissing;
if (!nMissing) break; // Early exit
}
}
}
}
if ((nBadEdges = returnReduce(badEdges.size(), sumOp<label>())) > 0)
{
FatalErrorInFunction
<< "Had " << nBadEdges
<< " boundary edges with missing or multiple edge connections"
<< abort(FatalError);
}
// Combine local bookkeeping into final list
badEdges.clear();
for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei)
{
const auto& pairing = patchPairings[bndEdgei];
const label nbrPatchi = pairing.second();
const label nbrPatchEdgei = pairing.patchEdgei_;
const label nbrMeshFacei = pairing.meshFacei_;
if (nbrMeshFacei >= 0)
{
// Add into primary bookkeeping
auto& tuple = bndEdgeConnections[bndEdgei].second();
tuple.procNo(Pstream::myProcNo());
tuple.patchi(nbrPatchi);
tuple.patchEdgei(nbrPatchEdgei);
tuple.meshFacei(nbrMeshFacei);
}
else if (!Pstream::parRun())
{
badEdges.insert(nInternalEdges + bndEdgei);
}
}
}
// ~~~~~~
// Serial - can return already
// ~~~~~~
if (!Pstream::parRun())
{
// Verbose report of missing edges - in serial
nBadEdges = badEdges.size();
if (nBadEdges)
{
edgeList dumpEdges(patch().edges(), badEdges.sortedToc());
vtk::lineWriter writer
(
patch().localPoints(),
dumpEdges,
fileName
(
mesh().time().globalPath()
/ ("faMesh-construct.invalidEdges")
)
);
writer.writeGeometry();
// CellData
writer.beginCellData();
writer.writeProcIDs();
FatalErrorInFunction
<< "Boundary edges with missing/invalid neighbours: "
<< nBadEdges << '/' << nBoundaryEdges << nl;
printPatchEdges
(
FatalError,
patch(),
badEdges.sortedToc()
);
InfoInFunction
<< "(debug) wrote " << writer.output().name() << nl;
FatalError << abort(FatalError);
}
// Globally consistent ordering
patchTuple::sort(bndEdgeConnections);
return bndEdgeConnections;
}
// ~~~~~~~~
// Parallel
// ~~~~~~~~
DebugInFunction
<< "Creating global coupling data" << endl;
const globalMeshData& globalData = mesh().globalData();
const indirectPrimitivePatch& cpp = globalData.coupledPatch();
const mapDistribute& map = globalData.globalEdgeSlavesMap();
const label nCoupledEdges = cpp.nEdges();
// Construct coupled edge usage with all data
List<bool> coupledEdgesUsed(map.constructSize(), false);
// Markup finiteArea boundary edges that are coupled across processors
for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei)
{
coupledEdgesUsed[cppEdgei] =
edgeToBoundaryIndex.found(cpp.meshEdge(cppEdgei));
}
DebugInFunction
<< "Starting sync of boundary edge topology" << endl;
globalMeshData::syncData
(
coupledEdgesUsed,
globalData.globalEdgeSlaves(),
globalData.globalEdgeTransformedSlaves(), // probably not used
map,
orEqOp<bool>()
);
if (debug)
{
label nAttached = 0;
for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei)
{
if (coupledEdgesUsed[cppEdgei])
{
++nAttached;
}
}
InfoInFunction
<< "Approx "
<< returnReduce(nAttached, sumOp<label>())
<< " connected boundary edges (total, some duplicates)" << endl;
}
// Combine information
// Normally 0 or 2 connections
List<DynamicList<patchTuple, 2>> gatheredConnections(map.constructSize());
// Map edges (mesh numbering) back to a coupled index in use
EdgeMap<label> edgeToCoupledIndex(2*nCoupledEdges);
// Pass 1
// Look for attached boundary edges
// - boundary edges from this side go into gathered connections
// - boundary edges connected from the other side are noted for later
for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei)
{
if (coupledEdgesUsed[cppEdgei])
{
const edge meshEdge(cpp.meshEdge(cppEdgei));
const label bndEdgei =
edgeToBoundaryIndex.lookup(meshEdge, -1);
if (bndEdgei != -1)
{
// A boundary finiteEdge edge (known from this side)
auto& gathered = gatheredConnections[cppEdgei];
gathered.setCapacity(2);
gathered.resize(1);
auto& tuple = gathered.last();
tuple = bndEdgeConnections[bndEdgei].first();
}
else
{
// Boundary edge connected from the other side
// - mark for it to be added later
edgeToCoupledIndex.insert(meshEdge, cppEdgei);
}
}
}
// Pass 2
// - add previously noted boundary edges (connected from other side)
// into gathered connections
badEdges.clear();
for (label patchi = 0; patchi < nNonProcessor; ++patchi)
{
if (edgeToCoupledIndex.empty()) break; // Early exit
const polyPatch& pp = pbm[patchi];
// Check boundary edges
for
(
label patchEdgei = pp.nInternalEdges();
patchEdgei < pp.nEdges();
++patchEdgei
)
{
const edge meshEdge(pp.meshEdge(patchEdgei));
const label cppEdgei =
edgeToCoupledIndex.lookup(meshEdge, -1);
if (cppEdgei != -1)
{
// A known connection
// The attached patch face. Should only be one!
const labelList& edgeFaces = pp.edgeFaces()[patchEdgei];
if (edgeFaces.size() != 1)
{
badEdges.insert(cppEdgei);
continue;
}
const label patchFacei = edgeFaces[0];
const label meshFacei = patchFacei + pp.start();
auto& gathered = gatheredConnections[cppEdgei];
gathered.setCapacity(2);
gathered.resize(1);
auto& tuple = gathered.last();
tuple.procNo(Pstream::myProcNo());
tuple.patchi(patchi);
tuple.patchEdgei(patchEdgei);
tuple.meshFacei(meshFacei);
// Do not consider again
edgeToCoupledIndex.erase(meshEdge);
if (edgeToCoupledIndex.empty()) break; // Early exit
}
}
}
if ((nBadEdges = returnReduce(badEdges.size(), sumOp<label>())) > 0)
{
FatalErrorInFunction
<< "Had " << nBadEdges << " coupled boundary edges"
<< " with missing or multiple edge connections"
<< abort(FatalError);
}
DebugInFunction
<< "Starting sync of boundary edge information" << endl;
globalMeshData::syncData
(
gatheredConnections,
globalData.globalEdgeSlaves(),
globalData.globalEdgeTransformedSlaves(), // probably not used
map,
ListOps::appendEqOp<patchTuple>()
);
DebugInFunction
<< "Collating sync information" << endl;
// Pick out gathered connections and add into primary bookkeeping
for (label cppEdgei = 0; cppEdgei < nCoupledEdges; ++cppEdgei)
{
const auto& gathered = gatheredConnections[cppEdgei];
const label bndEdgei =
edgeToBoundaryIndex.lookup(cpp.meshEdge(cppEdgei), -1);
if
(
// A boundary finiteEdge edge (known from this side)
bndEdgei != -1
// Gathered a one-to-one connection
&& gathered.size() == 2
)
{
const auto& a = gathered[0];
const auto& b = gathered[1];
// Copy second side of connection
auto& connection = bndEdgeConnections[bndEdgei];
connection.second() = (connection.first() == b) ? a : b;
}
}
// Check missing/invalid
badEdges.clear();
for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei)
{
const auto& connection = bndEdgeConnections[bndEdgei];
if (!connection.second().valid())
{
badEdges.insert(nInternalEdges + bndEdgei);
}
}
if (debug & 8)
{
// Boundary edges
{
vtk::lineWriter writer
(
patch().localPoints(),
patch().boundaryEdges(),
fileName
(
mesh().time().globalPath()
/ ("faMesh-construct.boundaryEdges")
)
);
writer.writeGeometry();
// CellData
writer.beginCellData();
writer.writeProcIDs();
// For each boundary edge - the associate neighbour patch
labelList neighProc(nBoundaryEdges);
labelList neighPatch(nBoundaryEdges);
for (label bndEdgei = 0; bndEdgei < nBoundaryEdges; ++bndEdgei)
{
const auto& connection = bndEdgeConnections[bndEdgei];
neighProc[bndEdgei] = connection.second().procNo();
neighPatch[bndEdgei] = connection.second().patchi();
}
writer.write("neighProc", neighProc);
writer.write("neighPatch", neighPatch);
}
// finiteArea
{
vtk::uindirectPatchWriter writer
(
patch(),
fileName
(
mesh().time().globalPath()
/ ("faMesh-construct.faPatch")
)
);
writer.writeGeometry();
// CellData
writer.beginCellData();
writer.writeProcIDs();
}
}
// Verbose report of missing edges
if ((nBadEdges = returnReduce(badEdges.size(), sumOp<label>())) > 0)
{
edgeList dumpEdges(patch().edges(), badEdges.sortedToc());
vtk::lineWriter writer
(
patch().localPoints(),
dumpEdges,
fileName
(
mesh().time().globalPath()
/ ("faMesh-construct.invalidEdges")
)
);
writer.writeGeometry();
// CellData
writer.beginCellData();
writer.writeProcIDs();
FatalErrorInFunction
<< "Boundary edges with missing/invalid neighbours: "
<< nBadEdges << '/' << nBoundaryEdges << nl;
printPatchEdges
(
FatalError,
patch(),
badEdges.sortedToc()
);
InfoInFunction
<< "(debug) wrote " << writer.output().name() << nl;
FatalError << abort(FatalError);
}
// Globally consistent ordering
patchTuple::sort(bndEdgeConnections);
DebugInFunction
<< "Return sorted list of boundary connections" << endl;
return bndEdgeConnections;
}
void Foam::faMesh::setBoundaryConnections
(
const List<Pair<patchTuple>>& bndEdgeConnections
) const
{
const label nInternalEdges = patch().nInternalEdges();
const label nBoundaryEdges = patch().nBoundaryEdges();
if (bndEdgeConnections.size() != nBoundaryEdges)
{
FatalErrorInFunction
<< "Sizing mismatch. Expected " << nBoundaryEdges
<< " boundary edge connections, but had "
<< bndEdgeConnections.size() << nl
<< abort(FatalError);
}
bndConnectPtr_.reset
(
new List<labelPair>(nBoundaryEdges, labelPair(-1,-1))
);
auto& bndConnect = *bndConnectPtr_;
for (const auto& connection : bndEdgeConnections)
{
const auto& a = connection.first();
const auto& b = connection.second();
if (a.is_finiteArea() && a.is_localProc())
{
const label bndEdgei = (a.patchEdgei() - nInternalEdges);
bndConnect[bndEdgei].first() = b.procNo();
bndConnect[bndEdgei].second() = b.meshFacei();
}
else if (b.is_finiteArea() && b.is_localProc())
{
const label bndEdgei = (b.patchEdgei() - nInternalEdges);
bndConnect[bndEdgei].first() = a.procNo();
bndConnect[bndEdgei].second() = a.meshFacei();
}
else
{
FatalErrorInFunction
<< "Unexpected pairing input " << connection
<< " ... programming error" << nl
<< abort(FatalError);
}
}
label nInvalid = 0;
for (const auto& connection : bndConnect)
{
if (connection.first() < 0 || connection.second() < 0)
{
++nInvalid;
}
}
if (Pstream::parRun())
{
reduce(nInvalid, sumOp<label>());
}
if (nInvalid)
{
FatalErrorInFunction
<< "Did not properly match " << nInvalid
<< " boundary edges" << nl
<< abort(FatalError);
}
}
void Foam::faMesh::calcBoundaryConnections() const
{
setBoundaryConnections(this->getBoundaryEdgeConnections());
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::labelList Foam::faMesh::boundaryProcs() const
{
const auto& connections = this->boundaryConnections();
labelHashSet procsUsed(2*Pstream::nProcs());
for (const labelPair& tuple : connections)
{
procsUsed.insert(tuple.first());
}
procsUsed.erase(-1); // placeholder value
procsUsed.erase(Pstream::myProcNo());
return procsUsed.sortedToc();
}
Foam::List<Foam::labelPair> Foam::faMesh::boundaryProcSizes() const
{
const auto& connections = this->boundaryConnections();
Map<label> procCount(2*Pstream::nProcs());
for (const labelPair& tuple : connections)
{
++procCount(tuple.first());
}
procCount.erase(-1); // placeholder value
procCount.erase(Pstream::myProcNo());
// Flatten as list
List<labelPair> output(procCount.size());
label count = 0;
for (const label proci : procCount.sortedToc())
{
output[count].first() = proci;
output[count].second() = procCount[proci]; // size
++count;
}
return output;
}
const Foam::faMeshBoundaryHalo& Foam::faMesh::boundaryHaloMap() const
{
if (!haloMapPtr_)
{
haloMapPtr_.reset(new faMeshBoundaryHalo(*this));
}
return *haloMapPtr_;
}
void Foam::faMesh::calcHaloFaceGeometry() const
{
if (haloFaceCentresPtr_ || haloFaceNormalsPtr_)
{
FatalErrorInFunction
<< "Halo centres/normals already calculated"
<< exit(FatalError);
}
DebugInFunction
<< "Calculating halo face centres/normals" << endl;
const faceList& faces = mesh().faces();
const pointField& points = mesh().points();
const faMeshBoundaryHalo& halo = boundaryHaloMap();
const labelList& inputFaceIds = halo.inputMeshFaces();
haloFaceCentresPtr_.reset(new pointField);
haloFaceNormalsPtr_.reset(new vectorField);
auto& centres = *haloFaceCentresPtr_;
auto& normals = *haloFaceNormalsPtr_;
centres.resize(inputFaceIds.size());
normals.resize(inputFaceIds.size());
// My values
forAll(inputFaceIds, i)
{
const face& f = faces[inputFaceIds[i]];
centres[i] = f.centre(points);
normals[i] = f.unitNormal(points);
}
// Swap information and resize
halo.distributeSparse(centres);
halo.distributeSparse(normals);
}
const Foam::pointField& Foam::faMesh::haloFaceCentres() const
{
if (!haloFaceCentresPtr_ || !haloFaceNormalsPtr_)
{
calcHaloFaceGeometry();
}
return *haloFaceCentresPtr_;
}
const Foam::vectorField& Foam::faMesh::haloFaceNormals() const
{
if (!haloFaceCentresPtr_ || !haloFaceNormalsPtr_)
{
calcHaloFaceGeometry();
}
return *haloFaceNormalsPtr_;
}
Foam::tmp<Foam::pointField>
Foam::faMesh::haloFaceCentres(const label patchi) const
{
if (patchi < 0 || patchi >= boundary().size())
{
FatalErrorInFunction
<< "Patch " << patchi << " is out-of-range 0.."
<< (boundary().size()-1) << nl
<< exit(FatalError);
}
return tmp<pointField>::New
(
List<point>
(
boundarySubset
(
this->haloFaceCentres(),
boundary()[patchi].edgeLabels()
)
)
);
}
Foam::tmp<Foam::vectorField>
Foam::faMesh::haloFaceNormals(const label patchi) const
{
if (patchi < 0 || patchi >= boundary().size())
{
FatalErrorInFunction
<< "Patch " << patchi << " is out-of-range 0.."
<< (boundary().size()-1) << nl
<< exit(FatalError);
}
return tmp<vectorField>::New
(
List<vector>
(
boundarySubset
(
this->haloFaceNormals(),
boundary()[patchi].edgeLabels()
)
)
);
}
// ************************************************************************* //

37
src/finiteArea/faMesh/faPatches/constraint/cyclic/cyclicFaPatch.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -47,34 +47,33 @@ const Foam::scalar Foam::cyclicFaPatch::matchTol_ = 1e-3;
void Foam::cyclicFaPatch::calcTransforms()
{
const label sizeby2 = this->size()/2;
if (size() > 0)
{
// const label sizeby2 = this->size()/2;
pointField half0Ctrs(size()/2);
pointField half1Ctrs(size()/2);
for (label i=0; i<size()/2; ++i)
pointField half0Ctrs(sizeby2);
pointField half1Ctrs(sizeby2);
for (label edgei=0; edgei < sizeby2; ++edgei)
{
half0Ctrs[i] = this->edgeCentres()[i];
half1Ctrs[i] = this->edgeCentres()[i+size()/2];
half0Ctrs[edgei] = this->edgeCentres()[edgei];
half1Ctrs[edgei] = this->edgeCentres()[edgei + sizeby2];
}
vectorField half0Normals(size()/2);
vectorField half1Normals(size()/2);
vectorField half0Normals(sizeby2);
vectorField half1Normals(sizeby2);
const vectorField eN(edgeNormals()*magEdgeLengths());
scalar maxMatchError = 0;
label errorEdge = -1;
for (label edgei = 0; edgei < size()/2; ++edgei)
for (label edgei = 0; edgei < sizeby2; ++edgei)
{
half0Normals[edgei] = eN[edgei];
label nbrEdgei = edgei + size()/2;
half1Normals[edgei] = eN[nbrEdgei];
half1Normals[edgei] = eN[edgei + sizeby2];
scalar magLe = mag(half0Normals[edgei]);
scalar nbrMagLe = mag(half1Normals[edgei]);
scalar avLe = (magLe + nbrMagLe)/2.0;
scalar avLe = 0.5*(magLe + nbrMagLe);
if (magLe < ROOTVSMALL && nbrMagLe < ROOTVSMALL)
{
@ -101,10 +100,10 @@ void Foam::cyclicFaPatch::calcTransforms()
// Check for error in edge matching
if (maxMatchError > matchTol_)
{
label nbrEdgei = errorEdge + size()/2;
label nbrEdgei = errorEdge + sizeby2;
scalar magLe = mag(half0Normals[errorEdge]);
scalar nbrMagLe = mag(half1Normals[errorEdge]);
scalar avLe = (magLe + nbrMagLe)/2.0;
scalar avLe = 0.5*(magLe + nbrMagLe);
FatalErrorInFunction
<< "edge " << errorEdge
@ -162,7 +161,7 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
for (label edgei = 0; edgei < sizeby2; ++edgei)
{
scalar avL = (magL[edgei] + magL[edgei + sizeby2])/2.0;
scalar avL = 0.5*(magL[edgei] + magL[edgei + sizeby2]);
if
(
@ -191,7 +190,7 @@ void Foam::cyclicFaPatch::makeWeights(scalarField& w) const
// Check for error in matching
if (maxMatchError > matchTol_)
{
scalar avL = (magL[errorEdge] + magL[errorEdge + sizeby2])/2.0;
scalar avL = 0.5*(magL[errorEdge] + magL[errorEdge + sizeby2]);
FatalErrorInFunction
<< "edge " << errorEdge << " and " << errorEdge + sizeby2
@ -313,7 +312,7 @@ Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::transfer
const label sizeby2 = this->size()/2;
for (label edgei=0; edgei<sizeby2; ++edgei)
for (label edgei=0; edgei < sizeby2; ++edgei)
{
pnf[edgei] = interfaceData[edgei + sizeby2];
pnf[edgei + sizeby2] = interfaceData[edgei];
@ -345,7 +344,7 @@ Foam::tmp<Foam::labelField> Foam::cyclicFaPatch::internalFieldTransfer
const label sizeby2 = this->size()/2;
for (label edgei=0; edgei<sizeby2; ++edgei)
for (label edgei=0; edgei < sizeby2; ++edgei)
{
pnf[edgei] = iF[edgeCells[edgei + sizeby2]];
pnf[edgei + sizeby2] = iF[edgeCells[edgei]];

7
src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.C
View File

@ -80,12 +80,7 @@ void Foam::processorFaPatch::makeNonGlobalPatchPoints() const
|| !boundaryMesh().mesh()().globalData().nGlobalPoints()
)
{
nonGlobalPatchPointsPtr_ = new labelList(nPoints());
labelList& ngpp = *nonGlobalPatchPointsPtr_;
forAll(ngpp, i)
{
ngpp[i] = i;
}
nonGlobalPatchPointsPtr_ = new labelList(identity(nPoints()));
}
else
{

19
src/finiteArea/faMesh/faPatches/constraint/processor/processorFaPatch.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2016-2017 Wikki Ltd
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -177,19 +177,19 @@ public:
// Member Functions
//- Return interface size
virtual label interfaceSize() const
virtual label interfaceSize() const noexcept
{
return size();
}
//- Return processor number
int myProcNo() const
int myProcNo() const noexcept
{
return myProcNo_;
}
//- Return neighbour processor number
int neighbProcNo() const
int neighbProcNo() const noexcept
{
return neighbProcNo_;
}
@ -197,18 +197,11 @@ public:
//- Return true if running parallel
virtual bool coupled() const
{
if (Pstream::parRun())
{
return true;
}
else
{
return false;
}
return Pstream::parRun();
}
//- Is this the master side?
virtual bool master() const
virtual bool master() const noexcept
{
return (myProcNo_ < neighbProcNo_);
}

332
src/finiteArea/faMesh/faPatches/faPatch/faPatch.C
View File

@ -32,6 +32,7 @@ License
#include "faMesh.H"
#include "areaFields.H"
#include "edgeFields.H"
#include "edgeHashes.H"
#include "polyMesh.H"
#include "demandDrivenData.H"
@ -116,12 +117,6 @@ Foam::faPatch::~faPatch()
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::faPatch::ngbPolyPatchIndex() const noexcept
{
return nbrPolyPatchId_;
}
const Foam::faBoundaryMesh& Foam::faPatch::boundaryMesh() const noexcept
{
return boundaryMesh_;
@ -134,6 +129,76 @@ Foam::label Foam::faPatch::start() const
}
Foam::List<Foam::labelPair> Foam::faPatch::boundaryConnections() const
{
const auto& connections = boundaryMesh().mesh().boundaryConnections();
const label nInternalEdges = boundaryMesh().mesh().nInternalEdges();
List<labelPair> output(this->nEdges());
// Like an IndirectList but removing the nInternalEdges offset
label count = 0;
for (const label patchEdgei : this->edgeLabels())
{
const label bndEdgei = (patchEdgei - nInternalEdges);
output[count] = connections[bndEdgei];
++count;
}
return output;
}
Foam::labelList Foam::faPatch::boundaryProcs() const
{
const auto& connections = boundaryMesh().mesh().boundaryConnections();
const label nInternalEdges = boundaryMesh().mesh().nInternalEdges();
labelHashSet procsUsed(2*Pstream::nProcs());
for (const label patchEdgei : this->edgeLabels())
{
const label bndEdgei = (patchEdgei - nInternalEdges);
const label proci = connections[bndEdgei].first();
procsUsed.insert(proci);
}
procsUsed.erase(-1); // placeholder value
procsUsed.erase(Pstream::myProcNo());
return procsUsed.sortedToc();
}
Foam::List<Foam::labelPair> Foam::faPatch::boundaryProcSizes() const
{
const auto& connections = boundaryMesh().mesh().boundaryConnections();
const label nInternalEdges = boundaryMesh().mesh().nInternalEdges();
Map<label> procCount(2*Pstream::nProcs());
for (const label patchEdgei : this->edgeLabels())
{
const label bndEdgei = (patchEdgei - nInternalEdges);
const label proci = connections[bndEdgei].first();
++procCount(proci);
}
procCount.erase(-1); // placeholder value
procCount.erase(Pstream::myProcNo());
// Flatten as list
List<labelPair> output(procCount.size());
label count = 0;
for (const label proci : procCount.sortedToc())
{
output[count].first() = proci;
output[count].second() = procCount[proci]; // size
++count;
}
return output;
}
const Foam::labelList& Foam::faPatch::pointLabels() const
{
if (!pointLabelsPtr_)
@ -145,88 +210,6 @@ const Foam::labelList& Foam::faPatch::pointLabels() const
}
void Foam::faPatch::calcPointLabels() const
{
SLList<label> labels;
UList<edge> edges = patchSlice(boundaryMesh().mesh().edges());
forAll(edges, edgeI)
{
bool existStart = false;
bool existEnd = false;
forAllIters(labels, iter)
{
if (*iter == edges[edgeI].start())
{
existStart = true;
}
if (*iter == edges[edgeI].end())
{
existEnd = true;
}
}
if (!existStart)
{
labels.append(edges[edgeI].start());
}
if (!existEnd)
{
labels.append(edges[edgeI].end());
}
}
pointLabelsPtr_ = new labelList(labels);
}
void Foam::faPatch::calcPointEdges() const
{
const labelList& points = pointLabels();
const edgeList::subList e = patchSlice(boundaryMesh().mesh().edges());
// set up storage for pointEdges
List<SLList<label>> pointEdgs(points.size());
forAll(e, edgeI)
{
const edge& curPoints = e[edgeI];
forAll(curPoints, pointI)
{
const label localPointIndex = points.find(curPoints[pointI]);
pointEdgs[localPointIndex].append(edgeI);
}
}
// sort out the list
pointEdgesPtr_ = new labelListList(pointEdgs.size());
labelListList& pEdges = *pointEdgesPtr_;
forAll(pointEdgs, pointI)
{
pEdges[pointI].setSize(pointEdgs[pointI].size());
label i = 0;
for
(
SLList<label>::iterator curEdgesIter = pointEdgs[pointI].begin();
curEdgesIter != pointEdgs[pointI].end();
++curEdgesIter, ++i
)
{
pEdges[pointI][i] = curEdgesIter();
}
}
}
const Foam::labelListList& Foam::faPatch::pointEdges() const
{
if (!pointEdgesPtr_)
@ -238,55 +221,91 @@ const Foam::labelListList& Foam::faPatch::pointEdges() const
}
Foam::labelList Foam::faPatch::ngbPolyPatchFaces() const
void Foam::faPatch::calcPointLabels() const
{
if (nbrPolyPatchId_ < 0)
const edgeList::subList edges = patchSlice(boundaryMesh().mesh().edges());
// Walk boundary edges.
// The edge orientation corresponds to the face orientation
// (outwards normal).
// Note: could combine this with calcPointEdges for more efficiency
// Map<label> markedPoints(4*edges.size());
labelHashSet markedPoints(4*edges.size());
DynamicList<label> dynEdgePoints(2*edges.size());
for (const edge& e : edges)
{
return labelList();
}
labelList ngbFaces(faPatch::size());
const faMesh& aMesh = boundaryMesh().mesh();
const polyMesh& pMesh = aMesh.mesh();
const auto& patch = aMesh.patch();
const labelListList& edgeFaces = pMesh.edgeFaces();
const labelList meshEdges
(
patch.meshEdges(pMesh.edges(), pMesh.pointEdges())
);
forAll(ngbFaces, edgeI)
{
ngbFaces[edgeI] = -1;
label curEdge = (*this)[edgeI];
label curPMeshEdge = meshEdges[curEdge];
forAll(edgeFaces[curPMeshEdge], faceI)
// if (markedPoints.insert(e.first(), markedPoints.size()))
if (markedPoints.insert(e.first()))
{
label curFace = edgeFaces[curPMeshEdge][faceI];
label curPatchID = pMesh.boundaryMesh().whichPatch(curFace);
if (curPatchID == nbrPolyPatchId_)
{
ngbFaces[edgeI] = curFace;
}
dynEdgePoints.append(e.first());
}
if (ngbFaces[edgeI] == -1)
// if (markedPoints.insert(e.second(), markedPoints.size()))
if (markedPoints.insert(e.second()))
{
WarningInFunction
<< "Problem with determination of edge ngb faces!"
<< endl;
dynEdgePoints.append(e.second());
}
}
return ngbFaces;
// Transfer to plain list (reuse storage)
pointLabelsPtr_ = new labelList(std::move(dynEdgePoints));
/// const labelList& edgePoints = *pointLabelsPtr_;
///
/// // Cannot use invertManyToMany - we have non-local edge numbering
///
/// // Intermediate storage for pointEdges.
/// // Points on the boundary will normally connect 1 or 2 edges only.
/// List<DynamicList<label,2>> dynPointEdges(edgePoints.size());
///
/// forAll(edges, edgei)
/// {
/// const edge& e = edges[edgei];
///
/// dynPointEdges[markedPoints[e.first()]].append(edgei);
/// dynPointEdges[markedPoints[e.second()]].append(edgei);
/// }
///
/// // Flatten to regular list
/// pointEdgesPtr_ = new labelListList(edgePoints.size());
/// auto& pEdges = *pointEdgesPtr_;
///
/// forAll(pEdges, pointi)
/// {
/// pEdges[pointi] = std::move(dynPointEdges[pointi]);
/// }
}
void Foam::faPatch::calcPointEdges() const
{
const edgeList::subList edges = patchSlice(boundaryMesh().mesh().edges());
const labelList& edgePoints = pointLabels();
// Cannot use invertManyToMany - we have non-local edge numbering
// Intermediate storage for pointEdges.
// Points on the boundary will normally connect 1 or 2 edges only.
List<DynamicList<label,2>> dynPointEdges(edgePoints.size());
forAll(edges, edgei)
{
const edge& e = edges[edgei];
dynPointEdges[edgePoints.find(e.first())].append(edgei);
dynPointEdges[edgePoints.find(e.second())].append(edgei);
}
// Flatten to regular list
pointEdgesPtr_ = new labelListList(edgePoints.size());
auto& pEdges = *pointEdgesPtr_;
forAll(pEdges, pointi)
{
pEdges[pointi] = std::move(dynPointEdges[pointi]);
}
}
@ -297,24 +316,7 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchFaceNormals() const
return tmp<vectorField>::New();
}
auto tfN = tmp<vectorField>::New();
auto& fN = tfN.ref();
fN.setSize(faPatch::size());
labelList ngbFaces = ngbPolyPatchFaces();
const polyMesh& pMesh = boundaryMesh().mesh()();
const faceList& faces = pMesh.faces();
const pointField& points = pMesh.points();
forAll(fN, faceI)
{
fN[faceI] = faces[ngbFaces[faceI]].unitNormal(points);
}
return tfN;
return boundaryMesh().mesh().haloFaceNormals(this->index());
}
@ -325,24 +327,31 @@ Foam::tmp<Foam::vectorField> Foam::faPatch::ngbPolyPatchPointNormals() const
return tmp<vectorField>::New();
}
// Unit normals for the neighbour patch faces
const vectorField faceNormals
(
boundaryMesh().mesh().haloFaceNormals(this->index())
);
const labelListList& pntEdges = pointEdges();
auto tpN = tmp<vectorField>::New(pntEdges.size(), Zero);
auto& pN = tpN.ref();
auto tpointNorm = tmp<vectorField>::New(pntEdges.size());
auto& pointNorm = tpointNorm.ref();
const vectorField faceNormals(ngbPolyPatchFaceNormals());
forAll(pN, pointI)
forAll(pointNorm, pointi)
{
forAll(pntEdges[pointI], edgeI)
vector& n = pointNorm[pointi];
n = Zero;
for (const label bndEdgei : pntEdges[pointi])
{
pN[pointI] += faceNormals[pntEdges[pointI][edgeI]];
n += faceNormals[bndEdgei];
}
n.normalise();
}
pN /= mag(pN);
return tpN;
return tpointNorm;
}
@ -380,11 +389,14 @@ const Foam::scalarField& Foam::faPatch::magEdgeLengths() const
Foam::tmp<Foam::vectorField> Foam::faPatch::edgeNormals() const
{
tmp<vectorField> eN(new vectorField(size()));
auto tedgeNorm = tmp<vectorField>::New(edgeLengths());
eN.ref() = edgeLengths()/magEdgeLengths();
for (vector& n : tedgeNorm.ref())
{
n.normalise();
}
return eN;
return tedgeNorm;
}

42
src/finiteArea/faMesh/faPatches/faPatch/faPatch.H
View File

@ -80,17 +80,14 @@ class faPatch
//- Reference to boundary mesh
const faBoundaryMesh& boundaryMesh_;
//- Demand-driven: edge-face addressing
mutable labelList::subList* edgeFacesPtr_;
// Demand-driven private data
//- Demand-driven: local points labels
mutable labelList* pointLabelsPtr_;
//- Edge-face addressing
mutable labelList::subList* edgeFacesPtr_;
//- Local points labels
mutable labelList* pointLabelsPtr_;
//- Point-edge addressing
mutable labelListList* pointEdgesPtr_;
//- Demand-driven: point-edge addressing
mutable labelListList* pointEdgesPtr_;
// Private Member Functions
@ -231,6 +228,7 @@ public:
return static_cast<const labelList&>(*this);
}
//- Define new list of edges
void edgeLabels(const UList<label>& newEdgeLabels);
//- Number of patch points
@ -245,8 +243,11 @@ public:
return labelList::size();
}
//- Return neighbour polyPatch index
label ngbPolyPatchIndex() const noexcept;
//- The neighbour polyPatch index
label ngbPolyPatchIndex() const noexcept
{
return nbrPolyPatchId_;
}
//- Return boundaryMesh reference
const faBoundaryMesh& boundaryMesh() const noexcept;
@ -284,6 +285,21 @@ public:
virtual void write(Ostream&) const;
// Topological information
//- List of proc/face for the boundary edge neighbours
//- in locally reordered edge numbering.
List<labelPair> boundaryConnections() const;
//- Boundary edge neighbour processors
//- (does not include own proc)
labelList boundaryProcs() const;
//- List of proc/size for the boundary edge neighbour processors
//- (does not include own proc)
List<labelPair> boundaryProcSizes() const;
// Access functions for geometrical data
//- Return patch point labels
@ -292,10 +308,8 @@ public:
//- Return patch point-edge addressing
const labelListList& pointEdges() const;
//- Return edge neighbour polyPatch faces
labelList ngbPolyPatchFaces() const;
//- Return normals of neighbour polyPatch faces
// Same as faMesh::haloFaceNormals()
tmp<vectorField> ngbPolyPatchFaceNormals() const;
//- Return normals of neighbour polyPatch joined points

19
src/finiteArea/faMesh/faPatches/faPatch/faPatchData.C
View File

@ -28,6 +28,7 @@ License
#include "faPatchData.H"
#include "dictionary.H"
#include "processorFaPatch.H"
#include "processorPolyPatch.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -108,6 +109,24 @@ void Foam::faPatchData::assign(const faPatch& fap)
}
bool Foam::faPatchData::assign_coupled(int ownProci, int neiProci)
{
clear();
if (ownProci == neiProci)
{
return false;
}
name_ = processorPolyPatch::newName(ownProci, neiProci);
type_ = processorFaPatch::typeName;
ownerProcId_ = ownProci;
neighProcId_ = neiProci;
return true;
}
int Foam::faPatchData::matchPatchPair
(
const labelPair& patchPair

3
src/finiteArea/faMesh/faPatches/faPatch/faPatchData.H
View File

@ -98,6 +98,9 @@ public:
//- Clear and populate with values from finiteArea patch
void assign(const faPatch& fap);
//- Set values consistent with a processor coupling
bool assign_coupled(int ownProci, int neiProci);
//- True if owner/neighbour processor ids are non-equal
bool coupled() const noexcept
{

15
src/parallel/decompose/faDecompose/faFieldDecomposer.C
View File

@ -293,7 +293,7 @@ void Foam::faFieldDecomposer::reset
new processorEdgePatchFieldDecomposer
(
procMesh_.boundary()[patchi].size(),
static_cast<const labelUList&>(localPatchSlice)
localPatchSlice
)
);
}
@ -309,6 +309,15 @@ void Foam::faFieldDecomposer::reset(const faMesh& completeMesh)
processorAreaPatchFieldDecomposerPtrs_.resize(nMappers);
processorEdgePatchFieldDecomposerPtrs_.resize(nMappers);
// Create weightings now - needed for proper parallel synchronization
(void)completeMesh.weights();
// faPatches don't have their own start() - so these are invariant
const labelList completePatchStarts
(
completeMesh.boundary().patchStarts()
);
forAll(boundaryAddressing_, patchi)
{
const label oldPatchi = boundaryAddressing_[patchi];
@ -324,7 +333,7 @@ void Foam::faFieldDecomposer::reset(const faMesh& completeMesh)
(
completeMesh.boundary()[oldPatchi].size(),
localPatchSlice,
completeMesh.boundary()[oldPatchi].start()
completePatchStarts[oldPatchi]
)
);
}
@ -346,7 +355,7 @@ void Foam::faFieldDecomposer::reset(const faMesh& completeMesh)
new processorEdgePatchFieldDecomposer
(
procMesh_.boundary()[patchi].size(),
static_cast<const labelUList&>(localPatchSlice)
localPatchSlice
)
);
}

4
src/parallel/decompose/faDecompose/faFieldDecomposerFields.C
View File

@ -219,9 +219,9 @@ void Foam::faFieldDecomposer::decomposeFields
const PtrList<GeoField>& fields
) const
{
forAll(fields, fieldI)
forAll(fields, fieldi)
{
decomposeField(fields[fieldI])().write();
decomposeField(fields[fieldi])().write();
}
}

98
src/parallel/reconstruct/faReconstruct/faMeshReconstructor.C
View File

@ -88,11 +88,85 @@ void Foam::faMeshReconstructor::calcAddressing
{
globalFaceNum.gather(faFaceProcAddr_, singlePatchFaceLabels_);
labelList order(Foam::sortedOrder(singlePatchFaceLabels_));
const labelList globalOrder(Foam::sortedOrder(singlePatchFaceLabels_));
singlePatchFaceLabels_ = labelList(singlePatchFaceLabels_, order);
singlePatchFaceLabels_ =
labelList(singlePatchFaceLabels_, globalOrder);
globalFaceNum.scatter(order, faFaceProcAddr_);
// Set first face to be zero relative to the finiteArea patch
// ie, local-face numbering with the first being 0 on any given patch
{
label patchFirstMeshfacei
(
singlePatchFaceLabels_.empty()
? 0
: singlePatchFaceLabels_.first()
);
Pstream::scatter(patchFirstMeshfacei);
for (label& facei : faFaceProcAddr_)
{
facei -= patchFirstMeshfacei;
}
}
PstreamBuffers pBufs(Pstream::commsTypes::nonBlocking);
if (Pstream::master())
{
// Determine the respective local portions of the global ordering
labelList procTargets(globalFaceNum.size());
for (const label proci : Pstream::allProcs())
{
labelList::subList
(
globalFaceNum.range(proci),
procTargets
) = proci;
}
labelList procStarts(globalFaceNum.offsets());
labelList procOrders(globalFaceNum.size());
for (const label globali : globalOrder)
{
const label proci = procTargets[globali];
procOrders[procStarts[proci]++] =
(globali - globalFaceNum.localStart(proci));
}
// Send the local portions
for (const int proci : Pstream::subProcs())
{
SubList<label> localOrder
(
procOrders,
globalFaceNum.range(proci)
);
UOPstream toProc(proci, pBufs);
toProc << localOrder;
}
SubList<label> localOrder(procOrders, globalFaceNum.range(0));
faFaceProcAddr_ = labelList(faFaceProcAddr_, localOrder);
}
pBufs.finishedSends();
if (!Pstream::master())
{
labelList localOrder;
UIPstream fromProc(Pstream::master(), pBufs);
fromProc >> localOrder;
faFaceProcAddr_ = labelList(faFaceProcAddr_, localOrder);
}
}
// Broadcast the same information everywhere
@ -201,12 +275,14 @@ void Foam::faMeshReconstructor::calcAddressing
tmpFaces = initialPatch.localFaces();
pointField tmpPoints(initialPatch.localPoints());
// The meshPointMap is contiguous, so flatten as linear list
/// Map<label> mpm(initialPatch.meshPointMap());
labelList mpm(initialPatch.nPoints());
forAllConstIters(initialPatch.meshPointMap(), iter)
// Equivalent to a flattened meshPointMap
labelList mpm(initialPatch.points().size(), -1);
{
mpm[iter.key()] = iter.val();
const labelList& mp = initialPatch.meshPoints();
forAll(mp, i)
{
mpm[mp[i]] = i;
}
}
Pstream::scatter(mpm);
@ -441,7 +517,13 @@ void Foam::faMeshReconstructor::createMesh()
);
}
// Serial mesh - no parallel communication
const bool oldParRun = Pstream::parRun(false);
completeMesh.addFaPatches(completePatches);
Pstream::parRun(oldParRun); // Restore parallel state
}