songtianlun/openfoam
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

blockMesh::blockDescriptor: Generate the edge points and weights on demand

Browse Source
This commit is contained in:
Henry Weller 2016-10-09 15:16:41 +01:00
parent 2997bfbdfd
commit e53520ecbf
6 changed files with 131 additions and 170 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

144
src/mesh/blockMesh/block/blockCreate.C
View File

@ -44,10 +44,10 @@ void Foam::block::createPoints() const
const point& p111 = blockPoint(6);
const point& p011 = blockPoint(7);
// List of edge point and weighting factors
const List<List<point>>& p = blockEdgePoints();
const scalarListList& w = blockEdgeWeights();
List<point> p[12];
scalarList w[12];
edgesPointsWeights(p, w);
//
// Generate vertices
@ -63,51 +63,34 @@ void Foam::block::createPoints() const
{
const label vertexNo = vtxLabel(i, j, k);
// 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];
vector edgex4 = p001 + (p101 - p001)*w[3][i];
vector edgey1 = p000 + (p010 - p000)*w[4][j];
vector edgey2 = p100 + (p110 - p100)*w[5][j];
vector edgey3 = p101 + (p111 - p101)*w[6][j];
vector edgey4 = p001 + (p011 - p001)*w[7][j];
vector edgez1 = p000 + (p001 - p000)*w[8][k];
vector edgez2 = p100 + (p101 - p100)*w[9][k];
vector edgez3 = p110 + (p111 - p110)*w[10][k];
vector edgez4 = p010 + (p011 - p010)*w[11][k];
// Calculate the importance factors for all edges
// x-direction
scalar impx1 =
(
(1.0 - w[0][i])*(1.0 - w[4][j])*(1.0 - w[8][k])
+ w[0][i]*(1.0 - w[5][j])*(1.0 - w[9][k])
(1 - w[0][i])*(1 - w[4][j])*(1 - w[8][k])
+ w[0][i]*(1 - w[5][j])*(1 - w[9][k])
);
scalar impx2 =
(
(1.0 - w[1][i])*w[4][j]*(1.0 - w[11][k])
+ w[1][i]*w[5][j]*(1.0 - w[10][k])
(1 - w[1][i])*w[4][j]*(1 - w[11][k])
+ w[1][i]*w[5][j]*(1 - w[10][k])
);
scalar impx3 =
(
(1.0 - w[2][i])*w[7][j]*w[11][k]
(1 - w[2][i])*w[7][j]*w[11][k]
+ w[2][i]*w[6][j]*w[10][k]
);
scalar impx4 =
(
(1.0 - w[3][i])*(1.0 - w[7][j])*w[8][k]
+ w[3][i]*(1.0 - w[6][j])*w[9][k]
(1 - w[3][i])*(1 - w[7][j])*w[8][k]
+ w[3][i]*(1 - w[6][j])*w[9][k]
);
scalar magImpx = impx1 + impx2 + impx3 + impx4;
const scalar magImpx = impx1 + impx2 + impx3 + impx4;
impx1 /= magImpx;
impx2 /= magImpx;
impx3 /= magImpx;
@ -117,29 +100,29 @@ void Foam::block::createPoints() const
// y-direction
scalar impy1 =
(
(1.0 - w[4][j])*(1.0 - w[0][i])*(1.0 - w[8][k])
+ w[4][j]*(1.0 - w[1][i])*(1.0 - w[11][k])
(1 - w[4][j])*(1 - w[0][i])*(1 - w[8][k])
+ w[4][j]*(1 - w[1][i])*(1 - w[11][k])
);
scalar impy2 =
(
(1.0 - w[5][j])*w[0][i]*(1.0 - w[9][k])
+ w[5][j]*w[1][i]*(1.0 - w[10][k])
(1 - w[5][j])*w[0][i]*(1 - w[9][k])
+ w[5][j]*w[1][i]*(1 - w[10][k])
);
scalar impy3 =
(
(1.0 - w[6][j])*w[3][i]*w[9][k]
(1 - w[6][j])*w[3][i]*w[9][k]
+ w[6][j]*w[2][i]*w[10][k]
);
scalar impy4 =
(
(1.0 - w[7][j])*(1.0 - w[3][i])*w[8][k]
+ w[7][j]*(1.0 - w[2][i])*w[11][k]
(1 - w[7][j])*(1 - w[3][i])*w[8][k]
+ w[7][j]*(1 - w[2][i])*w[11][k]
);
scalar magImpy = impy1 + impy2 + impy3 + impy4;
const scalar magImpy = impy1 + impy2 + impy3 + impy4;
impy1 /= magImpy;
impy2 /= magImpy;
impy3 /= magImpy;
@ -149,80 +132,75 @@ void Foam::block::createPoints() const
// z-direction
scalar impz1 =
(
(1.0 - w[8][k])*(1.0 - w[0][i])*(1.0 - w[4][j])
+ w[8][k]*(1.0 - w[3][i])*(1.0 - w[7][j])
(1 - w[8][k])*(1 - w[0][i])*(1 - w[4][j])
+ w[8][k]*(1 - w[3][i])*(1 - w[7][j])
);
scalar impz2 =
(
(1.0 - w[9][k])*w[0][i]*(1.0 - w[5][j])
+ w[9][k]*w[3][i]*(1.0 - w[6][j])
(1 - w[9][k])*w[0][i]*(1 - w[5][j])
+ w[9][k]*w[3][i]*(1 - w[6][j])
);
scalar impz3 =
(
(1.0 - w[10][k])*w[1][i]*w[5][j]
(1 - w[10][k])*w[1][i]*w[5][j]
+ w[10][k]*w[2][i]*w[6][j]
);
scalar impz4 =
(
(1.0 - w[11][k])*(1.0 - w[1][i])*w[4][j]
+ w[11][k]*(1.0 - w[2][i])*w[7][j]
(1 - w[11][k])*(1 - w[1][i])*w[4][j]
+ w[11][k]*(1 - w[2][i])*w[7][j]
);
scalar magImpz = impz1 + impz2 + impz3 + impz4;
const scalar magImpz = impz1 + impz2 + impz3 + impz4;
impz1 /= magImpz;
impz2 /= magImpz;
impz3 /= magImpz;
impz4 /= magImpz;
// 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);
vector corx4 = impx4*(p[3][i] - edgex4);
// Points on straight edges
const vector edgex1 = p000 + (p100 - p000)*w[0][i];
const vector edgex2 = p010 + (p110 - p010)*w[1][i];
const vector edgex3 = p011 + (p111 - p011)*w[2][i];
const vector edgex4 = p001 + (p101 - p001)*w[3][i];
vector cory1 = impy1*(p[4][j] - edgey1);
vector cory2 = impy2*(p[5][j] - edgey2);
vector cory3 = impy3*(p[6][j] - edgey3);
vector cory4 = impy4*(p[7][j] - edgey4);
vector corz1 = impz1*(p[8][k] - edgez1);
vector corz2 = impz2*(p[9][k] - edgez2);
vector corz3 = impz3*(p[10][k] - edgez3);
vector corz4 = impz4*(p[11][k] - edgez4);
// Multiply by the importance factor
// x-direction
edgex1 *= impx1;
edgex2 *= impx2;
edgex3 *= impx3;
edgex4 *= impx4;
// y-direction
edgey1 *= impy1;
edgey2 *= impy2;
edgey3 *= impy3;
edgey4 *= impy4;
// z-direction
edgez1 *= impz1;
edgez2 *= impz2;
edgez3 *= impz3;
edgez4 *= impz4;
const vector edgey1 = p000 + (p010 - p000)*w[4][j];
const vector edgey2 = p100 + (p110 - p100)*w[5][j];
const vector edgey3 = p101 + (p111 - p101)*w[6][j];
const vector edgey4 = p001 + (p011 - p001)*w[7][j];
const vector edgez1 = p000 + (p001 - p000)*w[8][k];
const vector edgez2 = p100 + (p101 - p100)*w[9][k];
const vector edgez3 = p110 + (p111 - p110)*w[10][k];
const vector edgez4 = p010 + (p011 - p010)*w[11][k];
// Add the contributions
points_[vertexNo] =
(
edgex1 + edgex2 + edgex3 + edgex4
+ edgey1 + edgey2 + edgey3 + edgey4
+ edgez1 + edgez2 + edgez3 + edgez4
) / 3.0;
impx1*edgex1 + impx2*edgex2 + impx3*edgex3 + impx4*edgex4
+ impy1*edgey1 + impy2*edgey2 + impy3*edgey3 + impy4*edgey4
+ impz1*edgez1 + impz2*edgez2 + impz3*edgez3 + impz4*edgez4
)/3.0;
// Calculate the correction vectors
const vector corx1 = impx1*(p[0][i] - edgex1);
const vector corx2 = impx2*(p[1][i] - edgex2);
const vector corx3 = impx3*(p[2][i] - edgex3);
const vector corx4 = impx4*(p[3][i] - edgex4);
const vector cory1 = impy1*(p[4][j] - edgey1);
const vector cory2 = impy2*(p[5][j] - edgey2);
const vector cory3 = impy3*(p[6][j] - edgey3);
const vector cory4 = impy4*(p[7][j] - edgey4);
const vector corz1 = impz1*(p[8][k] - edgez1);
const vector corz2 = impz2*(p[9][k] - edgez2);
const vector corz3 = impz3*(p[10][k] - edgez3);
const vector corz4 = impz4*(p[11][k] - edgez4);
points_[vertexNo] +=
(

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

@ -118,8 +118,6 @@ Foam::blockDescriptor::blockDescriptor
edges_(edges),
blockShape_(bshape),
density_(density),
edgePoints_(12),
edgeWeights_(12),
expand_(expand),
zoneName_(zoneName)
{
@ -129,9 +127,6 @@ Foam::blockDescriptor::blockDescriptor
<< "Unknown definition of expansion ratios"
<< exit(FatalError);
}
// Create a list of edges
makeBlockEdges();
}
@ -146,8 +141,6 @@ Foam::blockDescriptor::blockDescriptor
edges_(edges),
blockShape_(is),
density_(),
edgePoints_(12),
edgeWeights_(12),
expand_
(
12,
@ -238,9 +231,6 @@ Foam::blockDescriptor::blockDescriptor
}
check(is);
// Create a list of edges
makeBlockEdges();
}

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

@ -90,12 +90,6 @@ class blockDescriptor
//- The number of cells in the i,j,k directions
Vector<label> density_;
//- Block edge points
List<List<point>> edgePoints_;
//- Block edge weighting factors
scalarListList edgeWeights_;
//- Expansion ratios in all directions
List<gradingDescriptors> expand_;
@ -108,11 +102,16 @@ class blockDescriptor
//- Check block has outward-pointing faces
void check(const Istream& is);
//- Set the points/weights for all edges
void makeBlockEdges();
//- Set the edge points/weights
void setEdge(label edgei, label start, label end, label dim);
//- Calculate the points and weights for the specified edge
void edgePointsWeights
(
List<point> (&edgePoints)[12],
scalarList (&edgeWeights)[12],
const label edgei,
const label start,
const label end,
const label dim
) const;
// Private Member Functions
@ -156,12 +155,6 @@ public:
//- Return the block shape
inline const cellShape& blockShape() const;
//- Return the block points along each edge
inline const List<List<point>>& blockEdgePoints() const;
//- Return the weightings along each edge
inline const scalarListList& blockEdgeWeights() const;
//- Return the mesh density (number of cells) in the i,j,k directions
inline const Vector<label>& density() const;
@ -177,6 +170,13 @@ public:
//- Return block point for local label i
inline const point& blockPoint(const label i) const;
//- Calculate the points and weights for all edges
void edgesPointsWeights
(
List<point> (&edgePoints)[12],
scalarList (&edgeWeights)[12]
) const;
// IOstream Operators

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

@ -29,41 +29,15 @@ License
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
void Foam::blockDescriptor::makeBlockEdges()
{
const label ni = density_.x();
const label nj = density_.y();
const label nk = density_.z();
// These edges correspond to the "hex" cellModel
// X-direction
setEdge(0, 0, 1, ni);
setEdge(1, 3, 2, ni);
setEdge(2, 7, 6, ni);
setEdge(3, 4, 5, ni);
// Y-direction
setEdge(4, 0, 3, nj);
setEdge(5, 1, 2, nj);
setEdge(6, 5, 6, nj);
setEdge(7, 4, 7, nj);
// Z-direction
setEdge(8, 0, 4, nk);
setEdge(9, 1, 5, nk);
setEdge(10, 2, 6, nk);
setEdge(11, 3, 7, nk);
}
void Foam::blockDescriptor::setEdge
void Foam::blockDescriptor::edgePointsWeights
(
label edgei,
label start,
label end,
label nDiv
)
List<point> (&edgePoints)[12],
scalarList (&edgeWeights)[12],
const label edgei,
const label start,
const label end,
const label nDiv
) const
{
// Set reference to the list of labels defining the block
const labelList& blockLabels = blockShape_;
@ -78,7 +52,7 @@ void Foam::blockDescriptor::setEdge
{
const blockEdge& cedge = edges_[cedgei];
int cmp = cedge.compare(blockLabels[start], blockLabels[end]);
const int cmp = cedge.compare(blockLabels[start], blockLabels[end]);
if (cmp)
{
@ -87,29 +61,29 @@ void Foam::blockDescriptor::setEdge
// Curve has the same orientation
// Divide the line
lineDivide divEdge(cedge, nDiv, expand_[edgei]);
const lineDivide divEdge(cedge, nDiv, expand_[edgei]);
edgePoints_[edgei] = divEdge.points();
edgeWeights_[edgei] = divEdge.lambdaDivisions();
edgePoints[edgei] = divEdge.points();
edgeWeights[edgei] = divEdge.lambdaDivisions();
}
else
{
// Curve has the opposite orientation
// Divide the line
lineDivide divEdge(cedge, nDiv, expand_[edgei].inv());
const lineDivide divEdge(cedge, nDiv, expand_[edgei].inv());
const pointField& p = divEdge.points();
const scalarList& d = divEdge.lambdaDivisions();
edgePoints_[edgei].setSize(p.size());
edgeWeights_[edgei].setSize(d.size());
edgePoints[edgei].setSize(p.size());
edgeWeights[edgei].setSize(d.size());
label pMax = p.size() - 1;
forAll(p, pI)
{
edgePoints_[edgei][pI] = p[pMax - pI];
edgeWeights_[edgei][pI] = 1.0 - d[pMax - pI];
edgePoints[edgei][pI] = p[pMax - pI];
edgeWeights[edgei][pI] = 1.0 - d[pMax - pI];
}
}
@ -127,8 +101,41 @@ void Foam::blockDescriptor::setEdge
expand_[edgei]
);
edgePoints_[edgei] = divEdge.points();
edgeWeights_[edgei] = divEdge.lambdaDivisions();
edgePoints[edgei] = divEdge.points();
edgeWeights[edgei] = divEdge.lambdaDivisions();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::blockDescriptor::edgesPointsWeights
(
List<point> (&edgePoints)[12],
scalarList (&edgeWeights)[12]
) const
{
// These edges correspond to the "hex" cellModel
// X-direction
const label ni = density_.x();
edgePointsWeights(edgePoints, edgeWeights, 0, 0, 1, ni);
edgePointsWeights(edgePoints, edgeWeights, 1, 3, 2, ni);
edgePointsWeights(edgePoints, edgeWeights, 2, 7, 6, ni);
edgePointsWeights(edgePoints, edgeWeights, 3, 4, 5, ni);
// Y-direction
const label nj = density_.y();
edgePointsWeights(edgePoints, edgeWeights, 4, 0, 3, nj);
edgePointsWeights(edgePoints, edgeWeights, 5, 1, 2, nj);
edgePointsWeights(edgePoints, edgeWeights, 6, 5, 6, nj);
edgePointsWeights(edgePoints, edgeWeights, 7, 4, 7, nj);
// Z-direction
const label nk = density_.z();
edgePointsWeights(edgePoints, edgeWeights, 8, 0, 4, nk);
edgePointsWeights(edgePoints, edgeWeights, 9, 1, 5, nk);
edgePointsWeights(edgePoints, edgeWeights, 10, 2, 6, nk);
edgePointsWeights(edgePoints, edgeWeights, 11, 3, 7, nk);
}

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

@ -31,20 +31,6 @@ inline const Foam::cellShape& Foam::blockDescriptor::blockShape() const
}
inline const Foam::List<Foam::List<Foam::point>>&
Foam::blockDescriptor::blockEdgePoints() const
{
return edgePoints_;
}
inline const Foam::scalarListList&
Foam::blockDescriptor::blockEdgeWeights() const
{
return edgeWeights_;
}
inline const Foam::Vector<Foam::label>& Foam::blockDescriptor::density() const
{
return density_;

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

@ -48,7 +48,7 @@ Foam::blockMesh::blockMesh(const IOdictionary& dict, const word& regionName)
if (fastMerge)
{
calcMergeInfoFast();
calcMergeInfoFast();
}
else
{