diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C index 92cf26b0c4..0b6acc5ca0 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/cellSizeControlSurfaces/cellSizeControlSurfaces.C @@ -96,8 +96,9 @@ bool Foam::cellSizeControlSurfaces::evalCellSizeFunctions if (cellSizeFunctions_.size()) { - // Initialise to the last (lowest) priority - label previousPriority = cellSizeFunctions_.last().priority(); + // Maintain priority of current hit. Initialise so it always goes + // through at least once. + label previousPriority = -1; forAll(cellSizeFunctions_, i) { diff --git a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C index bf077cece0..2d850aff09 100644 --- a/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C +++ b/applications/utilities/mesh/generation/cvMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshConformToSurface.C @@ -196,7 +196,7 @@ void Foam::conformalVoronoiMesh::buildSurfaceConformation { if (vit->internalPoint() && !vit->nearBoundary()) { - const Foam::point& pt = topoint(vit->point()); + pointFromPoint pt = topoint(vit->point()); const scalar range = sqr(2.0*targetCellSize(pt)); pointIndexHit pHit; diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C index 6cd54bb461..1f72f1d221 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C +++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMesh.C @@ -238,6 +238,25 @@ int main(int argc, char *argv[]) dict.lookup("constructFrom") ); + // Any merging of small edges + const scalar mergeTol(dict.lookupOrDefault("mergeTol", 1e-4)); + + Info<< "Extruding from " << ExtrudeModeNames[mode] + << " using model " << model().type() << endl; + if (flipNormals) + { + Info<< "Flipping normals before extruding" << endl; + } + if (mergeTol > 0) + { + Info<< "Collapsing edges < " << mergeTol << " of bounding box" << endl; + } + else + { + Info<< "Not collapsing any edges after extrusion" << endl; + } + Info<< endl; + // Generated mesh (one of either) autoPtr meshFromMesh; @@ -715,7 +734,7 @@ int main(int argc, char *argv[]) const boundBox& bb = mesh.bounds(); const vector span = bb.span(); - const scalar mergeDim = 1e-4 * bb.minDim(); + const scalar mergeDim = mergeTol * bb.minDim(); Info<< "Mesh bounding box : " << bb << nl << " with span : " << span << nl @@ -726,6 +745,7 @@ int main(int argc, char *argv[]) // Collapse edges // ~~~~~~~~~~~~~~ + if (mergeDim > 0) { Info<< "Collapsing edges < " << mergeDim << " ..." << nl << endl; diff --git a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict index e075acbc4e..ecbc647797 100644 --- a/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict +++ b/applications/utilities/mesh/generation/extrude/extrudeMesh/extrudeMeshDict @@ -88,4 +88,8 @@ sigmaRadialCoeffs // degree wedges. mergeFaces false; //true; +// Merge small edges. Fraction of bounding box. +mergeTol 0; + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C index 6f3f040a0e..8104f9e7c3 100644 --- a/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C +++ b/applications/utilities/miscellaneous/foamFormatConvert/foamFormatConvert.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -206,6 +206,14 @@ int main(int argc, char *argv[]) runTime ); + // cvMesh vertices + writeMeshObject + ( + "internalDelaunayVertices", + regionPrefix, + runTime + ); + if (runTime.writeFormat() == IOstream::ASCII) { // Only do zones when converting from binary to ascii diff --git a/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C b/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C index 9e68e1f726..136dea3882 100644 --- a/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C +++ b/applications/utilities/parallelProcessing/redistributePar/loadOrCreateMesh.C @@ -51,9 +51,6 @@ Foam::autoPtr Foam::loadOrCreateMesh // Check who has a mesh const bool haveMesh = isDir(io.time().path()/io.instance()/meshSubDir); -Pout<< "meshpath:" << io.time().path()/io.instance()/meshSubDir << endl; -Pout<< "haveMesh:" << haveMesh << endl; - if (!haveMesh) { // Create dummy mesh. Only used on procs that don't have mesh. @@ -293,8 +290,10 @@ Pout<< "haveMesh:" << haveMesh << endl; if (!haveMesh) { // We created a dummy mesh file above. Delete it. - //Pout<< "Removing dummy mesh " << io.objectPath() << endl; - rmDir(io.objectPath()); + const fileName meshFiles = io.time().path()/io.instance()/meshSubDir; + //Pout<< "Removing dummy mesh " << meshFiles << endl; + mesh.removeFiles(); + rmDir(meshFiles); } // Force recreation of globalMeshData. diff --git a/bin/foamPack b/bin/foamPack index 9688cd039e..e3d090e2f2 100755 --- a/bin/foamPack +++ b/bin/foamPack @@ -3,7 +3,7 @@ # ========= | # \\ / F ield | OpenFOAM: The Open Source CFD Toolbox # \\ / O peration | -# \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation +# \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation # \\/ M anipulation | #------------------------------------------------------------------------------ # License @@ -107,15 +107,23 @@ then echo "pack manually" 1>&2 foamPackSource $packDir $packFile else - echo "pack with git-archive" 1>&2 - ( cd $packDir && git archive --format=tar --prefix=$packDir/ HEAD) > $packBase.tar + if [ ! -f $packDir/.build ] + then + echo "Error: $packDir/.build does not exists" 1>&2 + echo " Please update this by running e.g. 'wmakePrintBuild -update' in $packDir" 1>&2 + exit 2 + fi - echo "add in time-stamp and lnInclude directories" 1>&2 - tar cf $packBase.tar2 $packDir/.timeStamp $packDir/.build `find -H $packDir -type d -name lnInclude` - tar Af $packBase.tar $packBase.tar2 - echo "gzip tar file" 1>&2 - gzip -c9 $packBase.tar > $packFile + echo "pack with git-archive" 1>&2 && + ( cd $packDir && git archive --format=tar --prefix=$packDir/ HEAD) > $packBase.tar && + + echo "add in time-stamp and lnInclude directories" 1>&2 && + tar cf $packBase.tar2 $packDir/.timeStamp $packDir/.build `find -H $packDir -type d -name lnInclude` && + tar Af $packBase.tar $packBase.tar2 && + + echo "gzip tar file" 1>&2 && + gzip -c9 $packBase.tar > $packFile && rm -f $packBase.tar $packBase.tar2 2>/dev/null fi diff --git a/src/OpenFOAM/meshes/polyMesh/polyMesh.C b/src/OpenFOAM/meshes/polyMesh/polyMesh.C index e8bf1fc6b2..ac47fece91 100644 --- a/src/OpenFOAM/meshes/polyMesh/polyMesh.C +++ b/src/OpenFOAM/meshes/polyMesh/polyMesh.C @@ -1217,7 +1217,7 @@ const Foam::globalMeshData& Foam::polyMesh::globalData() const // Remove all files and some subdirs (eg, sets) void Foam::polyMesh::removeFiles(const fileName& instanceDir) const { - fileName meshFilesPath = thisDb().path()/instanceDir/meshDir(); + fileName meshFilesPath = thisDb().time().path()/instanceDir/meshDir(); rm(meshFilesPath/"points"); rm(meshFilesPath/"faces"); diff --git a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C index 701de4e23a..36d405d597 100644 --- a/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C +++ b/src/dummyThirdParty/ptscotchDecomp/dummyPtscotchDecomp.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -58,7 +58,7 @@ void Foam::ptscotchDecomp::check(const int retVal, const char* str) {} -Foam::label Foam::ptscotchDecomp::decomposeZeroDomains +Foam::label Foam::ptscotchDecomp::decompose ( const fileName& meshPath, const List& initxadj, @@ -72,6 +72,7 @@ Foam::label Foam::ptscotchDecomp::decomposeZeroDomains ( "label ptscotchDecomp::decompose" "(" + "onst fileName&," "const List&, " "const List&, " "const scalarField&, " @@ -84,8 +85,10 @@ Foam::label Foam::ptscotchDecomp::decomposeZeroDomains Foam::label Foam::ptscotchDecomp::decompose ( const fileName& meshPath, - const List& adjncy, - const List& xadj, + const int adjncySize, + const int adjncy[], + const int xadjSize, + const int xadj[], const scalarField& cWeights, List& finalDecomp ) const @@ -94,9 +97,12 @@ Foam::label Foam::ptscotchDecomp::decompose ( "label ptscotchDecomp::decompose" "(" - "const List&, " - "const List&, " - "const scalarField&, " + "const fileName&," + "const int," + "const int," + "const int," + "const int," + "const scalarField&," "List&" ")" ) << notImplementedMessage << exit(FatalError); diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C index c300cc53a7..abe3bc487c 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.C @@ -439,6 +439,86 @@ void Foam::autoLayerDriver::handleFeatureAngle } +//Foam::tmp Foam::autoLayerDriver::undistortedEdgeLength +//( +// const indirectPrimitivePatch& pp, +// const bool relativeSizes, +// const bool faceSize +//) +//{ +// const fvMesh& mesh = meshRefiner_.mesh(); +// +// tmp tfld(new scalarField()); +// scalarField& fld = tfld(); +// +// if (faceSize) +// { +// fld.setSize(pp.size()); +// } +// else +// { +// fld.setSize(pp.nPoints()); +// } +// +// +// if (relativeSizes) +// { +// const scalar edge0Len = meshRefiner_.meshCutter().level0EdgeLength(); +// const labelList& cellLevel = meshRefiner_.meshCutter().cellLevel(); +// +// if (faceSize) +// { +// forAll(pp, i) +// { +// label faceI = pp.addressing()[i]; +// label ownLevel = cellLevel[mesh.faceOwner()[faceI]]; +// fld[i] = edge0Len/(1<(), +// labelMin // null value +// ); +// +// +// forAll(maxPointLevel, pointI) +// { +// // Find undistorted edge size for this level. +// fld[i] = edge0Len/(1<(mesh).instance(), mesh.time(), // register with runTime - static_cast(mesh).readOpt(), + IOobject::NO_READ, static_cast(mesh).writeOpt() ), // io params from original mesh but new name mesh, // original mesh diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H index 6bbe54a708..cb720dc64f 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriver.H @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -460,7 +460,8 @@ class autoLayerDriver pointVectorField& dispVec, pointScalarField& medialRatio, - pointScalarField& medialDist + pointScalarField& medialDist, + pointVectorField& medialVec ) const; //- Main routine to shrink mesh @@ -481,6 +482,7 @@ class autoLayerDriver const pointVectorField& dispVec, const pointScalarField& medialRatio, const pointScalarField& medialDist, + const pointVectorField& medialVec, List& extrudeStatus, pointField& patchDisp, diff --git a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C index b6a9e338aa..35133f062d 100644 --- a/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C +++ b/src/mesh/autoMesh/autoHexMesh/autoHexMeshDriver/autoLayerDriverShrink.C @@ -689,7 +689,8 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo pointVectorField& dispVec, pointScalarField& medialRatio, - pointScalarField& medialDist + pointScalarField& medialDist, + pointVectorField& medialVec ) const { @@ -929,7 +930,7 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo forAll(pointMedialDist, pointI) { medialDist[pointI] = Foam::sqrt(pointMedialDist[pointI].distSqr()); - //medialVec[pointI] = pointMedialDist[pointI].origin(); + medialVec[pointI] = pointMedialDist[pointI].origin(); } } @@ -966,14 +967,14 @@ void Foam::autoLayerDriver::medialAxisSmoothingInfo << " : normalised direction of nearest displacement" << nl << " " << medialDist.name() << " : distance to medial axis" << nl - //<< " " << medialVec.name() - //<< " : nearest point on medial axis" << nl + << " " << medialVec.name() + << " : nearest point on medial axis" << nl << " " << medialRatio.name() << " : ratio of medial distance to wall distance" << nl << endl; dispVec.write(); medialDist.write(); - //medialVec.write(); + medialVec.write(); medialRatio.write(); } } @@ -996,7 +997,7 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance const pointVectorField& dispVec, const pointScalarField& medialRatio, const pointScalarField& medialDist, - //const pointVectorField& medialVec, + const pointVectorField& medialVec, List& extrudeStatus, pointField& patchDisp, @@ -1066,6 +1067,24 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance << str().name() << endl; } + autoPtr medialVecStr; + label medialVertI = 0; + if (debug) + { + medialVecStr.reset + ( + new OFstream + ( + mesh.time().path() + / "thicknessRatioExcludeMedialVec_" + + meshRefiner_.timeName() + + ".obj" + ) + ); + Info<< "Writing points with too large a extrusion distance to " + << medialVecStr().name() << endl; + } + forAll(meshPoints, patchPointI) { if (extrudeStatus[patchPointI] != NOEXTRUDE) @@ -1082,12 +1101,9 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance vector n = patchDisp[patchPointI] / (mag(patchDisp[patchPointI]) + VSMALL); - //vector mVec = mesh.points()[pointI]-medialVec[pointI]; - //scalar mDist = mag(mVec); - //scalar thicknessRatio = - // (n&mVec) - // *thickness[patchPointI] - // /(mDist+VSMALL); + vector mVec = mesh.points()[pointI]-medialVec[pointI]; + mVec /= mag(mVec)+VSMALL; + thicknessRatio *= (n&mVec); if (thicknessRatio > maxThicknessToMedialRatio) { @@ -1103,8 +1119,9 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance minThickness[patchPointI] +thickness[patchPointI] ) - //<< " since near medial at:" << medialVec[pointI] - //<< " with thicknessRatio:" << thicknessRatio + << " medial direction:" << mVec + << " extrusion direction:" << n + << " with thicknessRatio:" << thicknessRatio << endl; } @@ -1124,6 +1141,16 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance vertI++; str()<< "l " << vertI-1 << ' ' << vertI << nl; } + if (medialVecStr.valid()) + { + const point& pt = mesh.points()[pointI]; + meshTools::writeOBJ(medialVecStr(), pt); + medialVertI++; + meshTools::writeOBJ(medialVecStr(), medialVec[pointI]); + medialVertI++; + medialVecStr()<< "l " << medialVertI-1 + << ' ' << medialVertI << nl; + } } } } @@ -1228,6 +1255,15 @@ void Foam::autoLayerDriver::shrinkMeshMedialDistance ) )() ); + + // Above move will have changed the instance only on the points (which + // is correct). + // However the previous mesh written will be the mesh with layers + // (see autoLayerDriver.C) so we now have to force writing all files + // so we can easily step through time steps. Note that if you + // don't write the mesh with layers this is not necessary. + meshRefiner_.mesh().setInstance(meshRefiner_.timeName()); + meshRefiner_.write ( debug, diff --git a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C index e7c4995a36..81d3851837 100644 --- a/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C +++ b/src/parallel/decompose/ptscotchDecomp/ptscotchDecomp.C @@ -2,7 +2,7 @@ ========= | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | - \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\/ M anipulation | ------------------------------------------------------------------------------- License @@ -103,9 +103,9 @@ License - Note: writes out .dgr files for debugging. Run with e.g. + Note: writeGraph=true : writes out .dgr files for debugging. Run with e.g. - mpirun -np 4 dgpart 2 'processor%r.grf' + mpirun -np 4 dgpart 2 'region0_%r.dgr' - %r gets replaced by current processor rank - decompose into 2 domains @@ -167,192 +167,192 @@ void Foam::ptscotchDecomp::check(const int retVal, const char* str) } -//- Does prevention of 0 cell domains and calls ptscotch. -Foam::label Foam::ptscotchDecomp::decomposeZeroDomains -( - const fileName& meshPath, - const List& initadjncy, - const List& initxadj, - const scalarField& initcWeights, - - List& finalDecomp -) const -{ - globalIndex globalCells(initxadj.size()-1); - - bool hasZeroDomain = false; - for (label procI = 0; procI < Pstream::nProcs(); procI++) - { - if (globalCells.localSize(procI) == 0) - { - hasZeroDomain = true; - break; - } - } - - if (!hasZeroDomain) - { - return decompose - ( - meshPath, - initadjncy, - initxadj, - initcWeights, - finalDecomp - ); - } - - - if (debug) - { - Info<< "ptscotchDecomp : have graphs with locally 0 cells." - << " trickling down." << endl; - } - - // Make sure every domain has at least one cell - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // (scotch does not like zero sized domains) - // Trickle cells from processors that have them up to those that - // don't. - - - // Number of cells to send to the next processor - // (is same as number of cells next processor has to receive) - List nSendCells(Pstream::nProcs(), 0); - - for (label procI = nSendCells.size()-1; procI >=1; procI--) - { - label nLocalCells = globalCells.localSize(procI); - if (nLocalCells-nSendCells[procI] < 1) - { - nSendCells[procI-1] = nSendCells[procI]-nLocalCells+1; - } - } - - // First receive (so increasing the sizes of all arrays) - - Field xadj(initxadj); - Field adjncy(initadjncy); - scalarField cWeights(initcWeights); - - if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) - { - // Receive cells from previous processor - IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1); - - Field prevXadj(fromPrevProc); - Field prevAdjncy(fromPrevProc); - scalarField prevCellWeights(fromPrevProc); - - if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1]) - { - FatalErrorIn("ptscotchDecomp::decompose(..)") - << "Expected from processor " << Pstream::myProcNo()-1 - << " connectivity for " << nSendCells[Pstream::myProcNo()-1] - << " nCells but only received " << prevXadj.size() - << abort(FatalError); - } - - // Insert adjncy - prepend(prevAdjncy, adjncy); - // Adapt offsets and prepend xadj - xadj += prevAdjncy.size(); - prepend(prevXadj, xadj); - // Weights - prepend(prevCellWeights, cWeights); - } - - - // Send to my next processor - - if (nSendCells[Pstream::myProcNo()] > 0) - { - // Send cells to next processor - OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1); - - int nCells = nSendCells[Pstream::myProcNo()]; - int startCell = xadj.size()-1 - nCells; - int startFace = xadj[startCell]; - int nFaces = adjncy.size()-startFace; - - // Send for all cell data: last nCells elements - // Send for all face data: last nFaces elements - toNextProc - << Field::subField(xadj, nCells, startCell)-startFace - << Field::subField(adjncy, nFaces, startFace) - << - ( - cWeights.size() - ? static_cast - ( - scalarField::subField(cWeights, nCells, startCell) - ) - : scalarField(0) - ); - - // Remove data that has been sent - if (cWeights.size()) - { - cWeights.setSize(cWeights.size()-nCells); - } - adjncy.setSize(adjncy.size()-nFaces); - xadj.setSize(xadj.size() - nCells); - } - - - // Do decomposition as normal. Sets finalDecomp. - label result = decompose(meshPath, adjncy, xadj, cWeights, finalDecomp); - - - if (debug) - { - Info<< "ptscotchDecomp : have graphs with locally 0 cells." - << " trickling up." << endl; - } - - - // If we sent cells across make sure we undo it - // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - - // Receive back from next processor if I sent something - if (nSendCells[Pstream::myProcNo()] > 0) - { - IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1); - - List nextFinalDecomp(fromNextProc); - - if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()]) - { - FatalErrorIn("parMetisDecomp::decompose(..)") - << "Expected from processor " << Pstream::myProcNo()+1 - << " decomposition for " << nSendCells[Pstream::myProcNo()] - << " nCells but only received " << nextFinalDecomp.size() - << abort(FatalError); - } - - append(nextFinalDecomp, finalDecomp); - } - - // Send back to previous processor. - if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) - { - OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1); - - int nToPrevious = nSendCells[Pstream::myProcNo()-1]; - - toPrevProc << - SubList - ( - finalDecomp, - nToPrevious, - finalDecomp.size()-nToPrevious - ); - - // Remove locally what has been sent - finalDecomp.setSize(finalDecomp.size()-nToPrevious); - } - return result; -} +////- Does prevention of 0 cell domains and calls ptscotch. +//Foam::label Foam::ptscotchDecomp::decomposeZeroDomains +//( +// const fileName& meshPath, +// const List& initadjncy, +// const List& initxadj, +// const scalarField& initcWeights, +// +// List& finalDecomp +//) const +//{ +// globalIndex globalCells(initxadj.size()-1); +// +// bool hasZeroDomain = false; +// for (label procI = 0; procI < Pstream::nProcs(); procI++) +// { +// if (globalCells.localSize(procI) == 0) +// { +// hasZeroDomain = true; +// break; +// } +// } +// +// if (!hasZeroDomain) +// { +// return decompose +// ( +// meshPath, +// initadjncy, +// initxadj, +// initcWeights, +// finalDecomp +// ); +// } +// +// +// if (debug) +// { +// Info<< "ptscotchDecomp : have graphs with locally 0 cells." +// << " trickling down." << endl; +// } +// +// // Make sure every domain has at least one cell +// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// // (scotch does not like zero sized domains) +// // Trickle cells from processors that have them up to those that +// // don't. +// +// +// // Number of cells to send to the next processor +// // (is same as number of cells next processor has to receive) +// List nSendCells(Pstream::nProcs(), 0); +// +// for (label procI = nSendCells.size()-1; procI >=1; procI--) +// { +// label nLocalCells = globalCells.localSize(procI); +// if (nLocalCells-nSendCells[procI] < 1) +// { +// nSendCells[procI-1] = nSendCells[procI]-nLocalCells+1; +// } +// } +// +// // First receive (so increasing the sizes of all arrays) +// +// Field xadj(initxadj); +// Field adjncy(initadjncy); +// scalarField cWeights(initcWeights); +// +// if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) +// { +// // Receive cells from previous processor +// IPstream fromPrevProc(Pstream::blocking, Pstream::myProcNo()-1); +// +// Field prevXadj(fromPrevProc); +// Field prevAdjncy(fromPrevProc); +// scalarField prevCellWeights(fromPrevProc); +// +// if (prevXadj.size() != nSendCells[Pstream::myProcNo()-1]) +// { +// FatalErrorIn("ptscotchDecomp::decompose(..)") +// << "Expected from processor " << Pstream::myProcNo()-1 +// << " connectivity for " << nSendCells[Pstream::myProcNo()-1] +// << " nCells but only received " << prevXadj.size() +// << abort(FatalError); +// } +// +// // Insert adjncy +// prepend(prevAdjncy, adjncy); +// // Adapt offsets and prepend xadj +// xadj += prevAdjncy.size(); +// prepend(prevXadj, xadj); +// // Weights +// prepend(prevCellWeights, cWeights); +// } +// +// +// // Send to my next processor +// +// if (nSendCells[Pstream::myProcNo()] > 0) +// { +// // Send cells to next processor +// OPstream toNextProc(Pstream::blocking, Pstream::myProcNo()+1); +// +// int nCells = nSendCells[Pstream::myProcNo()]; +// int startCell = xadj.size()-1 - nCells; +// int startFace = xadj[startCell]; +// int nFaces = adjncy.size()-startFace; +// +// // Send for all cell data: last nCells elements +// // Send for all face data: last nFaces elements +// toNextProc +// << Field::subField(xadj, nCells, startCell)-startFace +// << Field::subField(adjncy, nFaces, startFace) +// << +// ( +// cWeights.size() +// ? static_cast +// ( +// scalarField::subField(cWeights, nCells, startCell) +// ) +// : scalarField(0) +// ); +// +// // Remove data that has been sent +// if (cWeights.size()) +// { +// cWeights.setSize(cWeights.size()-nCells); +// } +// adjncy.setSize(adjncy.size()-nFaces); +// xadj.setSize(xadj.size() - nCells); +// } +// +// +// // Do decomposition as normal. Sets finalDecomp. +// label result = decompose(meshPath, adjncy, xadj, cWeights, finalDecomp); +// +// +// if (debug) +// { +// Info<< "ptscotchDecomp : have graphs with locally 0 cells." +// << " trickling up." << endl; +// } +// +// +// // If we sent cells across make sure we undo it +// // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +// +// // Receive back from next processor if I sent something +// if (nSendCells[Pstream::myProcNo()] > 0) +// { +// IPstream fromNextProc(Pstream::blocking, Pstream::myProcNo()+1); +// +// List nextFinalDecomp(fromNextProc); +// +// if (nextFinalDecomp.size() != nSendCells[Pstream::myProcNo()]) +// { +// FatalErrorIn("parMetisDecomp::decompose(..)") +// << "Expected from processor " << Pstream::myProcNo()+1 +// << " decomposition for " << nSendCells[Pstream::myProcNo()] +// << " nCells but only received " << nextFinalDecomp.size() +// << abort(FatalError); +// } +// +// append(nextFinalDecomp, finalDecomp); +// } +// +// // Send back to previous processor. +// if (Pstream::myProcNo() >= 1 && nSendCells[Pstream::myProcNo()-1] > 0) +// { +// OPstream toPrevProc(Pstream::blocking, Pstream::myProcNo()-1); +// +// int nToPrevious = nSendCells[Pstream::myProcNo()-1]; +// +// toPrevProc << +// SubList +// ( +// finalDecomp, +// nToPrevious, +// finalDecomp.size()-nToPrevious +// ); +// +// // Remove locally what has been sent +// finalDecomp.setSize(finalDecomp.size()-nToPrevious); +// } +// return result; +//} // Call scotch with options from dictionary. @@ -362,13 +362,42 @@ Foam::label Foam::ptscotchDecomp::decompose const List& adjncy, const List& xadj, const scalarField& cWeights, + List& finalDecomp +) const +{ + List dummyAdjncy(1); + List dummyXadj(1); + dummyXadj[0] = 0; + + return decompose + ( + meshPath, + adjncy.size(), + (adjncy.size() ? adjncy.begin() : dummyAdjncy.begin()), + xadj.size(), + (xadj.size() ? xadj.begin() : dummyXadj.begin()), + cWeights, + finalDecomp + ); +} + + +// Call scotch with options from dictionary. +Foam::label Foam::ptscotchDecomp::decompose +( + const fileName& meshPath, + const int adjncySize, + const int adjncy[], + const int xadjSize, + const int xadj[], + const scalarField& cWeights, List& finalDecomp ) const { if (debug) { - Pout<< "ptscotchDecomp : entering with xadj:" << xadj.size() << endl; + Pout<< "ptscotchDecomp : entering with xadj:" << xadjSize << endl; } // Dump graph @@ -387,7 +416,7 @@ Foam::label Foam::ptscotchDecomp::decompose Pout<< "Dumping Scotch graph file to " << str.name() << endl << "Use this in combination with dgpart." << endl; - globalIndex globalCells(xadj.size()-1); + globalIndex globalCells(xadjSize-1); // Distributed graph file (.grf) label version = 2; @@ -400,17 +429,17 @@ Foam::label Foam::ptscotchDecomp::decompose // Total number of vertices (vertglbnbr) str << globalCells.size(); // Total number of connections (edgeglbnbr) - str << ' ' << returnReduce(xadj[xadj.size()-1], sumOp