This repository has been archived on 2023-11-20. You can view files and clone it, but cannot push or open issues or pull requests.
foam-extend4.1-coherent-io/applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C
2010-12-02 14:35:12 +00:00

764 lines
22 KiB
C

/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright held by original author
\\/ M anipulation |
-------------------------------------------------------------------------------
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 2 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, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Application
reconstructParMesh
Description
Reconstructs a mesh using geometric information only.
Writes point/face/cell procAddressing so afterwards reconstructPar can be
used to reconstruct fields.
Note:
- uses geometric matching tolerance (set with -mergeTol option)
If the parallel case does not have correct procBoundaries use the
-fullMatch option which will check all boundary faces (bit slower).
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "Time.H"
#include "IOobjectList.H"
#include "labelIOList.H"
#include "processorPolyPatch.H"
#include "mapAddedPolyMesh.H"
#include "polyMeshAdder.H"
#include "faceCoupleInfo.H"
#include "fvMeshAdder.H"
#include "directTopoChange.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Tolerance (as fraction of the bounding box). Needs to be fairly lax since
// usually meshes get written with limited precision (6 digits)
static const scalar defaultMergeTol = 1e-7;
static void renumber
(
const labelList& map,
labelList& elems
)
{
forAll(elems, i)
{
if (elems[i] >= 0)
{
elems[i] = map[elems[i]];
}
}
}
// Determine which faces are coupled. Uses geometric merge distance.
// Looks either at all boundaryFaces (fullMatch) or only at the
// procBoundaries for procI. Assumes that masterMesh contains already merged
// all the processors < procI.
autoPtr<faceCoupleInfo> determineCoupledFaces
(
const bool fullMatch,
const label procI,
const polyMesh& masterMesh,
const polyMesh& meshToAdd,
const scalar mergeDist
)
{
if (fullMatch || masterMesh.nCells() == 0)
{
return autoPtr<faceCoupleInfo>
(
new faceCoupleInfo
(
masterMesh,
meshToAdd,
mergeDist, // absolute merging distance
true // matching faces identical
)
);
}
else
{
// Pick up all patches on masterMesh ending in "toDDD" where DDD is
// the processor number procI.
const polyBoundaryMesh& masterPatches = masterMesh.boundaryMesh();
const string toProcString("to" + name(procI));
DynamicList<label> masterFaces
(
masterMesh.nFaces()
- masterMesh.nInternalFaces()
);
forAll(masterPatches, patchI)
{
const polyPatch& pp = masterPatches[patchI];
if
(
isA<processorPolyPatch>(pp)
&& (
pp.name().rfind(toProcString)
== (pp.name().size()-toProcString.size())
)
)
{
label meshFaceI = pp.start();
forAll(pp, i)
{
masterFaces.append(meshFaceI++);
}
}
}
masterFaces.shrink();
// Pick up all patches on meshToAdd ending in "procBoundaryDDDtoYYY"
// where DDD is the processor number procI and YYY is < procI.
const polyBoundaryMesh& addPatches = meshToAdd.boundaryMesh();
DynamicList<label> addFaces
(
meshToAdd.nFaces()
- meshToAdd.nInternalFaces()
);
forAll(addPatches, patchI)
{
const polyPatch& pp = addPatches[patchI];
if (isA<processorPolyPatch>(pp))
{
bool isConnected = false;
for (label mergedProcI = 0; mergedProcI < procI; mergedProcI++)
{
const string fromProcString
(
"procBoundary"
+ name(procI)
+ "to"
+ name(mergedProcI)
);
if (pp.name() == fromProcString)
{
isConnected = true;
break;
}
}
if (isConnected)
{
label meshFaceI = pp.start();
forAll(pp, i)
{
addFaces.append(meshFaceI++);
}
}
}
}
addFaces.shrink();
return autoPtr<faceCoupleInfo>
(
new faceCoupleInfo
(
masterMesh,
masterFaces,
meshToAdd,
addFaces,
mergeDist, // absolute merging distance
true, // matching faces identical?
false, // if perfectmatch are faces already ordered
// (e.g. processor patches)
false // are faces each on separate patch?
)
);
}
}
autoPtr<mapPolyMesh> mergeSharedPoints
(
const scalar mergeDist,
polyMesh& mesh,
labelListList& pointProcAddressing
)
{
// Find out which sets of points get merged and create a map from
// mesh point to unique point.
Map<label> pointToMaster
(
fvMeshAdder::findSharedPoints
(
mesh,
mergeDist
)
);
Info<< "mergeSharedPoints : detected " << pointToMaster.size()
<< " points that are to be merged." << endl;
if (returnReduce(pointToMaster.size(), sumOp<label>()) == 0)
{
return autoPtr<mapPolyMesh>(NULL);
}
directTopoChange meshMod(mesh);
fvMeshAdder::mergePoints(mesh, pointToMaster, meshMod);
// Change the mesh (no inflation). Note: parallel comms allowed.
autoPtr<mapPolyMesh> map = meshMod.changeMesh(mesh, false, true);
// Update fields. No inflation, parallel sync.
mesh.updateMesh(map);
// pointProcAddressing give indices into the master mesh so adapt them
// for changed point numbering.
// Adapt constructMaps for merged points.
forAll(pointProcAddressing, procI)
{
labelList& constructMap = pointProcAddressing[procI];
forAll(constructMap, i)
{
label oldPointI = constructMap[i];
// New label of point after changeMesh.
label newPointI = map().reversePointMap()[oldPointI];
if (newPointI < -1)
{
constructMap[i] = -newPointI-2;
}
else if (newPointI >= 0)
{
constructMap[i] = newPointI;
}
else
{
FatalErrorIn("fvMeshDistribute::mergeSharedPoints()")
<< "Problem. oldPointI:" << oldPointI
<< " newPointI:" << newPointI << abort(FatalError);
}
}
}
return map;
}
int main(int argc, char *argv[])
{
argList::noParallel();
argList::validOptions.insert("mergeTol", "relative merge distance");
argList::validOptions.insert("fullMatch", "");
# include "addTimeOptions.H"
# include "addRegionOption.H"
# include "setRootCase.H"
# include "createTime.H"
Info<< "This is an experimental tool which tries to merge"
<< " individual processor" << nl
<< "meshes back into one master mesh. Use it if the original"
<< " master mesh has" << nl
<< "been deleted or if the processor meshes have been modified"
<< " (topology change)." << nl
<< "This tool will write the resulting mesh to a new time step"
<< " and construct" << nl
<< "xxxxProcAddressing files in the processor meshes so"
<< " reconstructPar can be" << nl
<< "used to regenerate the fields on the master mesh."
<< endl;
word regionName = polyMesh::defaultRegion;
fileName regionPrefix = "";
if (args.optionFound("region"))
{
regionName = args.option("region");
regionPrefix = regionName;
Info<< "Operating on region " << regionName << nl << endl;
}
scalar mergeTol = defaultMergeTol;
args.optionReadIfPresent("mergeTol", mergeTol);
scalar writeTol = Foam::pow(10.0, -scalar(IOstream::defaultPrecision()));
Info<< "Merge tolerance : " << mergeTol << nl
<< "Write tolerance : " << writeTol << endl;
if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
{
FatalErrorIn(args.executable())
<< "Your current settings specify ASCII writing with "
<< IOstream::defaultPrecision() << " digits precision." << endl
<< "Your merging tolerance (" << mergeTol
<< ") is finer than this." << nl
<< "Please change your writeFormat to binary"
<< " or increase the writePrecision" << endl
<< "or adjust the merge tolerance (-mergeTol)."
<< exit(FatalError);
}
const bool fullMatch = args.optionFound("fullMatch");
if (fullMatch)
{
Info<< "Doing geometric matching on all boundary faces." << nl << endl;
}
else
{
Info<< "Doing geometric matching on correct procBoundaries only."
<< nl << "This assumes a correct decomposition." << endl;
}
int nProcs = 0;
while
(
isDir
(
args.rootPath()
/ args.caseName()
/ fileName(word("processor") + name(nProcs))
)
)
{
nProcs++;
}
Info<< "Found " << nProcs << " processor directories" << nl << endl;
// Read all databases.
PtrList<Time> databases(nProcs);
forAll (databases, procI)
{
Info<< "Reading database "
<< args.caseName()/fileName(word("processor") + name(procI))
<< endl;
databases.set
(
procI,
new Time
(
Time::controlDictName,
args.rootPath(),
args.caseName()/fileName(word("processor") + name(procI))
)
);
Time& procTime = databases[procI];
instantList Times = procTime.times();
// set startTime and endTime depending on -time and -latestTime options
# include "checkTimeOptions.H"
procTime.setTime(Times[startTime], startTime);
if (procI > 0 && databases[procI - 1].value() != procTime.value())
{
FatalErrorIn(args.executable())
<< "Time not equal on processors." << nl
<< "Processor:" << procI - 1
<< " time:" << databases[procI-1].value() << nl
<< "Processor:" << procI
<< " time:" << procTime.value()
<< exit(FatalError);
}
}
// Set master time
Info<< "Setting master time to " << databases[0].timeName() << nl << endl;
runTime.setTime(databases[0]);
// Read point on individual processors to determine merge tolerance
// (otherwise single cell domains might give problems)
boundBox bb = boundBox::invertedBox;
for (label procI = 0; procI < nProcs; procI++)
{
fileName pointsInstance
(
databases[procI].findInstance
(
regionPrefix/polyMesh::meshSubDir,
"points"
)
);
// Look for pointProcAddressing next to points instance
IOobject ppAddrHeader
(
"pointProcAddressing",
pointsInstance,
regionPrefix/polyMesh::meshSubDir,
databases[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
if (ppAddrHeader.headerOk())
{
Info<< "Found pointProcAddressing in " << pointsInstance
<< ". Skipping" << endl;
return 0;
}
Info<< "Reading points from "
<< databases[procI].caseName()
<< " for time = " << databases[procI].timeName()
<< nl << endl;
pointIOField points
(
IOobject
(
"points",
pointsInstance,
regionPrefix/polyMesh::meshSubDir,
databases[procI],
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
)
);
boundBox domainBb(points, false);
bb.min() = min(bb.min(), domainBb.min());
bb.max() = max(bb.max(), domainBb.max());
}
const scalar mergeDist = mergeTol*bb.mag();
Info<< "Overall mesh bounding box : " << bb << nl
<< "Relative tolerance : " << mergeTol << nl
<< "Absolute matching distance : " << mergeDist << nl
<< endl;
// Addressing from processor to reconstructed case
labelListList cellProcAddressing(nProcs);
labelListList faceProcAddressing(nProcs);
labelListList pointProcAddressing(nProcs);
labelListList boundaryProcAddressing(nProcs);
// Internal faces on the final reconstructed mesh
label masterInternalFaces;
// Owner addressing on the final reconstructed mesh
labelList masterOwner;
{
// Construct empty mesh.
Info<< "Constructing empty mesh to add to." << nl << endl;
polyMesh masterMesh
(
IOobject
(
regionName,
runTime.timeName(),
runTime,
IOobject::NO_READ
),
xferCopy(pointField()),
xferCopy(faceList()),
xferCopy(cellList())
);
for (label procI = 0; procI < nProcs; procI++)
{
Info<< "Reading mesh to add from "
<< databases[procI].caseName()
<< " for time = " << databases[procI].timeName()
<< nl << endl;
polyMesh meshToAdd
(
IOobject
(
regionName,
databases[procI].timeName(),
databases[procI]
)
);
// Initialize its addressing
cellProcAddressing[procI] = identity(meshToAdd.nCells());
faceProcAddressing[procI] = identity(meshToAdd.nFaces());
pointProcAddressing[procI] = identity(meshToAdd.nPoints());
boundaryProcAddressing[procI] =
identity(meshToAdd.boundaryMesh().size());
// Find geometrically shared points/faces.
autoPtr<faceCoupleInfo> couples = determineCoupledFaces
(
fullMatch,
procI,
masterMesh,
meshToAdd,
mergeDist
);
// Add elements to mesh
Info<< "Adding to master mesh" << nl << endl;
autoPtr<mapAddedPolyMesh> map = polyMeshAdder::add
(
masterMesh,
meshToAdd,
couples
);
// Update all addressing so xxProcAddressing points to correct item
// in masterMesh.
// Processors that were already in masterMesh
for (label mergedI = 0; mergedI < procI; mergedI++)
{
renumber(map().oldCellMap(), cellProcAddressing[mergedI]);
renumber(map().oldFaceMap(), faceProcAddressing[mergedI]);
renumber(map().oldPointMap(), pointProcAddressing[mergedI]);
// Note: boundary is special since can contain -1.
renumber(map().oldPatchMap(), boundaryProcAddressing[mergedI]);
}
// Added processor
renumber(map().addedCellMap(), cellProcAddressing[procI]);
renumber(map().addedFaceMap(), faceProcAddressing[procI]);
renumber(map().addedPointMap(), pointProcAddressing[procI]);
renumber(map().addedPatchMap(), boundaryProcAddressing[procI]);
Info<< endl;
}
// See if any points on the mastermesh have become connected
// because of connections through processor meshes.
mergeSharedPoints(mergeDist, masterMesh, pointProcAddressing);
// Save some properties on the reconstructed mesh
masterInternalFaces = masterMesh.nInternalFaces();
masterOwner = masterMesh.faceOwner();
Info<< "\nWriting merged mesh to "
<< runTime.path()/runTime.timeName()
<< nl << endl;
if (!masterMesh.write())
{
FatalErrorIn(args.executable())
<< "Failed writing polyMesh."
<< exit(FatalError);
}
}
// Write the addressing
Info<< "Reconstructing the addressing from the processor meshes"
<< " to the newly reconstructed mesh" << nl << endl;
forAll(databases, procI)
{
Info<< "Reading processor " << procI << " mesh from "
<< databases[procI].caseName() << endl;
polyMesh procMesh
(
IOobject
(
regionName,
databases[procI].timeName(),
databases[procI]
)
);
// From processor point to reconstructed mesh point
Info<< "Writing pointProcAddressing to "
<< databases[procI].caseName()
/procMesh.facesInstance()
/polyMesh::meshSubDir
<< endl;
labelIOList
(
IOobject
(
"pointProcAddressing",
procMesh.facesInstance(),
polyMesh::meshSubDir,
procMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
pointProcAddressing[procI]
).write();
// From processor face to reconstructed mesh face
Info<< "Writing faceProcAddressing to "
<< databases[procI].caseName()
/procMesh.facesInstance()
/polyMesh::meshSubDir
<< endl;
labelIOList faceProcAddr
(
IOobject
(
"faceProcAddressing",
procMesh.facesInstance(),
polyMesh::meshSubDir,
procMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
faceProcAddressing[procI]
);
// Now add turning index to faceProcAddressing.
// See reconstrurPar for meaning of turning index.
forAll(faceProcAddr, procFaceI)
{
label masterFaceI = faceProcAddr[procFaceI];
if
(
!procMesh.isInternalFace(procFaceI)
&& masterFaceI < masterInternalFaces
)
{
// proc face is now external but used to be internal face.
// Check if we have owner or neighbour.
label procOwn = procMesh.faceOwner()[procFaceI];
label masterOwn = masterOwner[masterFaceI];
if (cellProcAddressing[procI][procOwn] == masterOwn)
{
// No turning. Offset by 1.
faceProcAddr[procFaceI]++;
}
else
{
// Turned face.
faceProcAddr[procFaceI] =
-1 - faceProcAddr[procFaceI];
}
}
else
{
// No turning. Offset by 1.
faceProcAddr[procFaceI]++;
}
}
faceProcAddr.write();
// From processor cell to reconstructed mesh cell
Info<< "Writing cellProcAddressing to "
<< databases[procI].caseName()
/procMesh.facesInstance()
/polyMesh::meshSubDir
<< endl;
labelIOList
(
IOobject
(
"cellProcAddressing",
procMesh.facesInstance(),
polyMesh::meshSubDir,
procMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
cellProcAddressing[procI]
).write();
// From processor patch to reconstructed mesh patch
Info<< "Writing boundaryProcAddressing to "
<< databases[procI].caseName()
/procMesh.facesInstance()
/polyMesh::meshSubDir
<< endl;
labelIOList
(
IOobject
(
"boundaryProcAddressing",
procMesh.facesInstance(),
polyMesh::meshSubDir,
procMesh,
IOobject::NO_READ,
IOobject::NO_WRITE,
false // do not register
),
boundaryProcAddressing[procI]
).write();
Info<< endl;
}
Info<< "End.\n" << endl;
return 0;
}
// ************************************************************************* //