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

ENH: blockMesh enhancements

Browse Source 
- support non-uniform scaling, prescaling and cartesian coordinate
  transformations.

  Eg, stretch in one direction and then rotate
  ```
  prescale  (1.5 1 1);
  transform
  {
      origin   (0 0 0);
      rotation
      {
          type  axisAngle;
          axis  (0 0 1);
          angle 45;
      }
  }
  ```

- support "transformed" versions of blockMesh vertices, topology.
  With the additional of transformations etc, a simplistic application
  of a single scale parameter is no longer sufficient.

    new:  blMesh.vertices(true);
    old:  blMesh.vertices() * blMesh.scaleFactor();

    new:  blMesh.topology(true);
    old:  N/A

- add individual edge access for blockDescriptor.
  Saves copying and duplicate calculations.

- handle '(block face)' specification for curved faces,
  which is ok for external block faces, but likely somewhat
  questionable if used for internal block faces.
This commit is contained in:
Mark Olesen
2021-07-18 22:01:42 +02:00
parent 4eeff16735
commit d31f351c15
17 changed files with 1029 additions and 460 deletions
Download Patch File Download Diff File Expand all files Collapse all files

5
applications/utilities/mesh/generation/blockMesh/blockMeshOBJ.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@ -16,7 +16,8 @@ Description
\*---------------------------------------------------------------------------*/
{
const polyMesh& topoMesh = blocks.topology();
refPtr<polyMesh> topoMeshPtr(blocks.topology(true));
const polyMesh& topoMesh = topoMeshPtr();
// Write mesh as edges
{

6
applications/utilities/mesh/generation/blockMesh/blockMeshVTK.H
View File

@ -5,7 +5,7 @@
\\ / A nd | www.openfoam.com
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2020 OpenCFD Ltd.
Copyright (C) 2020-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@ -19,7 +19,9 @@ Description
// Non-legacy and ASCII (mesh is small, want readable output)
const vtk::outputOptions writeOpts = vtk::formatType::INLINE_ASCII;
const polyMesh& topoMesh = blocks.topology();
refPtr<polyMesh> topoMeshPtr(blocks.topology(true));
const polyMesh& topoMesh = topoMeshPtr();
const vtk::vtuCells topoCells(topoMesh, writeOpts);
vtk::internalMeshWriter writer

148
src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -31,6 +31,67 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::blockDescriptor::assignGradings
(
const UList<gradingDescriptors>& ratios
)
{
bool ok = true;
switch (ratios.size())
{
case 0:
{
expand_.resize(12);
expand_ = gradingDescriptors();
break;
}
case 1:
{
// Identical in x/y/z-directions
expand_.resize(12);
expand_ = ratios[0];
break;
}
case 3:
{
expand_.resize(12);
// x-direction
expand_[0] = ratios[0];
expand_[1] = ratios[0];
expand_[2] = ratios[0];
expand_[3] = ratios[0];
// y-direction
expand_[4] = ratios[1];
expand_[5] = ratios[1];
expand_[6] = ratios[1];
expand_[7] = ratios[1];
// z-direction
expand_[8] = ratios[2];
expand_[9] = ratios[2];
expand_[10] = ratios[2];
expand_[11] = ratios[2];
break;
}
case 12:
{
expand_ = ratios;
break;
}
default:
{
ok = false;
break;
}
}
return ok;
}
void Foam::blockDescriptor::check(const Istream& is)
{
for (const label pointi : blockShape_)
@ -38,8 +99,8 @@ void Foam::blockDescriptor::check(const Istream& is)
if (pointi < 0 || pointi >= vertices_.size())
{
FatalIOErrorInFunction(is)
<< "Point label " << pointi
<< " out of range 0.." << vertices_.size() - 1
<< "Point label (" << pointi
<< ") out of range 0.." << vertices_.size() - 1
<< " in block " << *this
<< exit(FatalIOError);
}
@ -100,7 +161,7 @@ void Foam::blockDescriptor::check(const Istream& is)
}
void Foam::blockDescriptor::findCurvedFaces()
void Foam::blockDescriptor::findCurvedFaces(const label blockIndex)
{
const faceList shapeFaces(blockShape().faces());
@ -108,17 +169,21 @@ void Foam::blockDescriptor::findCurvedFaces()
{
forAll(blockFaces_, facei)
{
const face& f = blockFaces_[facei].vertices();
// Accept (<block> <face>) face description
if
(
face::sameVertices
(
blockFaces_[facei].vertices(),
shapeFaces[shapeFacei]
f.size() == 2
&& f[0] == blockIndex
&& f[1] == shapeFacei
)
|| face::sameVertices(f, shapeFaces[shapeFacei])
)
{
curvedFaces_[shapeFacei] = facei;
nCurvedFaces_++;
++nCurvedFaces_;
break;
}
}
@ -144,19 +209,15 @@ Foam::blockDescriptor::blockDescriptor
blockEdges_(edges),
blockFaces_(faces),
blockShape_(bshape),
expand_(expand),
expand_(),
zoneName_(zoneName),
curvedFaces_(-1),
nCurvedFaces_(0)
{
if (expand_.empty())
{
expand_.resize(12, gradingDescriptors());
}
else if (expand_.size() != 12)
if (!assignGradings(expand))
{
FatalErrorInFunction
<< "Unknown definition of expansion ratios"
<< "Unknown definition of expansion ratios: " << expand
<< exit(FatalError);
}
@ -167,7 +228,7 @@ Foam::blockDescriptor::blockDescriptor
Foam::blockDescriptor::blockDescriptor
(
const dictionary& dict,
const label index,
const label blockIndex,
const pointField& vertices,
const blockEdgeList& edges,
const blockFaceList& faces,
@ -178,7 +239,8 @@ Foam::blockDescriptor::blockDescriptor
vertices_(vertices),
blockEdges_(edges),
blockFaces_(faces),
expand_(12, gradingDescriptors()),
blockShape_(),
expand_(),
zoneName_(),
curvedFaces_(-1),
nCurvedFaces_(0)
@ -218,7 +280,7 @@ Foam::blockDescriptor::blockDescriptor
else
{
FatalIOErrorInFunction(is)
<< "incorrect token while reading n, expected '(', found "
<< "Incorrect token while reading n, expected '(', found "
<< t.info()
<< exit(FatalIOError);
}
@ -226,6 +288,10 @@ Foam::blockDescriptor::blockDescriptor
else
{
// Old-style: read three labels
IOWarningInFunction(is)
<< "Encountered old-style specification of mesh divisions"
<< endl;
is >> ijkMesh::sizes().x()
>> ijkMesh::sizes().y()
>> ijkMesh::sizes().z();
@ -237,47 +303,18 @@ Foam::blockDescriptor::blockDescriptor
is.putBack(t);
}
List<gradingDescriptors> expRatios(is);
List<gradingDescriptors> expand(is);
if (expRatios.size() == 1)
if (!assignGradings(expand))
{
// Identical in x/y/z-directions
expand_ = expRatios[0];
}
else if (expRatios.size() == 3)
{
// x-direction
expand_[0] = expRatios[0];
expand_[1] = expRatios[0];
expand_[2] = expRatios[0];
expand_[3] = expRatios[0];
// y-direction
expand_[4] = expRatios[1];
expand_[5] = expRatios[1];
expand_[6] = expRatios[1];
expand_[7] = expRatios[1];
// z-direction
expand_[8] = expRatios[2];
expand_[9] = expRatios[2];
expand_[10] = expRatios[2];
expand_[11] = expRatios[2];
}
else if (expRatios.size() == 12)
{
expand_ = expRatios;
}
else
{
FatalIOErrorInFunction(is)
<< "Unknown definition of expansion ratios: " << expRatios
<< exit(FatalIOError);
FatalErrorInFunction
<< "Unknown definition of expansion ratios: " << expand
<< exit(FatalError);
}
check(is);
findCurvedFaces();
findCurvedFaces(blockIndex);
}
@ -403,7 +440,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const blockDescriptor& bd)
}
os << ' ' << bd.density()
<< " simpleGrading (";
<< " grading (";
const List<gradingDescriptors>& expand = bd.grading();
@ -445,8 +482,7 @@ Foam::Ostream& Foam::operator<<(Ostream& os, const blockDescriptor& bd)
}
}
os << ")";
os << ')';
return os;
}

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

@ -109,22 +109,26 @@ class blockDescriptor
// Private Member Functions
//- Assign edge grading.
// \return false for unsupported specification
bool assignGradings(const UList<gradingDescriptors>& ratios);
//- Check block has outward-pointing faces
void check(const Istream& is);
//- Calculate the points and weights for the specified edge.
// Return the number of curved edges
label edgePointsWeights
//- Calculate the points and weights for the specified hex edge.
// Return the number of curved edges (0-1)
int calcEdgePointsWeights
(
pointField (&edgePoints)[12],
scalarList (&edgeWeights)[12],
const label edgei,
pointField& edgePoints,
scalarList& edgeWeights,
const label start,
const label end,
const label dim
const label nDiv,
const gradingDescriptors& expand
) const;
void findCurvedFaces();
void findCurvedFaces(const label blockIndex = -1);
public:
@ -156,7 +160,7 @@ public:
blockDescriptor
(
const dictionary& dict,
const label index,
const label blockIndex,
const pointField& vertices,
const blockEdgeList& edges,
const blockFaceList& faces,
@ -167,28 +171,28 @@ public:
// Member Functions
//- Reference to point field defining the block mesh
inline const pointField& vertices() const;
inline const pointField& vertices() const noexcept;
//- Return reference to the list of curved faces
inline const blockFaceList& blockFaces() const;
inline const blockFaceList& blockFaces() const noexcept;
//- Return the block shape
inline const cellShape& blockShape() const;
inline const cellShape& blockShape() const noexcept;
//- Return the mesh density (number of cells) in the i,j,k directions
inline const labelVector& density() const;
inline const labelVector& density() const noexcept;
//- Expansion ratios in all directions
const List<gradingDescriptors>& grading() const;
inline const List<gradingDescriptors>& grading() const noexcept;
//- Return the (optional) zone name
inline const word& zoneName() const;
inline const word& zoneName() const noexcept;
//- Curved-face labels for each block-face (-1 for flat faces)
inline const FixedList<label, 6>& curvedFaces() const;
inline const FixedList<label, 6>& curvedFaces() const noexcept;
//- Number of curved faces in this block
inline label nCurvedFaces() const;
inline label nCurvedFaces() const noexcept;
//- Return block point for local label i
inline const point& blockPoint(const label i) const;
@ -209,11 +213,32 @@ public:
inline bool edge(const label i, const label j, const label k) const;
//- Calculate the points and weights for all edges.
// Return the number of curved edges
label edgesPointsWeights
// \return the number of curved edges (0-12)
int edgesPointsWeights
(
pointField (&edgePoints)[12],
scalarList (&edgeWeights)[12]
pointField (&edgesPoints)[12],
scalarList (&edgesWeights)[12]
) const;
//- Calculate points and weights for specified edge,
//- using the specified number of divisions and grading
// \return True if the edge is curved
bool edgePointsWeights
(
const label edgei,
pointField& edgePoints,
scalarList& edgeWeights,
const label nDiv,
const gradingDescriptors& gd = gradingDescriptors()
) const;
//- Calculate points and weights for specified edge.
// \return True if the edge is curved
bool edgePointsWeights
(
const label edgei,
pointField& edgePoints,
scalarList& edgeWeights
) const;
//- Return true if point i,j,k addresses a block flat face or edge
@ -235,7 +260,7 @@ public:
static void write(Ostream&, const label blocki, const dictionary&);
//- Return info proxy
inline InfoProxy<blockDescriptor> info() const
InfoProxy<blockDescriptor> info() const
{
return *this;
}

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

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -30,24 +30,29 @@ License
#include "lineEdge.H"
#include "lineDivide.H"
// * * * * * * * * * * * * * * Local Data Members * * * * * * * * * * * * * //
// Warning.
// Ordering of edges needs to be the same as hex cell shape model
static const int hexEdge0[12] = { 0, 3, 7, 4, 0, 1, 5, 4, 0, 1, 2, 3 };
static const int hexEdge1[12] = { 1, 2, 6, 5, 3, 2, 6, 7, 4, 5, 6, 7 };
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::label Foam::blockDescriptor::edgePointsWeights
int Foam::blockDescriptor::calcEdgePointsWeights
(
pointField (&edgePoints)[12],
scalarList (&edgeWeights)[12],
const label edgei,
pointField& edgePoints,
scalarList& edgeWeights,
const label start,
const label end,
const label nDiv
const label nDiv,
const gradingDescriptors& expand
) const
{
// Set reference to the list of labels defining the block
const labelList& blockLabels = blockShape_;
// Get list of points for this block
const pointField blockPoints = blockShape_.points(vertices_);
// Set the edge points/weights
// The edge is a straight-line if it is not in the list of blockEdges
@ -55,55 +60,57 @@ Foam::label Foam::blockDescriptor::edgePointsWeights
{
const int cmp = cedge.compare(blockLabels[start], blockLabels[end]);
if (cmp)
{
if (cmp > 0)
{
// Curve has the same orientation
// Divide the line
const lineDivide divEdge(cedge, nDiv, expand_[edgei]);
const lineDivide divEdge(cedge, nDiv, expand);
edgePoints[edgei] = divEdge.points();
edgeWeights[edgei] = divEdge.lambdaDivisions();
edgePoints = divEdge.points();
edgeWeights = divEdge.lambdaDivisions();
return 1; // Found curved-edge: done
}
else
else if (cmp < 0)
{
// Curve has the opposite orientation
// Divide the line
const lineDivide divEdge(cedge, nDiv, expand_[edgei].inv());
const lineDivide divEdge(cedge, nDiv, expand.inv());
const pointField& p = divEdge.points();
const scalarList& d = divEdge.lambdaDivisions();
edgePoints[edgei].setSize(p.size());
edgeWeights[edgei].setSize(d.size());
edgePoints.resize(p.size());
edgeWeights.resize(d.size());
label pn = p.size() - 1;
// Copy in reverse order
const label pn = (p.size() - 1);
forAll(p, pi)
{
edgePoints[edgei][pi] = p[pn - pi];
edgeWeights[edgei][pi] = 1 - d[pn - pi];
}
edgePoints[pi] = p[pn - pi];
edgeWeights[pi] = 1 - d[pn - pi];
}
// Found curved-edge: done
return 1;
return 1; // Found curved-edge: done
}
}
// Not curved-edge: divide the edge as a straight line
// Get list of points for this block
const pointField blockPoints(blockShape_.points(vertices_));
lineDivide divEdge
(
blockEdges::lineEdge(blockPoints, start, end),
nDiv,
expand_[edgei]
expand
);
edgePoints[edgei] = divEdge.points();
edgeWeights[edgei] = divEdge.lambdaDivisions();
edgePoints = divEdge.points();
edgeWeights = divEdge.lambdaDivisions();
return 0;
}
@ -111,37 +118,89 @@ Foam::label Foam::blockDescriptor::edgePointsWeights
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
Foam::label Foam::blockDescriptor::edgesPointsWeights
int Foam::blockDescriptor::edgesPointsWeights
(
pointField (&edgePoints)[12],
scalarList (&edgeWeights)[12]
pointField (&edgesPoints)[12],
scalarList (&edgesWeights)[12]
) const
{
const label ni = sizes().x();
const label nj = sizes().y();
const label nk = sizes().z();
int nCurved = 0;
label nCurvedEdges = 0;
for (label edgei = 0; edgei < 12; ++edgei)
{
nCurved += calcEdgePointsWeights
(
edgesPoints[edgei],
edgesWeights[edgei],
hexEdge0[edgei],
hexEdge1[edgei],
// X-direction
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 0, 0, 1, ni);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 1, 3, 2, ni);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 2, 7, 6, ni);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 3, 4, 5, ni);
sizes()[edgei/4], // 12 edges -> 3 components (x,y,z)
expand_[edgei]
);
}
// Y-direction
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 4, 0, 3, nj);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 5, 1, 2, nj);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 6, 5, 6, nj);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 7, 4, 7, nj);
return nCurved;
}
// Z-direction
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 8, 0, 4, nk);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 9, 1, 5, nk);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 10, 2, 6, nk);
nCurvedEdges += edgePointsWeights(edgePoints, edgeWeights, 11, 3, 7, nk);
return nCurvedEdges;
bool Foam::blockDescriptor::edgePointsWeights
(
const label edgei,
pointField& edgePoints,
scalarList& edgeWeights,
const label nDiv,
const gradingDescriptors& gd
) const
{
if (edgei < 0 || edgei >= 12)
{
FatalErrorInFunction
<< "Edge label " << edgei
<< " out of range 0..11"
<< exit(FatalError);
}
const int nCurved = calcEdgePointsWeights
(
edgePoints,
edgeWeights,
hexEdge0[edgei],
hexEdge1[edgei],
nDiv,
gd
);
return nCurved;
}
bool Foam::blockDescriptor::edgePointsWeights
(
const label edgei,
pointField& edgePoints,
scalarList& edgeWeights
) const
{
if (edgei < 0 || edgei >= 12)
{
FatalErrorInFunction
<< "Edge label " << edgei
<< " out of range 0..11"
<< exit(FatalError);
}
const int nCurved = calcEdgePointsWeights
(
edgePoints,
edgeWeights,
hexEdge0[edgei],
hexEdge1[edgei],
sizes()[edgei/4], // 12 edges -> 3 components (x,y,z)
expand_[edgei]
);
return nCurved;
}

24
src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2019 OpenCFD Ltd.
Copyright (C) 2019-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -28,51 +28,55 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::pointField& Foam::blockDescriptor::vertices() const
inline const Foam::pointField&
Foam::blockDescriptor::vertices() const noexcept
{
return vertices_;
}
inline const Foam::blockFaceList& Foam::blockDescriptor::blockFaces() const
inline const Foam::blockFaceList&
Foam::blockDescriptor::blockFaces() const noexcept
{
return blockFaces_;
}
inline const Foam::cellShape& Foam::blockDescriptor::blockShape() const
inline const Foam::cellShape&
Foam::blockDescriptor::blockShape() const noexcept
{
return blockShape_;
}
inline const Foam::labelVector& Foam::blockDescriptor::density() const
inline const Foam::labelVector&
Foam::blockDescriptor::density() const noexcept
{
return ijkMesh::sizes();
}
inline const Foam::List<Foam::gradingDescriptors>&
Foam::blockDescriptor::grading() const
Foam::blockDescriptor::grading() const noexcept
{
return expand_;
}
inline const Foam::word& Foam::blockDescriptor::zoneName() const
inline const Foam::word& Foam::blockDescriptor::zoneName() const noexcept
{
return zoneName_;
}
inline const Foam::FixedList<Foam::label, 6>&
Foam::blockDescriptor::curvedFaces() const
Foam::blockDescriptor::curvedFaces() const noexcept
{
return curvedFaces_;
}
inline Foam::label Foam::blockDescriptor::nCurvedFaces() const
inline Foam::label Foam::blockDescriptor::nCurvedFaces() const noexcept
{
return nCurvedFaces_;
}
@ -154,8 +158,10 @@ inline bool Foam::blockDescriptor::flatFaceOrEdge
{
if (i == 0 && curvedFaces_[0] == -1) return true;
if (i == sizes().x() && curvedFaces_[1] == -1) return true;
if (j == 0 && curvedFaces_[2] == -1) return true;
if (j == sizes().y() && curvedFaces_[3] == -1) return true;
if (k == 0 && curvedFaces_[4] == -1) return true;
if (k == sizes().z() && curvedFaces_[5] == -1) return true;

6
src/mesh/blockMesh/blockFaces/blockFace/blockFace.H
View File

@ -45,7 +45,7 @@ SourceFiles
namespace Foam
{
// Forward declarations
// Forward Declarations
class blockDescriptor;
class blockFace;
@ -91,7 +91,7 @@ public:
// Constructors
//- Construct from face vertices
blockFace(const face& vertices);
explicit blockFace(const face& vertices);
//- Construct from Istream
blockFace
@ -163,7 +163,7 @@ public:
void write(Ostream&, const dictionary&) const;
// Ostream operator
// Ostream Operator
friend Ostream& operator<<(Ostream&, const blockFace&);
};

318
src/mesh/blockMesh/blockMesh/blockMesh.C
View File

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2018-2020 OpenCFD Ltd.
Copyright (C) 2018-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -27,6 +27,7 @@ License
\*---------------------------------------------------------------------------*/
#include "blockMesh.H"
#include "transform.H"
#include "Time.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@ -47,6 +48,182 @@ Foam::blockMesh::strategyNames_
});
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
// Process dictionary entry with single scalar or vector quantity
// - return 0 if scaling is not needed. Eg, Not found or is unity
// - return 1 for uniform scaling
// - return 3 for non-uniform scaling
int readScaling(const entry *eptr, vector& scale)
{
int nCmpt = 0;
if (!eptr)
{
return nCmpt;
}
ITstream& is = eptr->stream();
if (is.peek().isNumber())
{
// scalar value
scalar val;
is >> val;
if ((val > 0) && !equal(val, 1))
{
// Uniform scaling
nCmpt = 1;
scale = vector::uniform(val);
}
}
else
{
// vector value
is >> scale;
bool nonUnity = false;
for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
{
if (scale[cmpt] <= 0)
{
scale[cmpt] = 1;
}
else if (!equal(scale[cmpt], 1))
{
nonUnity = true;
}
}
if (nonUnity)
{
if (equal(scale.x(), scale.y()) && equal(scale.x(), scale.z()))
{
// Uniform scaling
nCmpt = 1;
}
else
{
nCmpt = 3;
}
}
}
eptr->checkITstream(is);
return nCmpt;
}
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
bool Foam::blockMesh::readPointTransforms(const dictionary& dict)
{
transformType_ = transformTypes::NO_TRANSFORM;
// Optional cartesian coordinate system transform, since JUL-2021
if (const dictionary* dictptr = dict.findDict("transform"))
{
transform_ = coordSystem::cartesian(*dictptr);
constexpr scalar tol = VSMALL;
// Check for non-zero origin
const vector& o = transform_.origin();
if
(
(mag(o.x()) > tol)
|| (mag(o.y()) > tol)
|| (mag(o.z()) > tol)
)
{
transformType_ |= transformTypes::TRANSLATION;
}
// Check for non-identity rotation
const tensor& r = transform_.R();
if
(
(mag(r.xx() - 1) > tol)
|| (mag(r.xy()) > tol)
|| (mag(r.xz()) > tol)
|| (mag(r.yx()) > tol)
|| (mag(r.yy() - 1) > tol)
|| (mag(r.yz()) > tol)
|| (mag(r.zx()) > tol)
|| (mag(r.zy()) > tol)
|| (mag(r.zz() - 1) > tol)
)
{
transformType_ |= transformTypes::ROTATION;
}
}
else
{
transform_.clear();
}
// Optional 'prescale' factor.
{
prescaling_ = vector::uniform(1);
const int scaleType = readScaling
(
dict.findEntry("prescale", keyType::LITERAL),
prescaling_
);
if (scaleType == 1)
{
transformType_ |= transformTypes::PRESCALING;
}
else if (scaleType == 3)
{
transformType_ |= transformTypes::PRESCALING3;
}
}
// Optional 'scale' factor. Was 'convertToMeters' until OCT-2008
{
scaling_ = vector::uniform(1);
const int scaleType = readScaling
(
// Mark as changed from 2010 onwards
dict.findCompat
(
"scale",
{{"convertToMeters", 1012}},
keyType::LITERAL
),
scaling_
);
if (scaleType == 1)
{
transformType_ |= transformTypes::SCALING;
}
else if (scaleType == 3)
{
transformType_ |= transformTypes::SCALING3;
}
}
return bool(transformType_);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::blockMesh::blockMesh
@ -62,6 +239,8 @@ Foam::blockMesh::blockMesh
(
meshDict_.getOrDefault("checkFaceCorrespondence", true)
),
mergeStrategy_(strategy),
transformType_(transformTypes::NO_TRANSFORM),
geometry_
(
IOobject
@ -78,25 +257,27 @@ Foam::blockMesh::blockMesh
: dictionary(),
true
),
scaleFactor_(1),
blockVertices_
(
meshDict_.lookup("vertices"),
blockVertex::iNew(meshDict_, geometry_)
),
vertices_(Foam::vertices(blockVertices_)),
prescaling_(vector::uniform(1)),
scaling_(vector::uniform(1)),
transform_(),
topologyPtr_(createTopology(meshDict_, regionName))
{
// Command-line option has precedence over dictionary setting
if (strategy == mergeStrategy::DEFAULT_MERGE)
if (mergeStrategy_ == mergeStrategy::DEFAULT_MERGE)
{
strategyNames_.readIfPresent("mergeType", meshDict_, strategy);
strategyNames_.readIfPresent("mergeType", meshDict_, mergeStrategy_);
// Warn about fairly obscure old "fastMerge" option?
}
if (strategy == mergeStrategy::MERGE_POINTS)
if (mergeStrategy_ == mergeStrategy::MERGE_POINTS)
{
// MERGE_POINTS
calcGeometricalMerge();
@ -131,35 +312,38 @@ bool Foam::blockMesh::verbose(const bool on) noexcept
}
const Foam::pointField& Foam::blockMesh::vertices() const
const Foam::pointField& Foam::blockMesh::vertices() const noexcept
{
return vertices_;
}
const Foam::polyMesh& Foam::blockMesh::topology() const
Foam::tmp<Foam::pointField>
Foam::blockMesh::vertices(bool applyTransform) const
{
if (!topologyPtr_)
if (applyTransform && hasPointTransforms())
{
FatalErrorInFunction
<< "topologyPtr_ not allocated"
<< exit(FatalError);
auto tpts = tmp<pointField>::New(vertices_);
inplacePointTransforms(tpts.ref());
return tpts;
}
return *topologyPtr_;
return vertices_;
}
Foam::PtrList<Foam::dictionary> Foam::blockMesh::patchDicts() const
{
const polyPatchList& patchTopologies = topology().boundaryMesh();
const polyPatchList& topoPatches = topology().boundaryMesh();
PtrList<dictionary> patchDicts(patchTopologies.size());
PtrList<dictionary> patchDicts(topoPatches.size());
forAll(patchTopologies, patchi)
forAll(topoPatches, patchi)
{
OStringStream os;
patchTopologies[patchi].write(os);
topoPatches[patchi].write(os);
IStringStream is(os.str());
patchDicts.set(patchi, new dictionary(is));
}
@ -167,12 +351,6 @@ Foam::PtrList<Foam::dictionary> Foam::blockMesh::patchDicts() const
}
Foam::scalar Foam::blockMesh::scaleFactor() const
{
return scaleFactor_;
}
const Foam::pointField& Foam::blockMesh::points() const
{
if (points_.empty())
@ -242,4 +420,100 @@ Foam::label Foam::blockMesh::numZonedBlocks() const
}
bool Foam::blockMesh::hasPointTransforms() const noexcept
{
return bool(transformType_);
}
bool Foam::blockMesh::inplacePointTransforms(pointField& pts) const
{
if (!transformType_)
{
return false;
}
if (transformType_ & transformTypes::PRESCALING)
{
for (point& p : pts)
{
p *= prescaling_.x();
}
}
else if (transformType_ & transformTypes::PRESCALING3)
{
for (point& p : pts)
{
p = cmptMultiply(p, prescaling_);
}
}
if (transformType_ & transformTypes::ROTATION)
{
const tensor rot(transform_.R());
if (transformType_ & transformTypes::TRANSLATION)
{
const point origin(transform_.origin());
for (point& p : pts)
{
p = Foam::transform(rot, p) + origin;
}
}
else
{
for (point& p : pts)
{
p = Foam::transform(rot, p);
}
}
}
else if (transformType_ & transformTypes::TRANSLATION)
{
const point origin(transform_.origin());
for (point& p : pts)
{
p += origin;
}
}
if (transformType_ & transformTypes::SCALING)
{
for (point& p : pts)
{
p *= scaling_.x();
}
}
else if (transformType_ & transformTypes::SCALING3)
{
for (point& p : pts)
{
p = cmptMultiply(p, scaling_);
}
}
return true;
}
Foam::tmp<Foam::pointField>
Foam::blockMesh::globalPosition(const pointField& localPoints) const
{
if (hasPointTransforms())
{
auto tpts = tmp<pointField>::New(localPoints);
inplacePointTransforms(tpts.ref());
return tpts;
}
else
{
return localPoints;
}
}
// ************************************************************************* //

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

@ -6,7 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2017 OpenFOAM Foundation
Copyright (C) 2018-2020 OpenCFD Ltd.
Copyright (C) 2018-2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -33,7 +33,9 @@ Description
Dictionary controls
\table
Property | Description | Required | Default
scale | Point scaling | no | 1.0
prescale | Point scaling before transform | no | 1.0
scale | Point scaling after transform | no | 1.0
transform | Point transform (origin, rotation) | no |
vertices | | yes |
blocks | | yes |
edges | | no |
@ -48,6 +50,9 @@ Description
\endtable
Note
The \c prescale and \c scale can be a single scalar or a vector of
values.
The vertices, cells and patches for filling the blocks are demand-driven.
SourceFiles
@ -65,6 +70,7 @@ SourceFiles
#include "Enum.H"
#include "block.H"
#include "PtrList.H"
#include "cartesianCS.H"
#include "searchableSurfaces.H"
#include "polyMesh.H"
#include "IOdictionary.H"
@ -112,6 +118,21 @@ private:
static const Enum<mergeStrategy> strategyNames_;
// Data Types
//- Point transformations. Internal usage
enum transformTypes : unsigned
{
NO_TRANSFORM = 0, //!< No transformations
ROTATION = 0x1, //!< Local coordinate system rotation
TRANSLATION = 0x2, //!< Local coordinate system translation
PRESCALING = 0x4, //!< Uniform scale before rotate/translate
PRESCALING3 = 0x8, //!< Non-uniform scale before rotate/translate
SCALING = 0x10, //!< Uniform scale after rotate/translate
SCALING3 = 0x20 //!< Non-uniform scale after rotate/translate
};
// Private Data
//- Reference to mesh dictionary
@ -120,15 +141,18 @@ private:
//- Output verbosity
bool verbose_;
//- Switch checking face consistency (defaults to true)
//- Check face consistency (defaults to true)
bool checkFaceCorrespondence_;
//- The mesh merging strategy
mergeStrategy mergeStrategy_;
//- Types of point transformations requested
unsigned transformType_;
//- Optional searchable geometry to project face-points to
searchableSurfaces geometry_;
//- The scaling factor to convert to metres
scalar scaleFactor_;
//- The list of block vertices
blockVertexList blockVertices_;
@ -141,6 +165,15 @@ private:
//- The list of curved faces
blockFaceList faces_;
//- The scaling factor before rotate/translate
vector prescaling_;
//- The scaling factor after rotate/translate
vector scaling_;
//- Local coordinate system transformation
coordSystem::cartesian transform_;
//- The blocks themselves (the topology) as a polyMesh
autoPtr<polyMesh> topologyPtr_;
@ -150,12 +183,12 @@ private:
//- The sum of all cells in each block
label nCells_;
//- The point offset added to each block
labelList blockOffsets_;
//- The merge points information
labelList mergeList_;
//- The point offset added to each block. Offsets into mergeList_
labelList blockOffsets_;
mutable pointField points_;
mutable cellShapeList cells_;
@ -165,14 +198,9 @@ private:
// Private Member Functions
template<class Source>
void checkPatchLabels
(
const Source& source,
const word& patchName,
const pointField& points,
faceList& patchShapes
) const;
//- Get scaling and/or coordinate transforms
// \return True if scaling and/or transformations are needed
bool readPointTransforms(const dictionary& dict);
void readPatches
(
@ -191,7 +219,8 @@ private:
PtrList<dictionary>& patchDicts
);
void createCellShapes(cellShapeList& tmpBlockCells);
//- Topology blocks as cell shapes
cellShapeList getBlockShapes() const;
autoPtr<polyMesh> createTopology
(
@ -253,47 +282,77 @@ public:
// Member Functions
// Access
// General Access, Description
//- Access to input dictionary
const dictionary& meshDict() const
const dictionary& meshDict() const noexcept
{
return meshDict_;
}
//- Optional searchable geometry to project face-points to
const searchableSurfaces& geometry() const
const searchableSurfaces& geometry() const noexcept
{
return geometry_;
}
//- The curved edges
const blockEdgeList& edges() const noexcept
{
return edges_;
}
//- The curved faces
const blockFaceList& faces() const noexcept
{
return faces_;
}
//- True if the blockMesh topology exists
bool valid() const noexcept;
//- Reference to point field defining the blockMesh.
// These points are \b not scaled by scaleFactor
const pointField& vertices() const;
//- Return the patch names
wordList patchNames() const;
//- Return the blockMesh topology as a polyMesh
//- Number of blocks with specified zones
label numZonedBlocks() const;
// Point Transformations
//- True if scaling and/or transformations are needed
bool hasPointTransforms() const noexcept;
//- Apply coordinate transforms and scaling
bool inplacePointTransforms(pointField& pts) const;
//- Apply coordinate transforms and scaling
tmp<pointField> globalPosition(const pointField& localPoints) const;
// Block Topology
//- Reference to point field defining the blocks,
//- these points are \b unscaled and \b non-transformed
const pointField& vertices() const noexcept;
//- Point field defining the blocks,
//- optionally transformed and scaled
tmp<pointField> vertices(bool applyTransform) const;
//- The blockMesh topology as a polyMesh
//- \b unscaled and \b non-transformed
const polyMesh& topology() const;
//- Return the curved edges
const blockEdgeList& edges() const
{
return edges_;
}
//- The blockMesh topology as a polyMesh
//- optionally transformed and scaled
refPtr<polyMesh> topology(bool applyTransform) const;
//- Return the curved faces
const blockFaceList& faces() const
{
return faces_;
}
//- The scaling factor used to convert to metres
scalar scaleFactor() const;
// Detailed Mesh
//- The points for the entire mesh.
// These points \b are scaled by scaleFactor
//- These points \b are scaled and transformed
const pointField& points() const;
//- Return cell shapes list
@ -302,15 +361,9 @@ public:
//- Return the patch face lists
const faceListList& patches() const;
//- Get patch information from the topology mesh
//- Patch information from the topology mesh
PtrList<dictionary> patchDicts() const;
//- Return patch names
wordList patchNames() const;
//- Number of blocks with specified zones
label numZonedBlocks() const;
// Verbosity
@ -326,6 +379,16 @@ public:
//- Create polyMesh, with cell zones
autoPtr<polyMesh> mesh(const IOobject& io) const;
// Housekeeping
//- Old (v2106 and earlier) uniform scaling factor
// \deprecated use inplacePointTransforms or globalPosition instead
scalar scaleFactor() const
{
return scaling_.x();
}
};

26
src/mesh/blockMesh/blockMesh/blockMeshCheck.C
View File

@ -6,6 +6,7 @@
\\/ M anipulation |
-------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2021 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -105,20 +106,37 @@ void Foam::blockMesh::check(const polyMesh& bm, const dictionary& dict) const
}
// Check curved-face/block-face correspondence
forAll(faces_, cfi)
for (const blockFace& cface : faces_)
{
const face& cf = cface.vertices();
bool found = false;
forAll(faces, fi)
if (cf.size() == 2)
{
found = faces_[cfi].compare(faces[fi]) != 0;
const label bi = cf[0];
const label fi = cf[1];
found =
(
bi >= 0 && bi < blocks.size()
&& fi >= 0 && fi < blocks[bi].blockShape().nFaces()
);
}
if (!found)
{
for (const face& bf : faces)
{
found = cface.compare(bf) != 0;
if (found) break;
}
}
if (!found)
{
Info<< " Curved face ";
faces_[cfi].write(Info, dict);
cface.write(Info, dict);
Info<< " does not correspond to a block face." << endl;
ok = false;
}

203
src/mesh/blockMesh/blockMesh/blockMeshCreate.C
View File

@ -36,64 +36,85 @@ void Foam::blockMesh::createPoints() const
{
const blockList& blocks = *this;
const vector scaleTot
(
prescaling_.x() * scaling_.x(),
prescaling_.y() * scaling_.y(),
prescaling_.z() * scaling_.z()
);
if (verbose_)
{
Info<< "Creating points with scale " << scaleFactor_ << endl;
Info<< "Creating points with scale " << scaleTot << endl;
}
points_.setSize(nPoints_);
points_.resize(nPoints_);
forAll(blocks, blocki)
{
const pointField& blockPoints = blocks[blocki].points();
const labelSubList pointAddr
(
mergeList_,
blockPoints.size(),
blockOffsets_[blocki]
);
UIndirectList<point>(points_, pointAddr) = blockPoints;
if (verbose_)
{
const label nx = blocks[blocki].density().x();
const label ny = blocks[blocki].density().y();
const label nz = blocks[blocki].density().z();
Info<< " Block " << blocki << " cell size :" << nl;
label v0 = blocks[blocki].pointLabel(0, 0, 0);
label vi1 = blocks[blocki].pointLabel(1, 0, 0);
scalar diStart = mag(blockPoints[vi1] - blockPoints[v0]);
const label v0 = blocks[blocki].pointLabel(0, 0, 0);
label vinM1 = blocks[blocki].pointLabel(nx-1, 0, 0);
label vin = blocks[blocki].pointLabel(nx, 0, 0);
scalar diFinal = mag(blockPoints[vin] - blockPoints[vinM1]);
label vj1 = blocks[blocki].pointLabel(0, 1, 0);
scalar djStart = mag(blockPoints[vj1] - blockPoints[v0]);
label vjnM1 = blocks[blocki].pointLabel(0, ny-1, 0);
label vjn = blocks[blocki].pointLabel(0, ny, 0);
scalar djFinal = mag(blockPoints[vjn] - blockPoints[vjnM1]);
label vk1 = blocks[blocki].pointLabel(0, 0, 1);
scalar dkStart = mag(blockPoints[vk1] - blockPoints[v0]);
label vknM1 = blocks[blocki].pointLabel(0, 0, nz-1);
label vkn = blocks[blocki].pointLabel(0, 0, nz);
scalar dkFinal = mag(blockPoints[vkn] - blockPoints[vknM1]);
Info<< " Block " << blocki << " cell size :" << nl
<< " i : " << scaleFactor_*diStart << " .. "
<< scaleFactor_*diFinal << nl
<< " j : " << scaleFactor_*djStart << " .. "
<< scaleFactor_*djFinal << nl
<< " k : " << scaleFactor_*dkStart << " .. "
<< scaleFactor_*dkFinal << nl
<< endl;
}
forAll(blockPoints, blockPointi)
{
points_
[
mergeList_
[
blockOffsets_[blocki] + blockPointi
]
] = scaleFactor_ * blockPoints[blockPointi];
const label nx = blocks[blocki].density().x();
const label v1 = blocks[blocki].pointLabel(1, 0, 0);
const label vn = blocks[blocki].pointLabel(nx, 0, 0);
const label vn1 = blocks[blocki].pointLabel(nx-1, 0, 0);
const scalar cwBeg = mag(blockPoints[v1] - blockPoints[v0]);
const scalar cwEnd = mag(blockPoints[vn] - blockPoints[vn1]);
Info<< " i : "
<< cwBeg*scaleTot.x() << " .. "
<< cwEnd*scaleTot.x() << nl;
}
{
const label ny = blocks[blocki].density().y();
const label v1 = blocks[blocki].pointLabel(0, 1, 0);
const label vn = blocks[blocki].pointLabel(0, ny, 0);
const label vn1 = blocks[blocki].pointLabel(0, ny-1, 0);
const scalar cwBeg = mag(blockPoints[v1] - blockPoints[v0]);
const scalar cwEnd = mag(blockPoints[vn] - blockPoints[vn1]);
Info<< " j : "
<< cwBeg*scaleTot.y() << " .. "
<< cwEnd*scaleTot.y() << nl;
}
{
const label nz = blocks[blocki].density().z();
const label v1 = blocks[blocki].pointLabel(0, 0, 1);
const label vn = blocks[blocki].pointLabel(0, 0, nz);
const label vn1 = blocks[blocki].pointLabel(0, 0, nz-1);
const scalar cwBeg = mag(blockPoints[v1] - blockPoints[v0]);
const scalar cwEnd = mag(blockPoints[vn] - blockPoints[vn1]);
Info<< " k : "
<< cwBeg*scaleTot.z() << " .. "
<< cwEnd*scaleTot.z() << nl;
}
Info<< endl;
}
}
inplacePointTransforms(points_);
}
@ -107,24 +128,22 @@ void Foam::blockMesh::createCells() const
Info<< "Creating cells" << endl;
}
cells_.setSize(nCells_);
cells_.resize(nCells_);
label celli = 0;
labelList cellPoints(8); // Hex cells - 8 points
forAll(blocks, blocki)
{
const auto& blockCells = blocks[blocki].cells();
forAll(blockCells, blockCelli)
for (const hexCell& blockCell : blocks[blocki].cells())
{
labelList cellPoints(blockCells[blockCelli].size());
forAll(cellPoints, cellPointi)
{
cellPoints[cellPointi] =
mergeList_
[
blockCells[blockCelli][cellPointi]
blockCell[cellPointi]
+ blockOffsets_[blocki]
];
}
@ -170,7 +189,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
faceList patchFaces(nFaces);
face quadFace(4);
face quad(4);
label faceLabel = 0;
forAll(patchTopologyFaces, patchTopologyFaceLabel)
@ -195,7 +214,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
// Lookup the face points
// and collapse duplicate point labels
quadFace[0] =
quad[0] =
mergeList_
[
blockPatchFaces[blockFaceLabel][0]
@ -211,33 +230,33 @@ Foam::faceList Foam::blockMesh::createPatchFaces
facePointLabel++
)
{
quadFace[nUnique] =
quad[nUnique] =
mergeList_
[
blockPatchFaces[blockFaceLabel][facePointLabel]
+ blockOffsets_[blocki]
];
if (quadFace[nUnique] != quadFace[nUnique-1])
if (quad[nUnique] != quad[nUnique-1])
{
nUnique++;
}
}
if (quadFace[nUnique-1] == quadFace[0])
if (quad[nUnique-1] == quad[0])
{
nUnique--;
}
if (nUnique == 4)
{
patchFaces[faceLabel++] = quadFace;
patchFaces[faceLabel++] = quad;
}
else if (nUnique == 3)
{
patchFaces[faceLabel++] = face
(
labelList::subList(quadFace, 3)
labelList::subList(quad, 3)
);
}
// else the face has collapsed to an edge or point
@ -246,7 +265,7 @@ Foam::faceList Foam::blockMesh::createPatchFaces
}
}
patchFaces.setSize(faceLabel);
patchFaces.resize(faceLabel);
return patchFaces;
}
@ -261,7 +280,7 @@ void Foam::blockMesh::createPatches() const
Info<< "Creating patches" << endl;
}
patches_.setSize(topoPatches.size());
patches_.resize(topoPatches.size());
forAll(topoPatches, patchi)
{
@ -272,6 +291,66 @@ void Foam::blockMesh::createPatches() const
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::polyMesh& Foam::blockMesh::topology() const
{
if (!topologyPtr_)
{
FatalErrorInFunction
<< "topology not allocated"
<< exit(FatalError);
}
return *topologyPtr_;
}
Foam::refPtr<Foam::polyMesh>
Foam::blockMesh::topology(bool applyTransform) const
{
const polyMesh& origTopo = topology();
if (applyTransform && hasPointTransforms())
{
// Copy mesh components
IOobject newIO(origTopo, IOobject::NO_READ, IOobject::NO_WRITE);
newIO.registerObject(false);
pointField newPoints(origTopo.points());
inplacePointTransforms(newPoints);
auto ttopoMesh = refPtr<polyMesh>::New
(
newIO,
std::move(newPoints),
faceList(origTopo.faces()),
labelList(origTopo.faceOwner()),
labelList(origTopo.faceNeighbour())
);
auto& topoMesh = ttopoMesh.ref();
// Patches
const polyBoundaryMesh& pbmOld = origTopo.boundaryMesh();
const polyBoundaryMesh& pbmNew = topoMesh.boundaryMesh();
PtrList<polyPatch> newPatches(pbmOld.size());
forAll(pbmOld, patchi)
{
newPatches.set(patchi, pbmOld[patchi].clone(pbmNew));
}
topoMesh.addPatches(newPatches);
return ttopoMesh;
}
else
{
return origTopo;
}
}
Foam::autoPtr<Foam::polyMesh>
Foam::blockMesh::mesh(const IOobject& io) const
{
@ -285,7 +364,7 @@ Foam::blockMesh::mesh(const IOobject& io) const
auto meshPtr = autoPtr<polyMesh>::New
(
io,
pointField(blkMesh.points()), // Copy, could we re-use space?
pointField(blkMesh.points()), // Use a copy of the block points
blkMesh.cells(),
blkMesh.patches(),
blkMesh.patchNames(),
@ -381,9 +460,9 @@ Foam::blockMesh::mesh(const IOobject& io) const
);
}
pmesh.pointZones().resize(0);
pmesh.faceZones().resize(0);
pmesh.cellZones().resize(0);
pmesh.pointZones().clear();
pmesh.faceZones().clear();
pmesh.cellZones().clear();
pmesh.addZones(List<pointZone*>(), List<faceZone*>(), cz);
}

35
src/mesh/blockMesh/blockMesh/blockMeshMergeGeometrical.C
View File

@ -39,7 +39,7 @@ void Foam::blockMesh::calcGeometricalMerge()
Info<< "Creating block offsets" << endl;
}
blockOffsets_.setSize(blocks.size());
blockOffsets_.resize(blocks.size());
nPoints_ = 0;
nCells_ = 0;
@ -58,22 +58,29 @@ void Foam::blockMesh::calcGeometricalMerge()
Info<< "Creating merge list (geometric search).." << flush;
}
// set unused to -1
mergeList_.setSize(nPoints_);
// Set unused to -1
mergeList_.resize(nPoints_);
mergeList_ = -1;
const pointField& blockPoints = topology().points();
const cellList& blockCells = topology().cells();
const faceList& blockFaces = topology().faces();
const labelList& faceOwnerBlocks = topology().faceOwner();
// Block mesh topology
const polyMesh& topoMesh = topology();
const pointField& blockPoints = topoMesh.points();
const cellList& blockCells = topoMesh.cells();
const faceList& blockFaces = topoMesh.faces();
const labelList& faceOwnerBlocks = topoMesh.faceOwner();
const labelList& faceNeighbourBlocks = topoMesh.faceNeighbour();
const faceList::subList blockinternalFaces
(
blockFaces,
topoMesh.nInternalFaces()
);
// For efficiency, create merge pairs in the first pass
labelListListList glueMergePairs(blockFaces.size());
const labelList& faceNeighbourBlocks = topology().faceNeighbour();
forAll(blockFaces, blockFaceLabel)
{
label blockPlabel = faceOwnerBlocks[blockFaceLabel];
@ -179,7 +186,7 @@ void Foam::blockMesh::calcGeometricalMerge()
sqrMergeTol /= 10.0;
if (topology().isInternalFace(blockFaceLabel))
if (topoMesh.isInternalFace(blockFaceLabel))
{
label blockNlabel = faceNeighbourBlocks[blockFaceLabel];
const pointField& blockNpoints = blocks[blockNlabel].points();
@ -306,12 +313,6 @@ void Foam::blockMesh::calcGeometricalMerge()
}
const faceList::subList blockinternalFaces
(
blockFaces,
topology().nInternalFaces()
);
bool changedPointMerge = false;
label nPasses = 0;

16
src/mesh/blockMesh/blockMesh/blockMeshMergeTopological.C
View File

@ -313,7 +313,7 @@ void Foam::blockMesh::calcTopologicalMerge()
Info<< "Creating block offsets" << endl;
}
blockOffsets_.setSize(blocks.size());
blockOffsets_.resize(blocks.size());
nPoints_ = 0;
nCells_ = 0;
@ -335,11 +335,13 @@ void Foam::blockMesh::calcTopologicalMerge()
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();
const polyMesh& topoMesh = topology();
const pointField& topoPoints = topoMesh.points();
const cellList& topoCells = topoMesh.cells();
const faceList& topoFaces = topoMesh.faces();
const labelList& topoFaceOwn = topoMesh.faceOwner();
const labelList& topoFaceNei = topoMesh.faceNeighbour();
// Topological merging is only necessary for internal block faces
// Note edge and face collapse may apply to boundary faces
@ -347,7 +349,7 @@ void Foam::blockMesh::calcTopologicalMerge()
const faceList::subList topoInternalFaces
(
topoFaces,
topology().nInternalFaces()
topoMesh.nInternalFaces()
);
List<Pair<label>> mergeBlockP(topoInternalFaces.size());

207
src/mesh/blockMesh/blockMesh/blockMeshTopology.C
View File

@ -33,17 +33,27 @@ License
#include "emptyPolyPatch.H"
#include "cyclicPolyPatch.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
namespace Foam
{
// Replace (<block> <face>) face description
// with the corresponding block face
// - otherwise check point labels are in range
template<class Source>
void Foam::blockMesh::checkPatchLabels
static void rewritePatchLabels
(
const Source& source,
const word& patchName,
const pointField& points,
const PtrList<block>& blocks,
const label nPoints,
faceList& patchFaces
) const
)
{
const label nBlocks = blocks.size();
forAll(patchFaces, facei)
{
face& f = patchFaces[facei];
@ -55,40 +65,41 @@ void Foam::blockMesh::checkPatchLabels
const label bi = f[0];
const label fi = f[1];
if (bi >= size())
if (bi < 0 || bi >= nBlocks)
{
FatalIOErrorInFunction(source)
<< "Block index out of range for patch face " << f << nl
<< " Number of blocks = " << size()
<< ", index = " << f[0] << nl
<< "Block index out of range."
<< " Patch face (" << bi << ' ' << fi << ")\n"
<< " Number of blocks = " << nBlocks
<< ", block index = " << bi << nl
<< " on patch " << patchName << ", face " << facei
<< exit(FatalIOError);
}
else if (fi >= operator[](bi).blockShape().faces().size())
else if (fi >= blocks[bi].blockShape().nFaces())
{
FatalIOErrorInFunction(source)
<< "Block face index out of range for patch face " << f
<< nl
<< "Block face index out of range."
<< " Patch face (" << bi << ' ' << fi << ")\n"
<< " Number of block faces = "
<< operator[](bi).blockShape().faces().size()
<< ", index = " << f[1] << nl
<< blocks[bi].blockShape().nFaces()
<< ", face index = " << fi << nl
<< " on patch " << patchName << ", face " << facei
<< exit(FatalIOError);
}
else
{
f = operator[](bi).blockShape().faces()[fi];
f = blocks[bi].blockShape().face(fi);
}
}
else
{
for (const label pointi : f)
{
if (pointi < 0 || pointi >= points.size())
if (pointi < 0 || pointi >= nPoints)
{
FatalIOErrorInFunction(source)
<< "Point label " << pointi
<< " out of range 0.." << points.size() - 1 << nl
<< " out of range 0.." << (nPoints - 1) << nl
<< " on patch " << patchName << ", face " << facei
<< exit(FatalIOError);
}
@ -97,6 +108,10 @@ void Foam::blockMesh::checkPatchLabels
}
}
} // End namespace Foam
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::blockMesh::readPatches
(
@ -112,26 +127,20 @@ void Foam::blockMesh::readPatches
varDict.merge(meshDescription.subOrEmptyDict("namedBlocks"));
ITstream& patchStream(meshDescription.lookup("patches"));
ITstream& patchStream = meshDescription.lookup("patches");
// Read number of patches in mesh
label nPatches = 0;
token firstToken(patchStream);
if (firstToken.isLabel())
if (patchStream.peek().isLabel())
{
nPatches = firstToken.labelToken();
patchStream >> nPatches;
tmpBlocksPatches.setSize(nPatches);
patchNames.setSize(nPatches);
patchTypes.setSize(nPatches);
nbrPatchNames.setSize(nPatches);
}
else
{
patchStream.putBack(firstToken);
}
// Read beginning of blocks
patchStream.readBegin("patches");
@ -175,15 +184,16 @@ void Foam::blockMesh::readPatches
}
}
checkPatchLabels
rewritePatchLabels
(
patchStream,
patchNames[nPatches],
vertices_,
*this,
vertices_.size(),
tmpBlocksPatches[nPatches]
);
nPatches++;
++nPatches;
// Split old style cyclics
@ -298,67 +308,57 @@ void Foam::blockMesh::readBoundary
);
checkPatchLabels
rewritePatchLabels
(
patchInfo.dict(),
patchNames[patchi],
vertices_,
*this,
vertices_.size(),
tmpBlocksPatches[patchi]
);
}
}
void Foam::blockMesh::createCellShapes
(
cellShapeList& tmpBlockCells
)
Foam::cellShapeList Foam::blockMesh::getBlockShapes() const
{
const blockMesh& blocks = *this;
tmpBlockCells.setSize(blocks.size());
cellShapeList shapes(blocks.size());
forAll(blocks, blocki)
{
tmpBlockCells[blocki] = blocks[blocki].blockShape();
shapes[blocki] = blocks[blocki].blockShape();
}
return shapes;
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
Foam::autoPtr<Foam::polyMesh>
Foam::blockMesh::createTopology
(
const IOdictionary& meshDescription,
const word& regionName
)
{
blockList& blocks = *this;
word defaultPatchName = "defaultFaces";
word defaultPatchType = emptyPolyPatch::typeName;
// Read the names/types for the unassigned patch faces
// this is a bit heavy handed (and ugly), but there is currently
// no easy way to rename polyMesh patches subsequently
if (const dictionary* dictPtr = meshDescription.findDict("defaultPatch"))
if (const dictionary* dictptr = meshDescription.findDict("defaultPatch"))
{
dictPtr->readIfPresent("name", defaultPatchName);
dictPtr->readIfPresent("type", defaultPatchType);
dictptr->readIfPresent("name", defaultPatchName);
dictptr->readIfPresent("type", defaultPatchType);
}
// Optional 'scale' factor. Was 'convertToMeters' until OCT-2008
meshDescription.readIfPresentCompat
(
"scale",
{{"convertToMeters", 1012}}, // Mark as changed from 2010 onwards
scaleFactor_
);
// Treat (scale <= 0) as scaling == 1 (no scaling).
if (scaleFactor_ <= 0)
{
scaleFactor_ = 1.0;
}
// Scaling, transformations
readPointTransforms(meshDescription);
// Read the block edges
if (meshDescription.found("edges"))
@ -376,10 +376,14 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
edges_.transfer(edges);
}
else if (verbose_)
else
{
edges_.clear();
if (verbose_)
{
Info<< "No non-linear block edges defined" << endl;
}
}
// Read the block faces
@ -398,10 +402,14 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
faces_.transfer(faces);
}
else if (verbose_)
else
{
faces_.clear();
if (verbose_)
{
Info<< "No non-planar block faces defined" << endl;
}
}
// Create the blocks
@ -420,21 +428,25 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
}
autoPtr<polyMesh> blockMeshPtr;
// Create patch information
faceListList tmpBlocksPatches;
wordList patchNames;
PtrList<dictionary> patchDicts;
// Create the patches
if (verbose_)
{
Info<< "Creating topology patches" << endl;
Info<< nl << "Creating topology patches - ";
}
if (meshDescription.found("patches"))
{
Info<< nl << "Reading patches section" << endl;
if (verbose_)
{
Info<< "from patches section" << endl;
}
faceListList tmpBlocksPatches;
wordList patchNames;
wordList patchTypes;
wordList nbrPatchNames;
@ -447,16 +459,13 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
nbrPatchNames
);
Info<< nl << "Creating block mesh topology" << endl;
if (verbose_)
{
Info<< nl
<< "Reading physicalType from existing boundary file" << endl;
}
cellShapeList tmpBlockCells(blocks.size());
createCellShapes(tmpBlockCells);
Info<< nl << "Reading physicalType from existing boundary file" << endl;
PtrList<dictionary> patchDicts(patchNames.size());
word defaultFacesType;
patchDicts.resize(patchNames.size());
preservePatchTypes
(
@ -496,37 +505,18 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
}
// Override neighbour patch name
if (nbrPatchNames[patchi] != word::null)
if (!nbrPatchNames[patchi].empty())
{
dict.set("neighbourPatch", nbrPatchNames[patchi]);
}
}
blockMeshPtr = autoPtr<polyMesh>::New
(
IOobject
(
regionName,
meshDescription.time().constant(),
meshDescription.time(),
IOobject::NO_READ,
IOobject::NO_WRITE,
false
),
pointField(vertices_), // Copy these points, do NOT move
tmpBlockCells,
tmpBlocksPatches,
patchNames,
patchDicts,
defaultPatchName,
defaultPatchType
);
}
else if (meshDescription.found("boundary"))
{
wordList patchNames;
faceListList tmpBlocksPatches;
PtrList<dictionary> patchDicts;
if (verbose_)
{
Info<< "from boundary section" << endl;
}
readBoundary
(
@ -535,13 +525,27 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
tmpBlocksPatches,
patchDicts
);
}
else
{
if (verbose_)
{
Info<< "with default boundary only!!" << endl;
}
}
Info<< nl << "Creating block mesh topology" << endl;
cellShapeList tmpBlockCells(blocks.size());
createCellShapes(tmpBlockCells);
if (verbose_)
{
Info<< nl << "Creating block mesh topology";
if (hasPointTransforms())
{
Info<< " - scaling/transform applied later";
}
Info<< endl;
}
blockMeshPtr = autoPtr<polyMesh>::New
auto blockMeshPtr = autoPtr<polyMesh>::New
(
IOobject
(
@ -552,15 +556,14 @@ Foam::autoPtr<Foam::polyMesh> Foam::blockMesh::createTopology
IOobject::NO_WRITE,
false
),
pointField(vertices_), // Copy these points, do NOT move
tmpBlockCells,
pointField(vertices_), // Use a copy of vertices
getBlockShapes(),
tmpBlocksPatches,
patchNames,
patchDicts,
defaultPatchName,
defaultPatchType
);
}
check(*blockMeshPtr, meshDescription);

4
src/mesh/blockMesh/blocks/block/block.H
View File

@ -215,14 +215,14 @@ public:
// Access
//- The points for filling the block
inline const pointField& points() const;
inline const pointField& points() const noexcept;
//- The hex cells for filling the block
inline const List<hexCell>& cells() const;
//- The boundary patch faces for the block
inline const FixedList<List<FixedList<label, 4>>, 6>&
boundaryPatches() const;
boundaryPatches() const noexcept;
// Mesh Components

2
src/mesh/blockMesh/blocks/block/blockCreate.C
View File

@ -65,7 +65,7 @@ void Foam::block::createPoints()
// List of edge point and weighting factors
pointField p[12];
scalarList w[12];
label nCurvedEdges = edgesPointsWeights(p, w);
const int nCurvedEdges = edgesPointsWeights(p, w);
points_.resize(nPoints());

4
src/mesh/blockMesh/blocks/block/blockI.H
View File

@ -28,7 +28,7 @@ License
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
inline const Foam::pointField& Foam::block::points() const
inline const Foam::pointField& Foam::block::points() const noexcept
{
return points_;
}
@ -46,7 +46,7 @@ inline const Foam::List<Foam::hexCell>& Foam::block::cells() const
inline const Foam::FixedList<Foam::List<Foam::FixedList<Foam::label, 4>>, 6>&
Foam::block::boundaryPatches() const
Foam::block::boundaryPatches() const noexcept
{
return blockPatches_;
}