Merge branch 'master' of ssh://dm/home/dm4/OpenFOAM/OpenFOAM-dev

This commit is contained in:
Henry
2012-04-02 12:06:42 +01:00
26 changed files with 1531 additions and 1094 deletions

View File

@ -4,5 +4,7 @@ set -x
wclean libso conformalVoronoiMesh
wclean
wclean cvMeshSurfaceSimplify
wclean cvMeshBackgroundMesh
# ----------------------------------------------------------------- end-of-file

View File

@ -5,5 +5,6 @@ set -x
wmake libso conformalVoronoiMesh
wmake
wmake cvMeshSurfaceSimplify
wmake cvMeshBackgroundMesh
# ----------------------------------------------------------------- end-of-file

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -142,7 +142,7 @@ void Foam::backgroundMeshDecomposition::initialRefinement()
zeroGradientFvPatchScalarField::typeName
);
const conformationSurfaces& geometry = cvMesh_.geometryToConformTo();
const conformationSurfaces& geometry = geometryToConformTo_;
decompositionMethod& decomposer = decomposerPtr_();
@ -512,7 +512,7 @@ bool Foam::backgroundMeshDecomposition::refineCell
// Sample the box to find an estimate of the min size, and a volume
// estimate when overlapping == true.
const conformationSurfaces& geometry = cvMesh_.geometryToConformTo();
const conformationSurfaces& geometry = geometryToConformTo_;
treeBoundBox cellBb
(
@ -578,7 +578,7 @@ bool Foam::backgroundMeshDecomposition::refineCell
forAll(samplePoints, i)
{
scalar s = cvMesh_.cellSizeControl().cellSize
scalar s = cellSizeControl_.cellSize
(
hitInfo[i].hitPoint()
);
@ -693,7 +693,7 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
// Overall bb
treeBoundBox overallBb(boundaryFacesPtr_().localPoints());
Random& rnd = cvMesh_.rndGen();
Random& rnd = rndGen_;
bFTreePtr_.reset
(
@ -726,11 +726,11 @@ void Foam::backgroundMeshDecomposition::buildPatchAndTree()
octreeNearestDistances_ = bFTreePtr_().calcNearestDistance();
if (cvMesh_.cvMeshControls().objOutput())
if (cvMeshControls_.objOutput())
{
OFstream fStr
(
cvMesh_.time().path()
mesh_.time().path()
/"backgroundMeshDecomposition_proc_"
+ name(Pstream::myProcNo())
+ "_boundaryFaces.obj"
@ -777,15 +777,18 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
const conformalVoronoiMesh& cvMesh
)
:
coeffsDict_(coeffsDict),
cvMesh_(cvMesh),
runTime_(cvMesh.time()),
geometryToConformTo_(cvMesh.geometryToConformTo()),
cellSizeControl_(cvMesh.cellSizeControl()),
rndGen_(cvMesh.rndGen()),
cvMeshControls_(cvMesh.cvMeshControls()),
mesh_
(
IOobject
(
fvMesh::defaultRegion,
cvMesh_.time().timeName(),
cvMesh_.time(),
runTime_.timeName(),
runTime_,
IOobject::MUST_READ
)
),
@ -805,22 +808,22 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
IOobject
(
"decomposeParDict",
cvMesh_.time().system(),
cvMesh_.time(),
runTime_.system(),
runTime_,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
),
decomposerPtr_(decompositionMethod::New(decomposeDict_)),
mergeDist_(1e-6*mesh_.bounds().mag()),
spanScale_(readScalar(coeffsDict_.lookup("spanScale"))),
spanScale_(readScalar(coeffsDict.lookup("spanScale"))),
minCellSizeLimit_
(
coeffsDict_.lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
coeffsDict.lookupOrDefault<scalar>("minCellSizeLimit", 0.0)
),
minLevels_(readLabel(coeffsDict_.lookup("minLevels"))),
volRes_(readLabel(coeffsDict_.lookup("sampleResolution"))),
maxCellWeightCoeff_(readScalar(coeffsDict_.lookup("maxCellWeightCoeff")))
minLevels_(readLabel(coeffsDict.lookup("minLevels"))),
volRes_(readLabel(coeffsDict.lookup("sampleResolution"))),
maxCellWeightCoeff_(readScalar(coeffsDict.lookup("maxCellWeightCoeff")))
{
if (!Pstream::parRun())
{
@ -854,6 +857,74 @@ Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
}
Foam::backgroundMeshDecomposition::backgroundMeshDecomposition
(
const scalar spanScale,
const scalar minCellSizeLimit,
const label minLevels,
const label volRes,
const scalar maxCellWeightCoeff,
const Time& runTime,
const conformationSurfaces& geometryToConformTo,
const cellSizeControlSurfaces& cellSizeControl,
Random& rndGen,
const cvControls& cvMeshControls
)
:
runTime_(runTime),
geometryToConformTo_(geometryToConformTo),
cellSizeControl_(cellSizeControl),
rndGen_(rndGen),
cvMeshControls_(cvMeshControls),
mesh_
(
IOobject
(
fvMesh::defaultRegion,
runTime_.timeName(),
runTime_,
IOobject::MUST_READ
)
),
meshCutter_
(
mesh_,
labelList(mesh_.nCells(), 0),
labelList(mesh_.nPoints(), 0)
),
boundaryFacesPtr_(),
bFTreePtr_(),
octreeNearestDistances_(),
allBackgroundMeshBounds_(Pstream::nProcs()),
globalBackgroundBounds_(),
decomposeDict_
(
IOobject
(
"decomposeParDict",
runTime_.system(),
runTime_,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
),
decomposerPtr_(decompositionMethod::New(decomposeDict_)),
mergeDist_(1e-6*mesh_.bounds().mag()),
spanScale_(spanScale),
minCellSizeLimit_(minCellSizeLimit),
minLevels_(minLevels),
volRes_(volRes),
maxCellWeightCoeff_(maxCellWeightCoeff)
{
// Stand-alone operation
Info<< nl << "Building initial background mesh decomposition" << endl;
initialRefinement();
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::backgroundMeshDecomposition::~backgroundMeshDecomposition()

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -92,10 +92,22 @@ class backgroundMeshDecomposition
// Private data
//- Method details dictionary
dictionary coeffsDict_;
//dictionary coeffsDict_;
//- Reference to the conformalVoronoiMesh holding this object
const conformalVoronoiMesh& cvMesh_;
//- Reference to runtime
const Time& runTime_;
//- Reference to surface
const conformationSurfaces& geometryToConformTo_;
//- The cell size control object
const cellSizeControlSurfaces& cellSizeControl_;
//- Random number generator
Random& rndGen_;
//- Controls
const cvControls& cvMeshControls_;
//- Mesh stored on for this processor, specifiying the domain that it
// is responsible for.
@ -191,13 +203,29 @@ public:
// Constructors
//- Construct from components
//- Construct from components in cvMesh operation
backgroundMeshDecomposition
(
const dictionary& coeffsDict,
const conformalVoronoiMesh& cvMesh
);
//- Construct from components for standalone operation
backgroundMeshDecomposition
(
const scalar spanScale,
const scalar minCellSizeLimit,
const label minLevels,
const label volRes,
const scalar maxCellWeightCoeff,
const Time& runTime,
const conformationSurfaces& geometryToConformTo,
const cellSizeControlSurfaces& cellSizeControl,
Random& rndGen,
const cvControls& cvMeshControls
);
//- Destructor
~backgroundMeshDecomposition();
@ -299,6 +327,15 @@ public:
//- Return the boundBox of this processor
inline const treeBoundBox& procBounds() const;
//- Return the cell level of the underlying mesh
inline const labelList& cellLevel() const;
//- Return the point level of the underlying mesh
inline const labelList& pointLevel() const;
//- Return the current decomposition method
inline const decompositionMethod& decomposer() const;
};

View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -43,9 +43,30 @@ Foam::backgroundMeshDecomposition::octreeNearestDistances() const
}
const Foam::treeBoundBox& Foam::backgroundMeshDecomposition::procBounds() const
const Foam::treeBoundBox&
Foam::backgroundMeshDecomposition::procBounds() const
{
return allBackgroundMeshBounds_[Pstream::myProcNo()];
}
const Foam::labelList& Foam::backgroundMeshDecomposition::cellLevel() const
{
return meshCutter_.cellLevel();
}
const Foam::labelList& Foam::backgroundMeshDecomposition::pointLevel() const
{
return meshCutter_.pointLevel();
}
const Foam::decompositionMethod&
Foam::backgroundMeshDecomposition::decomposer() const
{
return decomposerPtr_();
}
// ************************************************************************* //

View File

@ -0,0 +1,3 @@
cvMeshBackgroundMesh.C
EXE = $(FOAM_APPBIN)/cvMeshBackgroundMesh

View File

@ -0,0 +1,29 @@
EXE_DEBUG = -DFULLDEBUG -g -O0
EXE_FROUNDING_MATH = -frounding-math
EXE_NDEBUG = -DNDEBUG
include $(GENERAL_RULES)/CGAL
EXE_INC = \
${EXE_FROUNDING_MATH} \
${EXE_NDEBUG} \
${CGAL_INC} \
-I../conformalVoronoiMesh/lnInclude \
-I$(LIB_SRC)/parallel/decompose/decompositionMethods/lnInclude \
-I$(LIB_SRC)/edgeMesh/lnInclude \
-I$(LIB_SRC)/sampling/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/dynamicMesh/lnInclude
EXE_LIBS = \
$(CGAL_LIBS) \
-lconformalVoronoiMesh \
-ldecompositionMethods /* -L$(FOAM_LIBBIN)/dummy -lscotchDecomp */ \
-ledgeMesh \
-lsampling \
-ltriSurface \
-lmeshTools \
-ldynamicMesh \
-lfiniteVolume

View File

@ -0,0 +1,770 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ 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 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/>.
Application
cvMeshBackGroundMesh
Description
Writes out background mesh as constructed by cvMesh and constructs
distanceSurface.
\*---------------------------------------------------------------------------*/
#include "PatchTools.H"
#include "argList.H"
#include "Time.H"
#include "triSurface.H"
#include "searchableSurfaces.H"
#include "conformationSurfaces.H"
#include "cellSizeControlSurfaces.H"
#include "backgroundMeshDecomposition.H"
#include "cellShape.H"
#include "cellModeller.H"
#include "DynamicField.H"
#include "isoSurfaceCell.H"
#include "vtkSurfaceWriter.H"
#include "syncTools.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-6;
// Get merging distance when matching face centres
scalar getMergeDistance
(
const argList& args,
const Time& runTime,
const boundBox& bb
)
{
scalar mergeTol = defaultMergeTol;
args.optionReadIfPresent("mergeTol", mergeTol);
scalar writeTol =
Foam::pow(scalar(10.0), -scalar(IOstream::defaultPrecision()));
Info<< "Merge tolerance : " << mergeTol << nl
<< "Write tolerance : " << writeTol << endl;
if (runTime.writeFormat() == IOstream::ASCII && mergeTol < writeTol)
{
FatalErrorIn("getMergeDistance")
<< "Your current settings specify ASCII writing with "
<< IOstream::defaultPrecision() << " digits precision." << endl
<< "Your merging tolerance (" << mergeTol << ") is finer than this."
<< endl
<< "Please change your writeFormat to binary"
<< " or increase the writePrecision" << endl
<< "or adjust the merge tolerance (-mergeTol)."
<< exit(FatalError);
}
scalar mergeDist = mergeTol * bb.mag();
Info<< "Overall meshes bounding box : " << bb << nl
<< "Relative tolerance : " << mergeTol << nl
<< "Absolute matching distance : " << mergeDist << nl
<< endl;
return mergeDist;
}
void printMeshData(const polyMesh& mesh)
{
// Collect all data on master
globalIndex globalCells(mesh.nCells());
labelListList patchNeiProcNo(Pstream::nProcs());
labelListList patchSize(Pstream::nProcs());
const labelList& pPatches = mesh.globalData().processorPatches();
patchNeiProcNo[Pstream::myProcNo()].setSize(pPatches.size());
patchSize[Pstream::myProcNo()].setSize(pPatches.size());
forAll(pPatches, i)
{
const processorPolyPatch& ppp = refCast<const processorPolyPatch>
(
mesh.boundaryMesh()[pPatches[i]]
);
patchNeiProcNo[Pstream::myProcNo()][i] = ppp.neighbProcNo();
patchSize[Pstream::myProcNo()][i] = ppp.size();
}
Pstream::gatherList(patchNeiProcNo);
Pstream::gatherList(patchSize);
// Print stats
globalIndex globalBoundaryFaces(mesh.nFaces()-mesh.nInternalFaces());
label maxProcCells = 0;
label totProcFaces = 0;
label maxProcPatches = 0;
label totProcPatches = 0;
label maxProcFaces = 0;
for (label procI = 0; procI < Pstream::nProcs(); procI++)
{
Info<< endl
<< "Processor " << procI << nl
<< " Number of cells = " << globalCells.localSize(procI)
<< endl;
label nProcFaces = 0;
const labelList& nei = patchNeiProcNo[procI];
forAll(patchNeiProcNo[procI], i)
{
Info<< " Number of faces shared with processor "
<< patchNeiProcNo[procI][i] << " = " << patchSize[procI][i]
<< endl;
nProcFaces += patchSize[procI][i];
}
Info<< " Number of processor patches = " << nei.size() << nl
<< " Number of processor faces = " << nProcFaces << nl
<< " Number of boundary faces = "
<< globalBoundaryFaces.localSize(procI) << endl;
maxProcCells = max(maxProcCells, globalCells.localSize(procI));
totProcFaces += nProcFaces;
totProcPatches += nei.size();
maxProcPatches = max(maxProcPatches, nei.size());
maxProcFaces = max(maxProcFaces, nProcFaces);
}
// Stats
scalar avgProcCells = scalar(globalCells.size())/Pstream::nProcs();
scalar avgProcPatches = scalar(totProcPatches)/Pstream::nProcs();
scalar avgProcFaces = scalar(totProcFaces)/Pstream::nProcs();
// In case of all faces on one processor. Just to avoid division by 0.
if (totProcPatches == 0)
{
avgProcPatches = 1;
}
if (totProcFaces == 0)
{
avgProcFaces = 1;
}
Info<< nl
<< "Number of processor faces = " << totProcFaces/2 << nl
<< "Max number of cells = " << maxProcCells
<< " (" << 100.0*(maxProcCells-avgProcCells)/avgProcCells
<< "% above average " << avgProcCells << ")" << nl
<< "Max number of processor patches = " << maxProcPatches
<< " (" << 100.0*(maxProcPatches-avgProcPatches)/avgProcPatches
<< "% above average " << avgProcPatches << ")" << nl
<< "Max number of faces between processors = " << maxProcFaces
<< " (" << 100.0*(maxProcFaces-avgProcFaces)/avgProcFaces
<< "% above average " << avgProcFaces << ")" << nl
<< endl;
}
// Return cellID
label cellLabel
(
const Vector<label>& nCells,
const label i,
const label j,
const label k
)
{
return i*nCells[1]*nCells[2]+j*nCells[2]+k;
}
label vtxLabel
(
const Vector<label>& nCells,
const label i,
const label j,
const label k
)
{
Vector<label> nPoints
(
nCells[0]+1,
nCells[1]+1,
nCells[2]+1
);
return i*nPoints[1]*nPoints[2]+j*nPoints[2]+k;
}
autoPtr<polyMesh> generateHexMesh
(
const IOobject& io,
const point& origin,
const vector& cellSize,
const Vector<label>& nCells
)
{
pointField points;
if (nCells[0]+nCells[1]+nCells[2] > 0)
{
points.setSize((nCells[0]+1)*(nCells[1]+1)*(nCells[2]+1));
// Generate points
for (label i = 0; i <= nCells[0]; i++)
{
for (label j = 0; j <= nCells[1]; j++)
{
for (label k = 0; k <= nCells[2]; k++)
{
point pt = origin;
pt.x() += i*cellSize[0];
pt.y() += j*cellSize[1];
pt.z() += k*cellSize[2];
points[vtxLabel(nCells, i, j, k)] = pt;
}
}
}
}
const cellModel& hex = *(cellModeller::lookup("hex"));
cellShapeList cellShapes(nCells[0]*nCells[1]*nCells[2]);
labelList hexPoints(8);
label cellI = 0;
for (label i = 0; i < nCells[0]; i++)
{
for (label j = 0; j < nCells[1]; j++)
{
for (label k = 0; k < nCells[2]; k++)
{
hexPoints[0] = vtxLabel(nCells, i, j, k);
hexPoints[1] = vtxLabel(nCells, i+1, j, k);
hexPoints[2] = vtxLabel(nCells, i+1, j+1, k);
hexPoints[3] = vtxLabel(nCells, i, j+1, k);
hexPoints[4] = vtxLabel(nCells, i, j, k+1);
hexPoints[5] = vtxLabel(nCells, i+1, j, k+1);
hexPoints[6] = vtxLabel(nCells, i+1, j+1, k+1);
hexPoints[7] = vtxLabel(nCells, i, j+1, k+1);
cellShapes[cellI++] = cellShape(hex, hexPoints);
}
}
}
faceListList boundary(0);
wordList patchNames(0);
wordList patchTypes(0);
word defaultFacesName = "defaultFaces";
word defaultFacesType = polyPatch::typeName;
wordList patchPhysicalTypes(0);
return autoPtr<polyMesh>
(
new polyMesh
(
io,
xferMoveTo<pointField>(points),
cellShapes,
boundary,
patchNames,
patchTypes,
defaultFacesName,
defaultFacesType,
patchPhysicalTypes
)
);
}
// Determine for every point a signed distance to the nearest surface
// (outside is positive)
tmp<scalarField> signedDistance
(
const scalarField& distSqr,
const pointField& points,
const searchableSurfaces& geometry,
const labelList& surfaces
)
{
tmp<scalarField> tfld(new scalarField(points.size(), Foam::sqr(GREAT)));
scalarField& fld = tfld();
// Find nearest
List<pointIndexHit> nearest;
labelList nearestSurfaces;
searchableSurfacesQueries::findNearest
(
geometry,
surfaces,
points,
scalarField(points.size(), Foam::sqr(GREAT)),//distSqr
nearestSurfaces,
nearest
);
// Determine sign of nearest. Sort by surface to do this.
DynamicField<point> surfPoints(points.size());
DynamicList<label> surfIndices(points.size());
forAll(surfaces, surfI)
{
// Extract points on this surface
surfPoints.clear();
surfIndices.clear();
forAll(nearestSurfaces, i)
{
if (nearestSurfaces[i] == surfI)
{
surfPoints.append(points[i]);
surfIndices.append(i);
}
}
// Calculate sideness of these surface points
label geomI = surfaces[surfI];
List<searchableSurface::volumeType> volType;
geometry[geomI].getVolumeType(surfPoints, volType);
// Push back to original
forAll(volType, i)
{
label pointI = surfIndices[i];
scalar dist = mag(points[pointI] - nearest[pointI].hitPoint());
searchableSurface::volumeType vT = volType[i];
if (vT == searchableSurface::OUTSIDE)
{
fld[pointI] = dist;
}
else if (vT == searchableSurface::INSIDE)
{
fld[i] = -dist;
}
else
{
FatalErrorIn("signedDistance()")
<< "getVolumeType failure, neither INSIDE or OUTSIDE"
<< exit(FatalError);
}
}
}
return tfld;
}
// Main program:
int main(int argc, char *argv[])
{
argList::addNote
(
"Generate cvMesh-consistent representation of surfaces"
);
argList::addBoolOption
(
"writeMesh",
"write the resulting mesh and distance fields"
);
argList::addOption
(
"mergeTol",
"scalar",
"specify the merge distance relative to the bounding box size "
"(default 1E-6)"
);
#include "setRootCase.H"
#include "createTime.H"
runTime.functionObjects().off();
const bool writeMesh = args.optionFound("writeMesh");
if (writeMesh)
{
Info<< "Writing resulting mesh and cellDistance, pointDistance fields."
<< nl << endl;
}
IOdictionary cvMeshDict
(
IOobject
(
"cvMeshDict",
runTime.system(),
runTime,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
// Define/load all geometry
searchableSurfaces allGeometry
(
IOobject
(
"cvSearchableSurfaces",
runTime.constant(),
"triSurface",
runTime,
IOobject::MUST_READ,
IOobject::NO_WRITE
),
cvMeshDict.subDict("geometry")
);
Random rndGen(64293*Pstream::myProcNo());
conformationSurfaces geometryToConformTo
(
runTime,
rndGen,
allGeometry,
cvMeshDict.subDict("surfaceConformation")
);
cellSizeControlSurfaces cellSizeControl
(
allGeometry,
cvMeshDict.subDict("motionControl")
);
// Generate starting block mesh
vector cellSize;
{
const treeBoundBox& bb = geometryToConformTo.globalBounds();
// Determine the number of cells in each direction.
const vector span = bb.span();
vector nScalarCells = span/cellSizeControl.defaultCellSize();
// Calculate initial cell size to be a little bit smaller than the
// defaultCellSize to avoid initial refinement triggering.
Vector<label> nCells = Vector<label>
(
label(nScalarCells.x())+2,
label(nScalarCells.y())+2,
label(nScalarCells.z())+2
);
cellSize = vector
(
span[0]/nCells[0],
span[1]/nCells[1],
span[2]/nCells[2]
);
Info<< "Generating initial hex mesh with" << nl
<< " bounding box : " << bb << nl
<< " nCells : " << nCells << nl
<< " cellSize : " << cellSize << nl
<< endl;
autoPtr<polyMesh> meshPtr
(
generateHexMesh
(
IOobject
(
polyMesh::defaultRegion,
runTime.constant(),
runTime
),
bb.min(),
cellSize,
(
Pstream::master()
? nCells
: Vector<label>(0, 0, 0)
)
)
);
Info<< "Writing initial hex mesh to " << meshPtr().instance() << nl
<< endl;
meshPtr().write();
}
// Distribute the initial mesh
if (Pstream::parRun())
{
# include "createMesh.H"
Info<< "Loaded mesh:" << endl;
printMeshData(mesh);
// Allocate a decomposer
IOdictionary decompositionDict
(
IOobject
(
"decomposeParDict",
runTime.system(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
autoPtr<decompositionMethod> decomposer
(
decompositionMethod::New
(
decompositionDict
)
);
labelList decomp = decomposer().decompose(mesh, mesh.cellCentres());
// Global matching tolerance
const scalar tolDim = getMergeDistance
(
args,
runTime,
mesh.bounds()
);
// Mesh distribution engine
fvMeshDistribute distributor(mesh, tolDim);
Info<< "Wanted distribution:"
<< distributor.countCells(decomp) << nl << endl;
// Do actual sending/receiving of mesh
autoPtr<mapDistributePolyMesh> map = distributor.distribute(decomp);
// Print some statistics
//Info<< "After distribution:" << endl;
//printMeshData(mesh);
mesh.setInstance(runTime.constant());
Info<< "Writing redistributed mesh" << nl << endl;
mesh.write();
}
Info<< "Refining backgroud mesh according to cell size specification" << nl
<< endl;
backgroundMeshDecomposition backgroundMesh
(
1.0, //spanScale,ratio of poly cell size v.s. hex cell size
0.0, //minCellSizeLimit
0, //minLevels
4, //volRes, check multiple points per cell
20.0, //maxCellWeightCoeff
runTime,
geometryToConformTo,
cellSizeControl,
rndGen,
cvMeshDict
);
if (writeMesh)
{
runTime++;
Info<< "Writing mesh to " << runTime.timeName() << endl;
backgroundMesh.mesh().write();
}
const scalar tolDim = getMergeDistance
(
args,
runTime,
backgroundMesh.mesh().bounds()
);
faceList isoFaces;
pointField isoPoints;
{
// Apply a distanceSurface to it.
const fvMesh& fvm = backgroundMesh.mesh();
volScalarField cellDistance
(
IOobject
(
"cellDistance",
fvm.time().timeName(),
fvm.time(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
fvm,
dimensionedScalar("zero", dimLength, 0)
);
const searchableSurfaces& geometry = geometryToConformTo.geometry();
const labelList& surfaces = geometryToConformTo.surfaces();
// Get maximum search size per cell
scalarField distSqr(cellDistance.size());
const labelList& cellLevel = backgroundMesh.cellLevel();
forAll(cellLevel, cellI)
{
// The largest edge of the cell will always be less than the
// span of the bounding box of the cell.
distSqr[cellI] = magSqr(cellSize)/pow(2, cellLevel[cellI]);
}
{
// Internal field
cellDistance.internalField() = signedDistance
(
distSqr,
fvm.C(),
geometry,
surfaces
);
// Patch fields
forAll(fvm.C().boundaryField(), patchI)
{
const pointField& cc = fvm.C().boundaryField()[patchI];
fvPatchScalarField& fld = cellDistance.boundaryField()[patchI];
scalarField patchDistSqr
(
fld.patch().patchInternalField(distSqr)
);
fld = signedDistance(patchDistSqr, cc, geometry, surfaces);
}
// On processor patches the fvm.C() will already be the cell centre
// on the opposite side so no need to swap cellDistance.
if (writeMesh)
{
cellDistance.write();
}
}
// Distance to points
pointScalarField pointDistance
(
IOobject
(
"pointDistance",
fvm.time().timeName(),
fvm.time(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
pointMesh::New(fvm),
dimensionedScalar("zero", dimLength, 0)
);
{
scalarField pointDistSqr(fvm.nPoints(), -sqr(GREAT));
for (label faceI = 0; faceI < fvm.nInternalFaces(); faceI++)
{
label own = fvm.faceOwner()[faceI];
label ownDistSqr = distSqr[own];
const face& f = fvm.faces()[faceI];
forAll(f, fp)
{
pointDistSqr[f[fp]] = max(pointDistSqr[f[fp]], ownDistSqr);
}
}
syncTools::syncPointList
(
fvm,
pointDistSqr,
maxEqOp<scalar>(),
-sqr(GREAT) // null value
);
pointDistance.internalField() = signedDistance
(
pointDistSqr,
fvm.points(),
geometry,
surfaces
);
if (writeMesh)
{
pointDistance.write();
}
}
isoSurfaceCell iso
(
fvm,
cellDistance,
pointDistance,
0, //distance,
false //regularise
);
isoFaces.setSize(iso.size());
forAll(isoFaces, i)
{
isoFaces[i] = iso[i].triFaceFace();
}
isoPoints = iso.points();
}
pointField mergedPoints;
faceList mergedFaces;
labelList pointMergeMap;
PatchTools::gatherAndMerge
(
tolDim,
primitivePatch
(
SubList<face>(isoFaces, isoFaces.size()),
isoPoints
),
mergedPoints,
mergedFaces,
pointMergeMap
);
vtkSurfaceWriter writer;
writer.write
(
runTime.path(),
"iso",
mergedPoints,
mergedFaces
);
Info<< "End\n" << endl;
return 0;
}
// ************************************************************************* //

View File

@ -1,3 +1,4 @@
loadOrCreateMesh.C
redistributePar.C
EXE = $(FOAM_APPBIN)/redistributePar

View File

@ -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
@ -54,6 +54,7 @@ Description
#include "mapDistributePolyMesh.H"
#include "IOobjectList.H"
#include "globalIndex.H"
#include "loadOrCreateMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -62,292 +63,292 @@ Description
static const scalar defaultMergeTol = 1E-6;
// Read mesh if available. Otherwise create empty mesh with same non-proc
// patches as proc0 mesh. Requires all processors to have all patches
// (and in same order).
autoPtr<fvMesh> createMesh
(
const Time& runTime,
const word& regionName,
const fileName& instDir,
const bool haveMesh
)
{
//Pout<< "Create mesh for time = "
// << runTime.timeName() << nl << endl;
IOobject io
(
regionName,
instDir,
runTime,
IOobject::MUST_READ
);
if (!haveMesh)
{
// Create dummy mesh. Only used on procs that don't have mesh.
IOobject noReadIO(io);
noReadIO.readOpt() = IOobject::NO_READ;
fvMesh dummyMesh
(
noReadIO,
xferCopy(pointField()),
xferCopy(faceList()),
xferCopy(labelList()),
xferCopy(labelList()),
false
);
// Add some dummy zones so upon reading it does not read them
// from the undecomposed case. Should be done as extra argument to
// regIOobject::readStream?
List<pointZone*> pz
(
1,
new pointZone
(
"dummyPointZone",
labelList(0),
0,
dummyMesh.pointZones()
)
);
List<faceZone*> fz
(
1,
new faceZone
(
"dummyFaceZone",
labelList(0),
boolList(0),
0,
dummyMesh.faceZones()
)
);
List<cellZone*> cz
(
1,
new cellZone
(
"dummyCellZone",
labelList(0),
0,
dummyMesh.cellZones()
)
);
dummyMesh.addZones(pz, fz, cz);
//Pout<< "Writing dummy mesh to " << dummyMesh.polyMesh::objectPath()
// << endl;
dummyMesh.write();
}
//Pout<< "Reading mesh from " << io.objectPath() << endl;
autoPtr<fvMesh> meshPtr(new fvMesh(io));
fvMesh& mesh = meshPtr();
// Sync patches
// ~~~~~~~~~~~~
if (Pstream::master())
{
// Send patches
for
(
int slave=Pstream::firstSlave();
slave<=Pstream::lastSlave();
slave++
)
{
OPstream toSlave(Pstream::scheduled, slave);
toSlave << mesh.boundaryMesh();
}
}
else
{
// Receive patches
IPstream fromMaster(Pstream::scheduled, Pstream::masterNo());
PtrList<entry> patchEntries(fromMaster);
if (haveMesh)
{
// Check master names against mine
const polyBoundaryMesh& patches = mesh.boundaryMesh();
forAll(patchEntries, patchI)
{
const entry& e = patchEntries[patchI];
const word type(e.dict().lookup("type"));
const word& name = e.keyword();
if (type == processorPolyPatch::typeName)
{
break;
}
if (patchI >= patches.size())
{
FatalErrorIn
(
"createMesh(const Time&, const fileName&, const bool)"
) << "Non-processor patches not synchronised."
<< endl
<< "Processor " << Pstream::myProcNo()
<< " has only " << patches.size()
<< " patches, master has "
<< patchI
<< exit(FatalError);
}
if
(
type != patches[patchI].type()
|| name != patches[patchI].name()
)
{
FatalErrorIn
(
"createMesh(const Time&, const fileName&, const bool)"
) << "Non-processor patches not synchronised."
<< endl
<< "Master patch " << patchI
<< " name:" << type
<< " type:" << type << endl
<< "Processor " << Pstream::myProcNo()
<< " patch " << patchI
<< " has name:" << patches[patchI].name()
<< " type:" << patches[patchI].type()
<< exit(FatalError);
}
}
}
else
{
// Add patch
List<polyPatch*> patches(patchEntries.size());
label nPatches = 0;
forAll(patchEntries, patchI)
{
const entry& e = patchEntries[patchI];
const word type(e.dict().lookup("type"));
const word& name = e.keyword();
if (type == processorPolyPatch::typeName)
{
break;
}
//Pout<< "Adding patch:" << nPatches
// << " name:" << name << " type:" << type << endl;
dictionary patchDict(e.dict());
patchDict.remove("nFaces");
patchDict.add("nFaces", 0);
patchDict.remove("startFace");
patchDict.add("startFace", 0);
patches[patchI] = polyPatch::New
(
name,
patchDict,
nPatches++,
mesh.boundaryMesh()
).ptr();
}
patches.setSize(nPatches);
mesh.addFvPatches(patches, false); // no parallel comms
//// Write empty mesh now we have correct patches
//meshPtr().write();
}
}
// Determine zones
// ~~~~~~~~~~~~~~~
wordList pointZoneNames(mesh.pointZones().names());
Pstream::scatter(pointZoneNames);
wordList faceZoneNames(mesh.faceZones().names());
Pstream::scatter(faceZoneNames);
wordList cellZoneNames(mesh.cellZones().names());
Pstream::scatter(cellZoneNames);
if (!haveMesh)
{
// Add the zones. Make sure to remove the old dummy ones first
mesh.pointZones().clear();
mesh.faceZones().clear();
mesh.cellZones().clear();
List<pointZone*> pz(pointZoneNames.size());
forAll(pointZoneNames, i)
{
pz[i] = new pointZone
(
pointZoneNames[i],
labelList(0),
i,
mesh.pointZones()
);
}
List<faceZone*> fz(faceZoneNames.size());
forAll(faceZoneNames, i)
{
fz[i] = new faceZone
(
faceZoneNames[i],
labelList(0),
boolList(0),
i,
mesh.faceZones()
);
}
List<cellZone*> cz(cellZoneNames.size());
forAll(cellZoneNames, i)
{
cz[i] = new cellZone
(
cellZoneNames[i],
labelList(0),
i,
mesh.cellZones()
);
}
mesh.addZones(pz, fz, cz);
}
if (!haveMesh)
{
// We created a dummy mesh file above. Delete it.
//Pout<< "Removing dummy mesh " << io.objectPath() << endl;
rmDir(io.objectPath());
}
// Force recreation of globalMeshData.
mesh.clearOut();
mesh.globalData();
// Do some checks.
// Check if the boundary definition is unique
mesh.boundaryMesh().checkDefinition(true);
// Check if the boundary processor patches are correct
mesh.boundaryMesh().checkParallelSync(true);
// Check names of zones are equal
mesh.cellZones().checkDefinition(true);
mesh.cellZones().checkParallelSync(true);
mesh.faceZones().checkDefinition(true);
mesh.faceZones().checkParallelSync(true);
mesh.pointZones().checkDefinition(true);
mesh.pointZones().checkParallelSync(true);
return meshPtr;
}
//// Read mesh if available. Otherwise create empty mesh with same non-proc
//// patches as proc0 mesh. Requires all processors to have all patches
//// (and in same order).
//autoPtr<fvMesh> createMesh
//(
// const Time& runTime,
// const word& regionName,
// const fileName& instDir,
// const bool haveMesh
//)
//{
// //Pout<< "Create mesh for time = "
// // << runTime.timeName() << nl << endl;
//
// IOobject io
// (
// regionName,
// instDir,
// runTime,
// IOobject::MUST_READ
// );
//
// if (!haveMesh)
// {
// // Create dummy mesh. Only used on procs that don't have mesh.
// IOobject noReadIO(io);
// noReadIO.readOpt() = IOobject::NO_READ;
// fvMesh dummyMesh
// (
// noReadIO,
// xferCopy(pointField()),
// xferCopy(faceList()),
// xferCopy(labelList()),
// xferCopy(labelList()),
// false
// );
// // Add some dummy zones so upon reading it does not read them
// // from the undecomposed case. Should be done as extra argument to
// // regIOobject::readStream?
// List<pointZone*> pz
// (
// 1,
// new pointZone
// (
// "dummyPointZone",
// labelList(0),
// 0,
// dummyMesh.pointZones()
// )
// );
// List<faceZone*> fz
// (
// 1,
// new faceZone
// (
// "dummyFaceZone",
// labelList(0),
// boolList(0),
// 0,
// dummyMesh.faceZones()
// )
// );
// List<cellZone*> cz
// (
// 1,
// new cellZone
// (
// "dummyCellZone",
// labelList(0),
// 0,
// dummyMesh.cellZones()
// )
// );
// dummyMesh.addZones(pz, fz, cz);
// //Pout<< "Writing dummy mesh to " << dummyMesh.polyMesh::objectPath()
// // << endl;
// dummyMesh.write();
// }
//
// //Pout<< "Reading mesh from " << io.objectPath() << endl;
// autoPtr<fvMesh> meshPtr(new fvMesh(io));
// fvMesh& mesh = meshPtr();
//
//
// // Sync patches
// // ~~~~~~~~~~~~
//
// if (Pstream::master())
// {
// // Send patches
// for
// (
// int slave=Pstream::firstSlave();
// slave<=Pstream::lastSlave();
// slave++
// )
// {
// OPstream toSlave(Pstream::scheduled, slave);
// toSlave << mesh.boundaryMesh();
// }
// }
// else
// {
// // Receive patches
// IPstream fromMaster(Pstream::scheduled, Pstream::masterNo());
// PtrList<entry> patchEntries(fromMaster);
//
// if (haveMesh)
// {
// // Check master names against mine
//
// const polyBoundaryMesh& patches = mesh.boundaryMesh();
//
// forAll(patchEntries, patchI)
// {
// const entry& e = patchEntries[patchI];
// const word type(e.dict().lookup("type"));
// const word& name = e.keyword();
//
// if (type == processorPolyPatch::typeName)
// {
// break;
// }
//
// if (patchI >= patches.size())
// {
// FatalErrorIn
// (
// "createMesh(const Time&, const fileName&, const bool)"
// ) << "Non-processor patches not synchronised."
// << endl
// << "Processor " << Pstream::myProcNo()
// << " has only " << patches.size()
// << " patches, master has "
// << patchI
// << exit(FatalError);
// }
//
// if
// (
// type != patches[patchI].type()
// || name != patches[patchI].name()
// )
// {
// FatalErrorIn
// (
// "createMesh(const Time&, const fileName&, const bool)"
// ) << "Non-processor patches not synchronised."
// << endl
// << "Master patch " << patchI
// << " name:" << type
// << " type:" << type << endl
// << "Processor " << Pstream::myProcNo()
// << " patch " << patchI
// << " has name:" << patches[patchI].name()
// << " type:" << patches[patchI].type()
// << exit(FatalError);
// }
// }
// }
// else
// {
// // Add patch
// List<polyPatch*> patches(patchEntries.size());
// label nPatches = 0;
//
// forAll(patchEntries, patchI)
// {
// const entry& e = patchEntries[patchI];
// const word type(e.dict().lookup("type"));
// const word& name = e.keyword();
//
// if (type == processorPolyPatch::typeName)
// {
// break;
// }
//
// //Pout<< "Adding patch:" << nPatches
// // << " name:" << name << " type:" << type << endl;
//
// dictionary patchDict(e.dict());
// patchDict.remove("nFaces");
// patchDict.add("nFaces", 0);
// patchDict.remove("startFace");
// patchDict.add("startFace", 0);
//
// patches[patchI] = polyPatch::New
// (
// name,
// patchDict,
// nPatches++,
// mesh.boundaryMesh()
// ).ptr();
// }
// patches.setSize(nPatches);
// mesh.addFvPatches(patches, false); // no parallel comms
//
// //// Write empty mesh now we have correct patches
// //meshPtr().write();
// }
// }
//
//
// // Determine zones
// // ~~~~~~~~~~~~~~~
//
// wordList pointZoneNames(mesh.pointZones().names());
// Pstream::scatter(pointZoneNames);
// wordList faceZoneNames(mesh.faceZones().names());
// Pstream::scatter(faceZoneNames);
// wordList cellZoneNames(mesh.cellZones().names());
// Pstream::scatter(cellZoneNames);
//
// if (!haveMesh)
// {
// // Add the zones. Make sure to remove the old dummy ones first
// mesh.pointZones().clear();
// mesh.faceZones().clear();
// mesh.cellZones().clear();
//
// List<pointZone*> pz(pointZoneNames.size());
// forAll(pointZoneNames, i)
// {
// pz[i] = new pointZone
// (
// pointZoneNames[i],
// labelList(0),
// i,
// mesh.pointZones()
// );
// }
// List<faceZone*> fz(faceZoneNames.size());
// forAll(faceZoneNames, i)
// {
// fz[i] = new faceZone
// (
// faceZoneNames[i],
// labelList(0),
// boolList(0),
// i,
// mesh.faceZones()
// );
// }
// List<cellZone*> cz(cellZoneNames.size());
// forAll(cellZoneNames, i)
// {
// cz[i] = new cellZone
// (
// cellZoneNames[i],
// labelList(0),
// i,
// mesh.cellZones()
// );
// }
// mesh.addZones(pz, fz, cz);
// }
//
//
// if (!haveMesh)
// {
// // We created a dummy mesh file above. Delete it.
// //Pout<< "Removing dummy mesh " << io.objectPath() << endl;
// rmDir(io.objectPath());
// }
//
// // Force recreation of globalMeshData.
// mesh.clearOut();
// mesh.globalData();
//
//
// // Do some checks.
//
// // Check if the boundary definition is unique
// mesh.boundaryMesh().checkDefinition(true);
// // Check if the boundary processor patches are correct
// mesh.boundaryMesh().checkParallelSync(true);
// // Check names of zones are equal
// mesh.cellZones().checkDefinition(true);
// mesh.cellZones().checkParallelSync(true);
// mesh.faceZones().checkDefinition(true);
// mesh.faceZones().checkParallelSync(true);
// mesh.pointZones().checkDefinition(true);
// mesh.pointZones().checkParallelSync(true);
//
// return meshPtr;
//}
// Get merging distance when matching face centres
@ -786,13 +787,24 @@ int main(int argc, char *argv[])
const bool allHaveMesh = (findIndex(haveMesh, false) == -1);
// Create mesh
autoPtr<fvMesh> meshPtr = createMesh
//autoPtr<fvMesh> meshPtr = createMesh
//(
// runTime,
// regionName,
// masterInstDir,
// haveMesh[Pstream::myProcNo()]
//);
autoPtr<fvMesh> meshPtr = loadOrCreateMesh
(
runTime,
regionName,
masterInstDir,
haveMesh[Pstream::myProcNo()]
IOobject
(
regionName,
masterInstDir,
runTime,
Foam::IOobject::MUST_READ
)
);
fvMesh& mesh = meshPtr();
// Print some statistics

View File

@ -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
@ -216,7 +216,6 @@ bool intersectSurfaces
int main(int argc, char *argv[])
{
argList::validArgs.clear();
argList::noParallel();
argList::validArgs.append("action");
argList::validArgs.append("surface file");

View File

@ -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
@ -106,6 +106,7 @@ AccessType ListListOps::combineOffset
return result;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -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
@ -92,10 +92,8 @@ SourceFiles
namespace Foam
{
//
//template <class T> class accessOp;
//template <class T> class offsetOp;
// Dummy access operator for combine()
// Dummy access operator for ListListOps::combine()
template <class T>
class accessOp
{
@ -108,7 +106,7 @@ public:
};
// Offset operator for combineOffset()
// Offset operator for ListListOps::combineOffset()
template <class T>
class offsetOp
{

View File

@ -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
@ -50,6 +50,7 @@ SourceFiles
#include "faceListFwd.H"
#include "intersection.H"
#include "pointHit.H"
#include "ListListOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -386,6 +387,30 @@ public:
};
//- Hash specialization to offset faces in ListListOps::combineOffset
template<>
class offsetOp<face>
{
public:
inline face operator()
(
const face& x,
const label offset
) const
{
face result(x.size());
forAll(x, xI)
{
result[xI] = x[xI] + offset;
}
return result;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -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
@ -46,6 +46,7 @@ SourceFiles
#include "intersection.H"
#include "pointField.H"
#include "triPointRef.H"
#include "ListListOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -249,6 +250,30 @@ template<>
inline bool contiguous<triFace>() {return true;}
//- Hash specialization to offset faces in ListListOps::combineOffset
template<>
class offsetOp<triFace>
{
public:
inline triFace operator()
(
const triFace& x,
const label offset
) const
{
triFace result(x);
forAll(x, xI)
{
result[xI] = x[xI] + offset;
}
return result;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam

View File

@ -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
@ -614,10 +614,17 @@ Foam::polyMesh::polyMesh
label nAllPatches = boundaryFaces.size();
if (nFaces > defaultPatchStart)
label nDefaultFaces = nFaces - defaultPatchStart;
if (syncPar)
{
reduce(nDefaultFaces, sumOp<label>());
}
if (nDefaultFaces > 0)
{
WarningIn("polyMesh::polyMesh(... construct from shapes...)")
<< "Found " << nFaces - defaultPatchStart
<< "Found " << nDefaultFaces
<< " undefined faces in mesh; adding to default patch." << endl;
// Check if there already exists a defaultFaces patch as last patch
@ -883,10 +890,16 @@ Foam::polyMesh::polyMesh
label nAllPatches = boundaryFaces.size();
if (nFaces > defaultPatchStart)
label nDefaultFaces = nFaces - defaultPatchStart;
if (syncPar)
{
reduce(nDefaultFaces, sumOp<label>());
}
if (nDefaultFaces > 0)
{
WarningIn("polyMesh::polyMesh(... construct from shapes...)")
<< "Found " << nFaces - defaultPatchStart
<< "Found " << nDefaultFaces
<< " undefined faces in mesh; adding to default patch." << endl;
// Check if there already exists a defaultFaces patch as last patch

View File

@ -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
@ -27,6 +27,7 @@ License
#include "PatchToolsCheck.C"
#include "PatchToolsEdgeOwner.C"
#include "PatchToolsGatherAndMerge.C"
#include "PatchToolsSearch.C"
#include "PatchToolsSortEdges.C"
#include "PatchToolsNormals.C"

View File

@ -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
@ -33,6 +33,7 @@ SourceFiles
PatchTools.C
PatchToolsCheck.C
PatchToolsEdgeOwner.C
PatchToolsGatherAndMerge.C
PatchToolsSearch.C
PatchToolsSortEdges.C
PatchToolsNormals.C
@ -236,6 +237,25 @@ public:
const labelList& meshFaces
);
//- Gather points and faces onto master and merge (geometrically) into
// single patch.
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
static void gatherAndMerge
(
const scalar mergeDist,
const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
Field<PointType>& mergedPoints,
List<Face>& mergedFaces,
labelList& pointMergeMap
);
};

View File

@ -0,0 +1,130 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
\\/ 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 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 "PatchTools.H"
#include "ListListOps.H"
#include "mergePoints.H"
#include "face.H"
#include "triFace.H"
#include "ListListOps.H"
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
template
<
class Face,
template<class> class FaceList,
class PointField,
class PointType
>
void Foam::PatchTools::gatherAndMerge
(
const scalar mergeDist,
const PrimitivePatch<Face, FaceList, PointField, PointType>& p,
Field<PointType>& mergedPoints,
List<Face>& mergedFaces,
labelList& pointMergeMap
)
{
// Collect points from all processors
labelList pointSizes;
{
List<Field<PointType> > gatheredPoints(Pstream::nProcs());
gatheredPoints[Pstream::myProcNo()] = p.localPoints();
Pstream::gatherList(gatheredPoints);
if (Pstream::master())
{
pointSizes = ListListOps::subSizes
(
gatheredPoints,
accessOp<Field<PointType> >()
);
mergedPoints = ListListOps::combine<Field<PointType> >
(
gatheredPoints,
accessOp<Field<PointType> >()
);
}
}
// Collect faces from all processors and renumber using sizes of
// gathered points
{
List<List<Face> > gatheredFaces(Pstream::nProcs());
gatheredFaces[Pstream::myProcNo()] = p.localFaces();
Pstream::gatherList(gatheredFaces);
if (Pstream::master())
{
mergedFaces = static_cast<const List<Face>&>
(
ListListOps::combineOffset<List<Face> >
(
gatheredFaces,
pointSizes,
accessOp<List<Face> >(),
offsetOp<Face>()
)
);
}
}
if (Pstream::master())
{
Field<PointType> newPoints;
labelList oldToNew;
bool hasMerged = mergePoints
(
mergedPoints,
mergeDist,
false, // verbosity
oldToNew,
newPoints
);
if (hasMerged)
{
// Store point mapping
pointMergeMap.transfer(oldToNew);
// Copy points
mergedPoints.transfer(newPoints);
// Relabel faces
List<Face>& faces = mergedFaces;
forAll(faces, faceI)
{
inplaceRenumber(pointMergeMap, faces[faceI]);
}
}
}
}
// ************************************************************************* //

View File

@ -340,10 +340,6 @@ $(snGradSchemes)/correctedSnGrad/correctedSnGrads.C
$(snGradSchemes)/limitedSnGrad/limitedSnGrads.C
$(snGradSchemes)/uncorrectedSnGrad/uncorrectedSnGrads.C
$(snGradSchemes)/orthogonalSnGrad/orthogonalSnGrads.C
/*
$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGradData.C
$(snGradSchemes)/quadraticFitSnGrad/quadraticFitSnGrads.C
*/
convectionSchemes = finiteVolume/convectionSchemes
$(convectionSchemes)/convectionScheme/convectionSchemes.C

View File

@ -1,159 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 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
quadraticFitSnGrad
Description
Simple central-difference snGrad scheme with quadratic fit correction from
a larger stencil.
SourceFiles
quadraticFitSnGrad.C
\*---------------------------------------------------------------------------*/
#ifndef quadraticFitSnGrad_H
#define quadraticFitSnGrad_H
#include "snGradScheme.H"
#include "quadraticFitSnGradData.H"
#include "extendedCellToFaceStencil.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace fv
{
/*---------------------------------------------------------------------------*\
Class quadraticFitSnGrad Declaration
\*---------------------------------------------------------------------------*/
template<class Type>
class quadraticFitSnGrad
:
public snGradScheme<Type>
{
// Private Data
//- weights for central stencil
const scalar centralWeight_;
// Private Member Functions
//- Disallow default bitwise assignment
void operator=(const quadraticFitSnGrad&);
public:
//- Runtime type information
TypeName("quadraticFit");
// Constructors
//- Construct from mesh and scheme data
quadraticFitSnGrad
(
const fvMesh& mesh,
const scalar centralWeight
)
:
snGradScheme<Type>(mesh),
centralWeight_(centralWeight)
{}
//- Construct from mesh and data stream
quadraticFitSnGrad(const fvMesh& mesh, Istream& is)
:
snGradScheme<Type>(mesh),
centralWeight_(readScalar(is))
{}
//- Destructor
virtual ~quadraticFitSnGrad() {}
// Member Functions
//- Return the interpolation weighting factors for the given field
virtual tmp<surfaceScalarField> deltaCoeffs
(
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().nonOrthDeltaCoeffs();
}
//- Return true if this scheme uses an explicit correction
virtual bool corrected() const
{
return true;
}
//- Return the explicit correction to the quadraticFitSnGrad
// for the given field
virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh> >
correction(const GeometricField<Type, fvPatchField, volMesh>& vf) const
{
const fvMesh& mesh = this->mesh();
const quadraticFitSnGradData& qfd = quadraticFitSnGradData::New
(
mesh,
centralWeight_
);
const extendedCellToFaceStencil& stencil = qfd.stencil();
const List<scalarList>& f = qfd.fit();
tmp<GeometricField<Type, fvsPatchField, surfaceMesh> > sft
= stencil.weightedSum(vf, f);
sft().dimensions() /= dimLength;
return sft;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace fv
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,321 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 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 "quadraticFitSnGradData.H"
#include "surfaceFields.H"
#include "volFields.H"
#include "SVD.H"
#include "syncTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(Foam::quadraticFitSnGradData, 0);
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::quadraticFitSnGradData::quadraticFitSnGradData
(
const fvMesh& mesh,
const scalar cWeight
)
:
MeshObject<fvMesh, quadraticFitSnGradData>(mesh),
centralWeight_(cWeight),
#ifdef SPHERICAL_GEOMETRY
dim_(2),
#else
dim_(mesh.nGeometricD()),
#endif
minSize_
(
dim_ == 1 ? 3 :
dim_ == 2 ? 6 :
dim_ == 3 ? 9 : 0
),
stencil_(mesh),
fit_(mesh.nInternalFaces())
{
if (debug)
{
Info<< "Contructing quadraticFitSnGradData" << endl;
}
// check input
if (centralWeight_ < 1 - SMALL)
{
FatalErrorIn("quadraticFitSnGradData::quadraticFitSnGradData")
<< "centralWeight requested = " << centralWeight_
<< " should not be less than one"
<< exit(FatalError);
}
if (minSize_ == 0)
{
FatalErrorIn("quadraticFitSnGradData")
<< " dimension must be 1,2 or 3, not" << dim_ << exit(FatalError);
}
// store the polynomial size for each face to write out
surfaceScalarField snGradPolySize
(
IOobject
(
"quadraticFitSnGradPolySize",
"constant",
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh,
dimensionedScalar("quadraticFitSnGradPolySize", dimless, scalar(0))
);
// Get the cell/face centres in stencil order.
// Centred face stencils no good for triangles of tets. Need bigger stencils
List<List<point> > stencilPoints(stencil_.stencil().size());
stencil_.collectData
(
mesh.C(),
stencilPoints
);
// find the fit coefficients for every face in the mesh
for (label faci = 0; faci < mesh.nInternalFaces(); faci++)
{
snGradPolySize[faci] = calcFit(stencilPoints[faci], faci);
}
if (debug)
{
snGradPolySize.write();
Info<< "quadraticFitSnGradData::quadraticFitSnGradData() :"
<< "Finished constructing polynomialFit data"
<< endl;
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::quadraticFitSnGradData::findFaceDirs
(
vector& idir, // value changed in return
vector& jdir, // value changed in return
vector& kdir, // value changed in return
const fvMesh& mesh,
const label faci
)
{
idir = mesh.Sf()[faci];
idir /= mag(idir);
#ifndef SPHERICAL_GEOMETRY
if (mesh.nGeometricD() <= 2) // find the normal direcion
{
if (mesh.geometricD()[0] == -1)
{
kdir = vector(1, 0, 0);
}
else if (mesh.geometricD()[1] == -1)
{
kdir = vector(0, 1, 0);
}
else
{
kdir = vector(0, 0, 1);
}
}
else // 3D so find a direction in the plane of the face
{
const face& f = mesh.faces()[faci];
kdir = mesh.points()[f[0]] - mesh.points()[f[1]];
}
#else
// Spherical geometry so kdir is the radial direction
kdir = mesh.Cf()[faci];
#endif
if (mesh.nGeometricD() == 3)
{
// Remove the idir component from kdir and normalise
kdir -= (idir & kdir)*idir;
scalar magk = mag(kdir);
if (magk < SMALL)
{
FatalErrorIn("findFaceDirs") << " calculated kdir = zero"
<< exit(FatalError);
}
else
{
kdir /= magk;
}
}
jdir = kdir ^ idir;
}
Foam::label Foam::quadraticFitSnGradData::calcFit
(
const List<point>& C,
const label faci
)
{
vector idir(1,0,0);
vector jdir(0,1,0);
vector kdir(0,0,1);
findFaceDirs(idir, jdir, kdir, mesh(), faci);
scalarList wts(C.size(), scalar(1));
wts[0] = centralWeight_;
wts[1] = centralWeight_;
point p0 = mesh().faceCentres()[faci];
scalar scale = 0;
// calculate the matrix of the polynomial components
scalarRectangularMatrix B(C.size(), minSize_, scalar(0));
forAll(C, ip)
{
const point& p = C[ip];
scalar px = (p - p0)&idir;
scalar py = (p - p0)&jdir;
#ifdef SPHERICAL_GEOMETRY
scalar pz = mag(p) - mag(p0);
#else
scalar pz = (p - p0)&kdir;
#endif
if (ip == 0) scale = max(max(mag(px), mag(py)), mag(pz));
px /= scale;
py /= scale;
pz /= scale;
label is = 0;
B[ip][is++] = wts[0]*wts[ip];
B[ip][is++] = wts[0]*wts[ip]*px;
B[ip][is++] = wts[ip]*sqr(px);
if (dim_ >= 2)
{
B[ip][is++] = wts[ip]*py;
B[ip][is++] = wts[ip]*px*py;
B[ip][is++] = wts[ip]*sqr(py);
}
if (dim_ == 3)
{
B[ip][is++] = wts[ip]*pz;
B[ip][is++] = wts[ip]*px*pz;
//B[ip][is++] = wts[ip]*py*pz;
B[ip][is++] = wts[ip]*sqr(pz);
}
}
// Set the fit
label stencilSize = C.size();
fit_[faci].setSize(stencilSize);
scalarList singVals(minSize_);
label nSVDzeros = 0;
const scalar deltaCoeff = deltaCoeffs()[faci];
bool goodFit = false;
for (int iIt = 0; iIt < 10 && !goodFit; iIt++)
{
SVD svd(B, SMALL);
scalar fit0 = wts[0]*wts[0]*svd.VSinvUt()[1][0]/scale;
scalar fit1 = wts[0]*wts[1]*svd.VSinvUt()[1][1]/scale;
goodFit =
fit0 < 0 && fit1 > 0
&& mag(fit0 + deltaCoeff) < 0.5*deltaCoeff
&& mag(fit1 - deltaCoeff) < 0.5*deltaCoeff;
if (goodFit)
{
fit_[faci][0] = fit0;
fit_[faci][1] = fit1;
for (label i = 2; i < stencilSize; i++)
{
fit_[faci][i] = wts[0]*wts[i]*svd.VSinvUt()[1][i]/scale;
}
singVals = svd.S();
nSVDzeros = svd.nZeros();
}
else // (not good fit so increase weight in the centre and for linear)
{
wts[0] *= 10;
wts[1] *= 10;
for (label i = 0; i < B.n(); i++)
{
B[i][0] *= 10;
B[i][1] *= 10;
}
for (label j = 0; j < B.m(); j++)
{
B[0][j] *= 10;
B[1][j] *= 10;
}
}
}
if (goodFit)
{
// remove the uncorrected snGradScheme coefficients
fit_[faci][0] += deltaCoeff;
fit_[faci][1] -= deltaCoeff;
}
else
{
Pout<< "quadratifFitSnGradData could not fit face " << faci
<< " fit_[faci][0] = " << fit_[faci][0]
<< " fit_[faci][1] = " << fit_[faci][1]
<< " deltaCoeff = " << deltaCoeff << endl;
fit_[faci] = 0;
}
return minSize_ - nSVDzeros;
}
bool Foam::quadraticFitSnGradData::movePoints()
{
notImplemented("quadraticFitSnGradData::movePoints()");
return true;
}
// ************************************************************************* //

View File

@ -1,129 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 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
quadraticFitSnGradData
Description
Data for the quadratic fit correction snGrad scheme
SourceFiles
quadraticFitSnGradData.C
\*---------------------------------------------------------------------------*/
#ifndef quadraticFitSnGradData_H
#define quadraticFitSnGradData_H
#include "MeshObject.H"
#include "fvMesh.H"
#include "extendedCellToFaceStencil.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
/*---------------------------------------------------------------------------*\
Class quadraticFitSnGradData Declaration
\*---------------------------------------------------------------------------*/
class quadraticFitSnGradData
:
public MeshObject<fvMesh, quadraticFitSnGradData>
{
// Private data
const scalar centralWeight_;
const label dim_;
//- minimum stencil size
const label minSize_;
//- Extended stencil addressing
extendedCellToFaceStencil stencil_;
//- For each cell in the mesh store the values which multiply the
// values of the stencil to obtain the gradient for each direction
List<scalarList> fit_;
// Private Member Functions
//- Find the normal direction and i, j and k directions for face faci
static void findFaceDirs
(
vector& idir, // value changed in return
vector& jdir, // value changed in return
vector& kdir, // value changed in return
const fvMesh& mesh,
const label faci
);
label calcFit(const List<point>&, const label faci);
public:
TypeName("quadraticFitSnGradData");
// Constructors
explicit quadraticFitSnGradData
(
const fvMesh& mesh,
const scalar cWeight
);
//- Destructor
virtual ~quadraticFitSnGradData()
{};
// Member functions
//- Return reference to the stencil
const extendedCellToFaceStencil& stencil() const
{
return stencil_;
}
//- Return reference to fit coefficients
const List<scalarList>& fit() const { return fit_; }
//- Delete the data when the mesh moves not implemented
virtual bool movePoints();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -1,43 +0,0 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
\\/ 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 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/>.
Description
Simple central-difference snGrad scheme with quadratic fit correction from
a larger stencil.
\*---------------------------------------------------------------------------*/
#include "quadraticFitSnGrad.H"
#include "fvMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace fv
{
makeSnGradScheme(quadraticFitSnGrad)
}
}
// ************************************************************************* //

View File

@ -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
@ -28,9 +28,8 @@ License
#include "dictionary.H"
#include "Time.H"
#include "IOmanip.H"
#include "ListListOps.H"
#include "mergePoints.H"
#include "volPointInterpolation.H"
#include "PatchTools.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -38,34 +37,6 @@ defineTypeNameAndDebug(Foam::sampledSurfaces, 0);
bool Foam::sampledSurfaces::verbose_ = false;
Foam::scalar Foam::sampledSurfaces::mergeTol_ = 1e-10;
namespace Foam
{
//- Used to offset faces in Pstream::combineOffset
template <>
class offsetOp<face>
{
public:
face operator()
(
const face& x,
const label offset
) const
{
face result(x.size());
forAll(x, xI)
{
result[xI] = x[xI] + offset;
}
return result;
}
};
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::sampledSurfaces::writeGeometry() const
@ -379,80 +350,18 @@ bool Foam::sampledSurfaces::update()
continue;
}
// Collect points from all processors
List<pointField> gatheredPoints(Pstream::nProcs());
gatheredPoints[Pstream::myProcNo()] = s.points();
Pstream::gatherList(gatheredPoints);
if (Pstream::master())
{
mergeList_[surfI].points = ListListOps::combine<pointField>
(
gatheredPoints,
accessOp<pointField>()
);
}
// Collect faces from all processors and renumber using sizes of
// gathered points
List<faceList> gatheredFaces(Pstream::nProcs());
gatheredFaces[Pstream::myProcNo()] = s.faces();
Pstream::gatherList(gatheredFaces);
if (Pstream::master())
{
mergeList_[surfI].faces = static_cast<const faceList&>
(
ListListOps::combineOffset<faceList>
(
gatheredFaces,
ListListOps::subSizes
(
gatheredPoints,
accessOp<pointField>()
),
accessOp<faceList>(),
offsetOp<face>()
)
);
}
pointField newPoints;
labelList oldToNew;
bool hasMerged = mergePoints
PatchTools::gatherAndMerge
(
mergeList_[surfI].points,
mergeDim,
false, // verbosity
oldToNew,
newPoints
primitivePatch
(
SubList<face>(s.faces(), s.faces().size()),
s.points()
),
mergeList_[surfI].points,
mergeList_[surfI].faces,
mergeList_[surfI].pointsMap
);
if (hasMerged)
{
// Store point mapping
mergeList_[surfI].pointsMap.transfer(oldToNew);
// Copy points
mergeList_[surfI].points.transfer(newPoints);
// Relabel faces
faceList& faces = mergeList_[surfI].faces;
forAll(faces, faceI)
{
inplaceRenumber(mergeList_[surfI].pointsMap, faces[faceI]);
}
if (Pstream::master() && debug)
{
Pout<< "For surface " << surfI << " merged from "
<< mergeList_[surfI].pointsMap.size() << " points down to "
<< mergeList_[surfI].points.size() << " points" << endl;
}
}
}
return updated;

View File

@ -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
@ -36,6 +36,7 @@ SourceFiles
#define labelledTri_H
#include "triFace.H"
#include "ListListOps.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -115,6 +116,30 @@ template<>
inline bool contiguous<labelledTri>() {return true;}
//- Hash specialization to offset faces in ListListOps::combineOffset
template<>
class offsetOp<labelledTri>
{
public:
inline labelledTri operator()
(
const labelledTri& x,
const label offset
) const
{
labelledTri result(x);
forAll(x, xI)
{
result[xI] = x[xI] + offset;
}
return result;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam