research/OpenFOAM-2.4.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

blockMesh: Updated to latest version from OpenFOAM-dev:

Browse source 
Improved multi-grading so that the divisions are independent of section ordering.
Corrected multi-grading support for polyLine.
Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1809

Changed default location of blockMeshDict from constant/polyMesh to system
    For multi-region cases the default location of blockMeshDict is now system/<region name>

    If the blockMeshDict is not found in system then the constant directory
    is also checked providing backward-compatibility

Added experimental fast-merge algorithm

    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-08-06 09:08:58 +01:00
parent a6f90b1782
commit dfcf5d948e
16 changed files with 767 additions and 99 deletions
Show stats 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
// ************************************************************************* //

57
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,20 +28,9 @@ 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
// Set local variables for mesh specification
const label ni = meshDensity().x();
const label nj = meshDensity().y();
const label nk = meshDensity().z();
@ -57,12 +46,12 @@ void Foam::block::createPoints() const
const point& p011 = blockPoint(7);
// list of edge point and weighting factors
// List of edge point and weighting factors
const List< List<point> >& p = blockEdgePoints();
const scalarListList& w = blockEdgeWeights();
//
// generate vertices
// Generate vertices
//
vertices_.clear();
vertices_.setSize(nPoints());
@ -75,7 +64,7 @@ void Foam::block::createPoints() const
{
const label vertexNo = vtxLabel(i, j, k);
// points on edges
// Points on edges
vector edgex1 = p000 + (p100 - p000)*w[0][i];
vector edgex2 = p010 + (p110 - p010)*w[1][i];
vector edgex3 = p011 + (p111 - p011)*w[2][i];
@ -92,7 +81,7 @@ void Foam::block::createPoints() const
vector edgez4 = p010 + (p011 - p010)*w[11][k];
// calculate the importance factors for all edges
// Calculate the importance factors for all edges
// x-direction
scalar impx1 =
@ -190,7 +179,7 @@ void Foam::block::createPoints() const
impz4 /= magImpz;
// calculate the correction vectors
// Calculate the correction vectors
vector corx1 = impx1*(p[0][i] - edgex1);
vector corx2 = impx2*(p[1][i] - edgex2);
vector corx3 = impx3*(p[2][i] - edgex3);
@ -207,7 +196,7 @@ void Foam::block::createPoints() const
vector corz4 = impz4*(p[11][k] - edgez4);
// multiply by the importance factor
// Multiply by the importance factor
// x-direction
edgex1 *= impx1;
@ -228,7 +217,7 @@ void Foam::block::createPoints() const
edgez4 *= impz4;
// add the contributions
// Add the contributions
vertices_[vertexNo] =
(
edgex1 + edgex2 + edgex3 + edgex4
@ -255,7 +244,7 @@ void Foam::block::createCells() const
const label nk = meshDensity().z();
//
// generate cells
// Generate cells
//
cells_.clear();
cells_.setSize(nCells());
@ -292,7 +281,7 @@ void Foam::block::createBoundary() const
const label nk = meshDensity().z();
//
// generate boundaries on each side of the hex
// Generate boundaries on each side of the hex
//
boundaryPatches_.clear();
boundaryPatches_.setSize(6);
@ -311,7 +300,7 @@ void Foam::block::createBoundary() const
{
boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
// set the points
// Set the points
boundaryPatches_[wallLabel][wallCellLabel][0] =
vtxLabel(0, j, k);
boundaryPatches_[wallLabel][wallCellLabel][1] =
@ -321,7 +310,7 @@ void Foam::block::createBoundary() const
boundaryPatches_[wallLabel][wallCellLabel][3] =
vtxLabel(0, j + 1, k);
// update the counter
// Update the counter
wallCellLabel++;
}
}
@ -338,7 +327,7 @@ void Foam::block::createBoundary() const
{
boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
// set the points
// Set the points
boundaryPatches_[wallLabel][wallCellLabel][0] =
vtxLabel(ni, j, k);
boundaryPatches_[wallLabel][wallCellLabel][1] =
@ -348,7 +337,7 @@ void Foam::block::createBoundary() const
boundaryPatches_[wallLabel][wallCellLabel][3] =
vtxLabel(ni, j, k+1);
// update the counter
// Update the counter
wallCellLabel++;
}
}
@ -366,7 +355,7 @@ void Foam::block::createBoundary() const
{
boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
// set the points
// Set the points
boundaryPatches_[wallLabel][wallCellLabel][0] =
vtxLabel(i, 0, k);
boundaryPatches_[wallLabel][wallCellLabel][1] =
@ -376,7 +365,7 @@ void Foam::block::createBoundary() const
boundaryPatches_[wallLabel][wallCellLabel][3] =
vtxLabel(i, 0, k + 1);
// update the counter
// Update the counter
wallCellLabel++;
}
}
@ -393,7 +382,7 @@ void Foam::block::createBoundary() const
{
boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
// set the points
// Set the points
boundaryPatches_[wallLabel][wallCellLabel][0] =
vtxLabel(i, nj, k);
boundaryPatches_[wallLabel][wallCellLabel][1] =
@ -403,7 +392,7 @@ void Foam::block::createBoundary() const
boundaryPatches_[wallLabel][wallCellLabel][3] =
vtxLabel(i + 1, nj, k);
// update the counter
// Update the counter
wallCellLabel++;
}
}
@ -422,7 +411,7 @@ void Foam::block::createBoundary() const
{
boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
// set the points
// Set the points
boundaryPatches_[wallLabel][wallCellLabel][0] =
vtxLabel(i, j, 0);
boundaryPatches_[wallLabel][wallCellLabel][1] =
@ -432,7 +421,7 @@ void Foam::block::createBoundary() const
boundaryPatches_[wallLabel][wallCellLabel][3] =
vtxLabel(i + 1, j, 0);
// update the counter
// Update the counter
wallCellLabel++;
}
}
@ -449,7 +438,7 @@ void Foam::block::createBoundary() const
{
boundaryPatches_[wallLabel][wallCellLabel].setSize(4);
// set the points
// Set the points
boundaryPatches_[wallLabel][wallCellLabel][0] =
vtxLabel(i, j, nk);
boundaryPatches_[wallLabel][wallCellLabel][1] =
@ -459,7 +448,7 @@ void Foam::block::createBoundary() const
boundaryPatches_[wallLabel][wallCellLabel][3] =
vtxLabel(i, j + 1, nk);
// update the counter
// Update the counter
wallCellLabel++;
}
}

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;
}
// ************************************************************************* //

12
src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
View file

@ -59,7 +59,7 @@ Foam::blockDescriptor::blockDescriptor
<< exit(FatalError);
}
// create a list of edges
// Create a list of edges
makeBlockEdges();
}
@ -99,7 +99,7 @@ Foam::blockDescriptor::blockDescriptor
if (t.isPunctuation())
{
// new-style: read a list of 3 values
// New-style: read a list of 3 values
if (t.pToken() == token::BEGIN_LIST)
{
is >> meshDensity_;
@ -118,7 +118,7 @@ Foam::blockDescriptor::blockDescriptor
}
else
{
// old-style: read three labels
// Old-style: read three labels
is >> meshDensity_.x()
>> meshDensity_.y()
>> meshDensity_.z();
@ -134,7 +134,7 @@ Foam::blockDescriptor::blockDescriptor
if (expRatios.size() == 1)
{
// identical in x/y/z-directions
// Identical in x/y/z-directions
expand_ = expRatios[0];
}
else if (expRatios.size() == 3)
@ -171,7 +171,7 @@ Foam::blockDescriptor::blockDescriptor
<< exit(FatalError);
}
// create a list of edges
// Create a list of edges
makeBlockEdges();
}
@ -279,7 +279,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const blockDescriptor& bd)
const List<gradingDescriptors>& expand = bd.expand_;
// can we use a compact notation?
// Can we use a compact notation?
if
(
// x-direction

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

19
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))
{
calcMergeInfo();
Switch fastMerge(dict.lookupOrDefault<Switch>("fastMerge", false));
if (fastMerge)
{
calcMergeInfoFast();
}
else
{
calcMergeInfo();
}
}

26
src/mesh/blockMesh/blockMesh/blockMesh.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
@ -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,18 +136,26 @@ 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;
void createCells() const;
void createPatches() const;
//- as copy (not implemented)
//- As copy (not implemented)
blockMesh(const 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
// ************************************************************************* //

17
src/mesh/blockMesh/blockMesh/blockMeshMerge.C
View file

@ -86,11 +86,7 @@ void Foam::blockMesh::calcMergeInfo()
blockPfaceLabel++
)
{
if
(
blockFaces[blockPfaces[blockPfaceLabel]]
== blockFaces[blockFaceLabel]
)
if (blockPfaces[blockPfaceLabel] == blockFaceLabel)
{
foundFace = true;
break;
@ -345,7 +341,6 @@ void Foam::blockMesh::calcMergeInfo()
}
}
// FIXME? - there seems to be some logic missing here
label blockNfaceLabel;
for
@ -365,8 +360,6 @@ void Foam::blockMesh::calcMergeInfo()
}
}
// FIXME? - there seems to be some logic missing here
const labelListList& blockPfaceFaces =
blocks[blockPlabel].boundaryPatches()[blockPfaceLabel];
@ -416,7 +409,7 @@ void Foam::blockMesh::calcMergeInfo()
{
FatalErrorIn("blockMesh::calcMergeInfo()")
<< "Point merging failed after max number of passes."
<< abort(FatalError);
<< exit(FatalError);
}
}
while (changedPointMerge);
@ -551,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;
@ -560,7 +553,8 @@ void Foam::blockMesh::calcMergeInfo()
if (mergeList_[pointLabel] > pointLabel)
{
FatalErrorIn("blockMesh::calcMergeInfo()")
<< "ouch" << exit(FatalError);
<< "Merge list contains point index out of range"
<< exit(FatalError);
}
if
@ -579,7 +573,6 @@ void Foam::blockMesh::calcMergeInfo()
}
nPoints_ = newPointLabel;
}

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;
}
// ************************************************************************* //

6
src/mesh/blockMesh/curvedEdges/curvedEdge.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
@ -117,11 +117,11 @@ Foam::pointField Foam::curvedEdge::appendEndPoints
{
pointField allKnots(otherKnots.size() + 2);
// start/end knots
// Start/end knots
allKnots[0] = points[start];
allKnots[otherKnots.size() + 1] = points[end];
// intermediate knots
// Intermediate knots
forAll(otherKnots, knotI)
{
allKnots[knotI+1] = otherKnots[knotI];

34
src/mesh/blockMesh/curvedEdges/lineDivide.C
View file

@ -59,17 +59,38 @@ Foam::lineDivide::lineDivide
// Check that there are more divisions than sections
if (nDiv >= gd.size())
{
// Calculate distribution of divisions to be independent
// of the order of the sections
labelList secnDivs(gd.size());
label sumSecnDivs = 0;
label secnMaxDivs = 0;
forAll(gd, sectioni)
{
scalar nDivFrac = gd[sectioni].nDivFraction();
secnDivs[sectioni] = label(nDivFrac*nDiv + 0.5);
sumSecnDivs += secnDivs[sectioni];
// Find the section with the largest number of divisions
if (nDivFrac > gd[secnMaxDivs].nDivFraction())
{
secnMaxDivs = sectioni;
}
}
// Adjust the number of divisions on the section with the largest
// number so that the total is nDiv
if (sumSecnDivs != nDiv)
{
secnDivs[secnMaxDivs] += (nDiv - sumSecnDivs);
}
forAll(gd, sectioni)
{
scalar blockFrac = gd[sectioni].blockFraction();
scalar nDivFrac = gd[sectioni].nDivFraction();
scalar expRatio = gd[sectioni].expansionRatio();
label secnDiv = label(nDivFrac*nDiv + 0.5);
if (sectioni == gd.size() - 1)
{
secnDiv = nDiv - secnStart + 1;
}
label secnDiv = secnDivs[sectioni];
label secnEnd = secnStart + secnDiv;
// Calculate the spacing
@ -109,7 +130,6 @@ Foam::lineDivide::lineDivide
}
}
// Calculate the points
for (label i = 0; i <= nDiv; i++)
{

32
src/mesh/blockMesh/curvedEdges/polyLine.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
@ -23,7 +23,6 @@ License
\*---------------------------------------------------------------------------*/
#include "error.H"
#include "polyLine.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -41,7 +40,7 @@ void Foam::polyLine::calcParam()
param_[i] = param_[i-1] + mag(points_[i] - points_[i-1]);
}
// normalize on the interval 0-1
// Normalize on the interval 0-1
lineLength_ = param_.last();
for (label i=1; i < param_.size() - 1; i++)
{
@ -85,7 +84,7 @@ Foam::label Foam::polyLine::nSegments() const
Foam::label Foam::polyLine::localParameter(scalar& lambda) const
{
// check endpoints
// Check endpoints
if (lambda < SMALL)
{
lambda = 0;
@ -97,22 +96,20 @@ Foam::label Foam::polyLine::localParameter(scalar& lambda) const
return nSegments();
}
// search table of cumulative distances to find which line-segment
// we are on. Check the upper bound.
// Search table of cumulative distances to find which line-segment
// we are on.
// Check the upper bound.
label segmentI = 1;
while (param_[segmentI] < lambda)
{
segmentI++;
}
segmentI--; // we want the corresponding lower bound
segmentI--; // We want the corresponding lower bound
// the local parameter [0-1] on this line segment
// The local parameter [0-1] on this line segment
lambda =
(
( lambda - param_[segmentI] )
/ ( param_[segmentI+1] - param_[segmentI] )
);
(lambda - param_[segmentI])/(param_[segmentI+1] - param_[segmentI]);
return segmentI;
}
@ -120,7 +117,7 @@ Foam::label Foam::polyLine::localParameter(scalar& lambda) const
Foam::point Foam::polyLine::position(const scalar mu) const
{
// check endpoints
// Check end-points
if (mu < SMALL)
{
return points_.first();
@ -130,7 +127,6 @@ Foam::point Foam::polyLine::position(const scalar mu) const
return points_.last();
}
scalar lambda = mu;
label segment = localParameter(lambda);
return position(segment, lambda);
@ -143,7 +139,7 @@ Foam::point Foam::polyLine::position
const scalar mu
) const
{
// out-of-bounds
// Out-of-bounds
if (segment < 0)
{
return points_.first();
@ -156,7 +152,7 @@ Foam::point Foam::polyLine::position
const point& p0 = points()[segment];
const point& p1 = points()[segment+1];
// special cases - no calculation needed
// Special cases - no calculation needed
if (mu <= 0.0)
{
return p0;
@ -167,8 +163,8 @@ Foam::point Foam::polyLine::position
}
else
{
// linear interpolation
return points_[segment] + mu * (p1 - p0);
// Linear interpolation
return points_[segment] + mu*(p1 - p0);
}
}

2
src/mesh/blockMesh/gradingDescriptor/gradingDescriptors.C
View file

@ -47,7 +47,7 @@ Foam::gradingDescriptors Foam::gradingDescriptors::inv() const
forAll(ret, i)
{
ret[i] = operator[](i).inv();
ret[i] = operator[](ret.size() - i - 1).inv();
}
return ret;