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

blockMesh: added experimental fast-merge algorithm

Browse Source 
The standard merge-algorithm is N^2 over the face-points and uses a
geometric proximity test for the merge.  These are both choices for
implementation simplicity and are rather inefficient for large meshes.
I have now implemented an experimental linear topological merge
algorithm which is VERY fast and effective for meshes of any size.
Currently it will merge internal faces on meshes of arbitrary complexity
but does not yet handle edge or face collapse needed for wedges and
other degenerate blocks.

The new fast-merge algorithm may be selected using the optional
"fastMerge" entry:

fastMerge yes;

and if not present the standard N^2 algorithm will be used.

Henry G. Weller
CFD Direct
This commit is contained in:
Henry Weller
2015-07-12 20:32:25 +01:00
parent b2968f3bda
commit 171c25ab76
11 changed files with 688 additions and 29 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
src/mesh/blockMesh/Make/files
View File

@ -24,5 +24,6 @@ blockMesh/blockMeshCreate.C
blockMesh/blockMeshTopology.C
blockMesh/blockMeshCheck.C
blockMesh/blockMeshMerge.C
blockMesh/blockMeshMergeFast.C
LIB = $(FOAM_LIBBIN)/libblockMesh

2
src/mesh/blockMesh/block/block.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-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License

13
src/mesh/blockMesh/block/block.H
View File

@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -126,13 +126,10 @@ public:
// Access
//- Return the block definition
inline const blockDescriptor& blockDef() const
{
return *this;
}
inline const blockDescriptor& blockDef() const;
//- Vertex label offset for a particular i,j,k position
label vtxLabel(label i, label j, label k) const;
inline label vtxLabel(label i, label j, label k) const;
//- Return the points for filling the block
const pointField& points() const;
@ -161,6 +158,10 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#include "blockI.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

13
src/mesh/blockMesh/block/blockCreate.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-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -28,17 +28,6 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::block::vtxLabel(label i, label j, label k) const
{
return
(
i
+ j * (meshDensity().x() + 1)
+ k * (meshDensity().x() + 1) * (meshDensity().y() + 1)
);
}
void Foam::block::createPoints() const
{
// set local variables for mesh specification

47
src/mesh/blockMesh/block/blockI.H Normal file
View File

@ -0,0 +1,47 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 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/>.
\*---------------------------------------------------------------------------*/
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
inline Foam::label Foam::block::vtxLabel(label i, label j, label k) const
{
return
(
i
+ j * (meshDensity().x() + 1)
+ k * (meshDensity().x() + 1) * (meshDensity().y() + 1)
);
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::blockDescriptor& Foam::block::blockDef() const
{
return *this;
}
// ************************************************************************* //

6
src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
View File

@ -142,6 +142,12 @@ public:
// Access
//- Return the number of cells in the i,j,k directions
const Vector<label>& density() const
{
return meshDensity_;
}
//- Reference to point field defining the block mesh
const pointField& blockPointField() const;

2
src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
View File

@ -68,7 +68,7 @@ void Foam::blockDescriptor::setEdge
// Set reference to the list of labels defining the block
const labelList& blockLabels = blockShape_;
// Set reference to global list of points
// Get list of points for this block
const pointField blockPoints = blockShape_.points(blockPointField_);
// Set the edge points/weights

17
src/mesh/blockMesh/blockMesh/blockMesh.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-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@ -24,11 +24,17 @@ License
\*---------------------------------------------------------------------------*/
#include "blockMesh.H"
#include "Switch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
bool Foam::blockMesh::blockMesh::verboseOutput(false);
namespace Foam
{
defineDebugSwitch(blockMesh, 0);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@ -38,7 +44,16 @@ Foam::blockMesh::blockMesh(const IOdictionary& dict, const word& regionName)
scaleFactor_(1.0),
topologyPtr_(createTopology(dict, regionName))
{
Switch fastMerge(dict.lookupOrDefault<Switch>("fastMerge", false));
if (fastMerge)
{
calcMergeInfoFast();
}
else
{
calcMergeInfo();
}
}

22
src/mesh/blockMesh/blockMesh/blockMesh.H
View File

@ -61,6 +61,8 @@ class blockMesh
public blockList
{
// Private data
//- Switch for verbose output
static bool verboseOutput;
//- Point field defining the block mesh (corners)
@ -134,6 +136,9 @@ class blockMesh
//- Determine the merge info and the final number of cells/points
void calcMergeInfo();
//- Determine the merge info and the final number of cells/points
void calcMergeInfoFast();
faceList createPatchFaces(const polyPatch& patchTopologyFaces) const;
void createPoints() const;
@ -146,6 +151,11 @@ class blockMesh
public:
// Static data members
ClassName("blockMesh");
// Constructors
//- Construct from IOdictionary
@ -164,8 +174,10 @@ public:
// these points have not been scaled by scaleFactor
const pointField& blockPointField() const;
//- Return the blockMesh topology as a polyMesh
const polyMesh& topology() const;
//- Return the curved edges
const curvedEdgeList& edges() const
{
return edges_;
@ -178,20 +190,18 @@ public:
// these points have been scaled by scaleFactor
const pointField& points() const;
//- Return cell shapes list
const cellShapeList& cells() const;
//- Return the patch face lists
const faceListList& patches() const;
//- Get patch information from the topology mesh
PtrList<dictionary> patchDicts() const;
//- Return patch names
wordList patchNames() const;
// wordList patchTypes() const;
//
// wordList patchPhysicalTypes() const;
//- Number of blocks with specified zones
label numZonedBlocks() const;
@ -204,6 +214,7 @@ public:
//- Enable/disable verbose information about the progress
static void verbose(const bool on=true);
// Write
//- Writes edges of blockMesh in OBJ format.
@ -220,4 +231,3 @@ public:
#endif
// ************************************************************************* //

2
src/mesh/blockMesh/blockMesh/blockMeshMerge.C
View File

@ -544,7 +544,7 @@ void Foam::blockMesh::calcMergeInfo()
}
// sort merge list to return new point label (in new shorter list)
// Sort merge list to return new point label (in new shorter list)
// given old point label
label newPointLabel = 0;

590
src/mesh/blockMesh/blockMesh/blockMeshMergeFast.C Normal file
View File

@ -0,0 +1,590 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2015 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 "blockMesh.H"
// * * * * * * * * * * * * * * * Private Functions * * * * * * * * * * * * * //
namespace Foam
{
// faces
// 6
// (
// 4(0 4 7 3) // x-min
// 4(1 2 6 5) // x-max
// 4(0 1 5 4) // y-min
// 4(3 7 6 2) // y-max
// 4(0 3 2 1) // z-min
// 4(4 5 6 7) // z-max
// );
// Face-edge directions
static const int faceEdgeDirs[6][4] =
{
{2, 1, -2, -1},
{1, 2, -1, -2},
{1, 2, -1, -2},
{2, 1, -2, -1},
{2, 1, -2, -1},
{1, 2, -1, -2}
};
// The face-face-rotation direction correspondence map
static Pair<int> faceFaceRotMap[6][6][4];
// Generate the face-face-rotation direction correspondence map
void genFaceFaceRotMap()
{
for(int facePi=0; facePi<6; facePi++)
{
for(int faceNi=0; faceNi<6; faceNi++)
{
for(int rot=0; rot<4; rot++)
{
Pair<int>& map = faceFaceRotMap[facePi][faceNi][rot];
for(int Pp=0; Pp<2; Pp++)
{
int Pdir = faceEdgeDirs[facePi][Pp];
int Np = (3 - Pp + rot)%4;
int Ndir = faceEdgeDirs[faceNi][Np];
map[Pdir-1] = -Ndir;
}
// Handle sign change due to match-face transpose
if (mag(map[0]) == 2 && map[0]*map[1] < 0)
{
map[0] = -map[0];
map[1] = -map[1];
}
}
}
}
}
// Return the direction map for the merge-faces
Pair<int> faceMap
(
const label facePi,
const face& faceP,
const label faceNi,
const face& faceN
)
{
// Search for the point on faceN corresponding to the 0-point on faceP
for(int rot=0; rot<4; rot++)
{
if (faceN[rot] == faceP[0])
{
return faceFaceRotMap[facePi][faceNi][rot];
}
}
FatalErrorIn
(
"faceMap(const label facePi, const face& faceP, "
"const label faceNi, const face& faceN)"
) << "Cannot find point correspondance for faces "
<< faceP << " and " << faceN
<< exit(FatalError);
return Pair<int>(0, 0);
}
// Set the block and face indices for all the merge faces
void setBlockFaceCorrespondence
(
const cellList& topoCells,
const faceList::subList& topoInternalFaces,
const labelList& topoFaceCell,
List<Pair<label> >& mergeBlock
)
{
forAll(topoInternalFaces, topoFacei)
{
label topoPi = topoFaceCell[topoFacei];
const labelList& topoPfaces = topoCells[topoPi];
bool foundFace = false;
label topoPfacei;
for
(
topoPfacei = 0;
topoPfacei < topoPfaces.size();
topoPfacei++
)
{
if (topoPfaces[topoPfacei] == topoFacei)
{
foundFace = true;
break;
}
}
if (!foundFace)
{
FatalErrorIn("setBlockFaceCorrespondence()")
<< "Cannot find merge face for block " << topoPi
<< exit(FatalError);
}
mergeBlock[topoFacei].first() = topoPi;
mergeBlock[topoFacei].second() = topoPfacei;
}
}
// Return the number of divisions in each direction for the face
Pair<label> faceNij(const label facei, const block& block)
{
Pair<label> fnij;
int i = facei/2;
if (i == 0)
{
fnij.first() = block.meshDensity().y() + 1;
fnij.second() = block.meshDensity().z() + 1;
}
else if (i == 1)
{
fnij.first() = block.meshDensity().x() + 1;
fnij.second() = block.meshDensity().z() + 1;
}
else if (i == 2)
{
fnij.first() = block.meshDensity().x() + 1;
fnij.second() = block.meshDensity().y() + 1;
}
return fnij;
}
// Sign the index corresponding to the map
inline label signIndex(const int map, const label i)
{
return map < 0 ? -i-1 : i;
}
// Reverse a signed index with the number of divisions
inline label unsignIndex(const label i, const label ni)
{
return i >= 0 ? i : ni + i + 1;
}
// Return the mapped index
inline label mapij(const int map, const label i, const label j)
{
return signIndex(map, mag(map) == 1 ? i : j);
}
// Return the face point index
inline label facePoint
(
const int facei,
const block& block,
const label i,
const label j
)
{
switch (facei)
{
case 0:
return block.vtxLabel(0, i, j);
case 1:
return block.vtxLabel(block.meshDensity().x(), i, j);
case 2:
return block.vtxLabel(i, 0, j);
case 3:
return block.vtxLabel(i, block.meshDensity().y(), j);
case 4:
return block.vtxLabel(i, j, 0);
case 5:
return block.vtxLabel(i, j, block.meshDensity().z());
default:
return -1;
}
}
// Return the neighbour face point from the signed indices
inline label facePointN
(
const block& block,
const label i,
const label j,
const label k
)
{
return block.vtxLabel
(
unsignIndex(i, block.meshDensity().x()),
unsignIndex(j, block.meshDensity().y()),
unsignIndex(k, block.meshDensity().z())
);
}
// Return the neighbour face point from the mapped indices
inline label facePointN
(
const int facei,
const block& block,
const label i,
const label j
)
{
switch (facei)
{
case 0:
return facePointN(block, 0, i, j);
case 1:
return facePointN(block, block.meshDensity().x(), i, j);
case 2:
return facePointN(block, i, 0, j);
case 3:
return facePointN(block, i, block.meshDensity().y(), j);
case 4:
return facePointN(block, i, j, 0);
case 5:
return facePointN(block, i, j, block.meshDensity().z());
default:
return -1;
}
}
// Return the neighbour face point using the map
inline label facePointN
(
const int facei,
const Pair<int>& fmap,
const block& block,
const label i,
const label j
)
{
return facePointN(facei, block, mapij(fmap[0], i, j), mapij(fmap[1], i, j));
}
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::blockMesh::calcMergeInfoFast()
{
// Generate the static face-face map
genFaceFaceRotMap();
const blockList& blocks = *this;
if (verboseOutput)
{
Info<< "Creating block offsets" << endl;
}
blockOffsets_.setSize(blocks.size());
nPoints_ = 0;
nCells_ = 0;
forAll(blocks, blockI)
{
blockOffsets_[blockI] = nPoints_;
nPoints_ += blocks[blockI].nPoints();
nCells_ += blocks[blockI].nCells();
}
if (verboseOutput)
{
Info<< "Creating merge list using the fast topological search"
<< flush;
}
// Size merge list and initialize to -1
mergeList_.setSize(nPoints_, -1);
// Block mesh topology
const pointField& topoPoints = topology().points();
const cellList& topoCells = topology().cells();
const faceList& topoFaces = topology().faces();
const labelList& topoFaceOwn = topology().faceOwner();
const labelList& topoFaceNei = topology().faceNeighbour();
// Topological merging is only necessary for internal block faces
// Note edge and face collapse may apply to boundary faces
// but is not yet supported in the "fast" algorithm
const faceList::subList topoInternalFaces
(
topoFaces,
topology().nInternalFaces()
);
List<Pair<label> > mergeBlockP(topoInternalFaces.size());
setBlockFaceCorrespondence
(
topoCells,
topoInternalFaces,
topoFaceOwn,
mergeBlockP
);
List<Pair<label> > mergeBlockN(topoInternalFaces.size());
setBlockFaceCorrespondence
(
topoCells,
topoInternalFaces,
topoFaceNei,
mergeBlockN
);
if (debug)
{
Info<< endl;
}
forAll(topoInternalFaces, topoFacei)
{
if (debug)
{
Info<< "Processing face " << topoFacei << endl;
}
label blockPi = mergeBlockP[topoFacei].first();
label blockPfacei = mergeBlockP[topoFacei].second();
label blockNi = mergeBlockN[topoFacei].first();
label blockNfacei = mergeBlockN[topoFacei].second();
Pair<int> fmap
(
faceMap
(
blockPfacei,
blocks[blockPi].blockShape().faces()[blockPfacei],
blockNfacei,
blocks[blockNi].blockShape().faces()[blockNfacei]
)
);
if (debug)
{
Info<< " Face map for faces "
<< blocks[blockPi].blockShape().faces()[blockPfacei] << " "
<< blocks[blockNi].blockShape().faces()[blockNfacei] << ": "
<< fmap << endl;
}
const pointField& blockPpoints = blocks[blockPi].points();
const pointField& blockNpoints = blocks[blockNi].points();
Pair<label> Pnij(faceNij(blockPfacei, blocks[blockPi]));
// Check block subdivision correspondence
{
Pair<label> Nnij(faceNij(blockNfacei, blocks[blockNi]));
Pair<label> NPnij;
NPnij[0] = Nnij[mag(fmap[0]) - 1];
NPnij[1] = Nnij[mag(fmap[1]) - 1];
if (Pnij != NPnij)
{
FatalErrorIn("blockMesh::calcMergeInfoFast()")
<< "Sub-division mismatch between face "
<< blockPfacei << " of block " << blockPi << Pnij
<< " and face "
<< blockNfacei << " of block " << blockNi << Nnij
<< exit(FatalError);
}
}
// Calculate a suitable test distance from the bounding box of the face.
// Note this is used only as a sanity check and for diagnostics in
// case there is a grading inconsistency.
const boundBox bb(topoCells[blockPi].points(topoFaces, topoPoints));
const scalar testSqrDist = magSqr(1e-6*bb.span());
// Accumulate the maximum merge distance for diagnostics
scalar maxSqrDist = 0;
for (label j=0; j<Pnij.second(); j++)
{
for (label i=0; i<Pnij.first(); i++)
{
label blockPpointi =
facePoint(blockPfacei, blocks[blockPi], i, j);
label blockNpointi =
facePointN(blockNfacei, fmap, blocks[blockNi], i, j);
scalar sqrDist
(
magSqr
(
blockPpoints[blockPpointi]
- blockNpoints[blockNpointi]
)
);
if (sqrDist > testSqrDist)
{
FatalErrorIn("blockMesh::calcMergeInfoFast()")
<< "Point merge failure between face "
<< blockPfacei << " of block " << blockPi
<< " and face "
<< blockNfacei << " of block " << blockNi
<< endl
<< " This may be due to inconsistent grading."
<< exit(FatalError);
}
maxSqrDist = max(maxSqrDist, sqrDist);
label Ppointi = blockPpointi + blockOffsets_[blockPi];
label Npointi = blockNpointi + blockOffsets_[blockNi];
label minPNi = min(Ppointi, Npointi);
if (mergeList_[Ppointi] != -1)
{
minPNi = min(minPNi, mergeList_[Ppointi]);
}
if (mergeList_[Npointi] != -1)
{
minPNi = min(minPNi, mergeList_[Npointi]);
}
mergeList_[Ppointi] = mergeList_[Npointi] = minPNi;
}
}
if (debug)
{
Info<< " Max distance between merge points: "
<< sqrt(maxSqrDist) << endl;
}
}
bool changedPointMerge = false;
label nPasses = 0;
do
{
changedPointMerge = false;
nPasses++;
forAll(topoInternalFaces, topoFacei)
{
label blockPi = mergeBlockP[topoFacei].first();
label blockPfacei = mergeBlockP[topoFacei].second();
label blockNi = mergeBlockN[topoFacei].first();
label blockNfacei = mergeBlockN[topoFacei].second();
Pair<int> fmap
(
faceMap
(
blockPfacei,
blocks[blockPi].blockShape().faces()[blockPfacei],
blockNfacei,
blocks[blockNi].blockShape().faces()[blockNfacei]
)
);
Pair<label> Pnij(faceNij(blockPfacei, blocks[blockPi]));
for (label j=0; j<Pnij.second(); j++)
{
for (label i=0; i<Pnij.first(); i++)
{
label blockPpointi =
facePoint(blockPfacei, blocks[blockPi], i, j);
label blockNpointi =
facePointN(blockNfacei, fmap, blocks[blockNi], i, j);
label Ppointi =
blockPpointi + blockOffsets_[blockPi];
label Npointi =
blockNpointi + blockOffsets_[blockNi];
if (mergeList_[Ppointi] != mergeList_[Npointi])
{
changedPointMerge = true;
mergeList_[Ppointi]
= mergeList_[Npointi]
= min(mergeList_[Ppointi], mergeList_[Npointi]);
}
}
}
}
if (verboseOutput)
{
Info<< "." << flush;
}
if (nPasses > 100)
{
FatalErrorIn("blockMesh::calcMergeInfoFast()")
<< "Point merging failed after 100 passes."
<< exit(FatalError);
}
} while (changedPointMerge);
// Sort merge list and count number of unique points
label nUniqPoints = 0;
forAll(mergeList_, pointi)
{
if (mergeList_[pointi] > pointi)
{
FatalErrorIn("blockMesh::calcMergeInfoFast()")
<< "Merge list contains point index out of range"
<< exit(FatalError);
}
if (mergeList_[pointi] == -1 || mergeList_[pointi] == pointi)
{
mergeList_[pointi] = nUniqPoints++;
}
else
{
mergeList_[pointi] = mergeList_[mergeList_[pointi]];
}
}
// Correct number of points in mesh
nPoints_ = nUniqPoints;
}
// ************************************************************************* //