diff --git a/applications/utilities/mesh/generation/blockMesh/Make/options b/applications/utilities/mesh/generation/blockMesh/Make/options index 2e07c477d..9d13a25de 100644 --- a/applications/utilities/mesh/generation/blockMesh/Make/options +++ b/applications/utilities/mesh/generation/blockMesh/Make/options @@ -1,9 +1,11 @@ EXE_INC = \ -I$(LIB_SRC)/mesh/blockMesh/lnInclude \ -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/fileFormats/lnInclude \ -I$(LIB_SRC)/dynamicMesh/lnInclude EXE_LIBS = \ -lblockMesh \ -lmeshTools \ + -lfileFormats \ -ldynamicMesh diff --git a/applications/utilities/mesh/generation/blockMesh/blockMesh.C b/applications/utilities/mesh/generation/blockMesh/blockMesh.C index 2be39caae..a7860c8f5 100644 --- a/applications/utilities/mesh/generation/blockMesh/blockMesh.C +++ b/applications/utilities/mesh/generation/blockMesh/blockMesh.C @@ -176,8 +176,6 @@ int main(int argc, char *argv[]) Info<< "Creating block mesh from\n " << meshDictIO.objectPath() << endl; - blockMesh::verbose(true); - IOdictionary meshDict(meshDictIO); blockMesh blocks(meshDict, regionName); @@ -286,8 +284,8 @@ int main(int argc, char *argv[]) forAll(blocks, blockI) { const block& b = blocks[blockI]; - const labelListList& blockCells = b.cells(); - const word& zoneName = b.blockDef().zoneName(); + const List> blockCells = b.cells(); + const word& zoneName = b.zoneName(); if (zoneName.size()) { diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMesh.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMesh.C index a3ccb7589..3aea81a88 100644 --- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMesh.C +++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMesh.C @@ -75,7 +75,7 @@ void Foam::vtkPV3blockMesh::updateInfoBlocks const int nBlocks = blkMesh.size(); for (int blockI = 0; blockI < nBlocks; ++blockI) { - const blockDescriptor& blockDef = blkMesh[blockI].blockDef(); + const blockDescriptor& blockDef = blkMesh[blockI]; word partName = Foam::name(blockI); diff --git a/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMeshConvert.C b/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMeshConvert.C index 876fc47a3..8dce33a98 100644 --- a/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMeshConvert.C +++ b/applications/utilities/postProcessing/graphics/PV3Readers/PV3blockMeshReader/vtkPV3blockMesh/vtkPV3blockMeshConvert.C @@ -77,7 +77,7 @@ void Foam::vtkPV3blockMesh::convertMeshBlocks continue; } - const blockDescriptor& blockDef = blkMesh[blockI].blockDef(); + const blockDescriptor& blockDef = blkMesh[blockI]; vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New(); @@ -168,7 +168,7 @@ void Foam::vtkPV3blockMesh::convertMeshEdges // search each block forAll(blkMesh, blockI) { - const blockDescriptor& blockDef = blkMesh[blockI].blockDef(); + const blockDescriptor& blockDef = blkMesh[blockI]; edgeList blkEdges = blockDef.blockShape().edges(); diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/Make/options b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/Make/options index ccbea978f..9dcba4b39 100644 --- a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/Make/options +++ b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/Make/options @@ -1,5 +1,6 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/fileFormats/lnInclude \ -I$(LIB_SRC)/mesh/blockMesh/lnInclude \ -I$(ParaView_INCLUDE_DIR) \ -I$(ParaView_INCLUDE_DIR)/vtkkwiml \ @@ -8,6 +9,7 @@ EXE_INC = \ LIB_LIBS = \ -lmeshTools \ + -lfileFormats \ -lblockMesh \ -L$(FOAM_LIBBIN) -lvtkPVReaders \ $(GLIBS) diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMesh.C b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMesh.C index 22cdf2db5..a9932ffcc 100644 --- a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMesh.C +++ b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMesh.C @@ -75,7 +75,7 @@ void Foam::vtkPVblockMesh::updateInfoBlocks const int nBlocks = blkMesh.size(); for (int blockI = 0; blockI < nBlocks; ++blockI) { - const blockDescriptor& blockDef = blkMesh[blockI].blockDef(); + const blockDescriptor& blockDef = blkMesh[blockI]; word partName = Foam::name(blockI); diff --git a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMeshConvert.C b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMeshConvert.C index a287be2dd..9c7f8a06b 100644 --- a/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMeshConvert.C +++ b/applications/utilities/postProcessing/graphics/PVReaders/PVblockMeshReader/vtkPVblockMesh/vtkPVblockMeshConvert.C @@ -77,7 +77,7 @@ void Foam::vtkPVblockMesh::convertMeshBlocks continue; } - const blockDescriptor& blockDef = blkMesh[blockI].blockDef(); + const blockDescriptor& blockDef = blkMesh[blockI]; vtkUnstructuredGrid* vtkmesh = vtkUnstructuredGrid::New(); @@ -168,12 +168,12 @@ void Foam::vtkPVblockMesh::convertMeshEdges // search each block forAll(blkMesh, blockI) { - const blockDescriptor& blockDef = blkMesh[blockI].blockDef(); + const blockDescriptor& blockDef = blkMesh[blockI]; edgeList blkEdges = blockDef.blockShape().edges(); // List of edge point and weighting factors - List edgesPoints[12]; + pointField edgesPoints[12]; scalarList edgesWeights[12]; blockDef.edgesPointsWeights(edgesPoints, edgesWeights); diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files index 27bed6559..98e1c9cf6 100644 --- a/src/mesh/blockMesh/Make/files +++ b/src/mesh/blockMesh/Make/files @@ -1,14 +1,16 @@ -blockEdges/BSpline.C -blockEdges/CatmullRomSpline.C -blockEdges/polyLine.C +blockEdges/blockEdge/blockEdge.C +blockEdges/lineDivide/lineDivide.C +blockEdges/lineEdge/lineEdge.C +blockEdges/polyLineEdge/polyLine.C +blockEdges/polyLineEdge/polyLineEdge.C +blockEdges/arcEdge/arcEdge.C +blockEdges/BSplineEdge/BSpline.C +blockEdges/BSplineEdge/BSplineEdge.C +blockEdges/splineEdge/CatmullRomSpline.C +blockEdges/splineEdge/splineEdge.C -blockEdges/arcEdge.C -blockEdges/blockEdge.C -blockEdges/lineEdge.C -blockEdges/polyLineEdge.C -blockEdges/lineDivide.C -blockEdges/BSplineEdge.C -blockEdges/splineEdge.C +blockFaces/blockFace/blockFace.C +blockFaces/projectFace/projectFace.C gradingDescriptor/gradingDescriptor.C gradingDescriptor/gradingDescriptors.C diff --git a/src/mesh/blockMesh/Make/options b/src/mesh/blockMesh/Make/options index 1618ab57e..3be93bd21 100644 --- a/src/mesh/blockMesh/Make/options +++ b/src/mesh/blockMesh/Make/options @@ -1,7 +1,7 @@ EXE_INC = \ -I$(LIB_SRC)/meshTools/lnInclude \ - -I$(LIB_SRC)/dynamicMesh/lnInclude + -I$(LIB_SRC)/fileFormats/lnInclude LIB_LIBS = \ -lmeshTools \ - -ldynamicMesh + -lfileFormats diff --git a/src/mesh/blockMesh/block/block.C b/src/mesh/blockMesh/block/block.C index 0b345511f..1722b0095 100644 --- a/src/mesh/blockMesh/block/block.C +++ b/src/mesh/blockMesh/block/block.C @@ -31,63 +31,23 @@ Foam::block::block ( const pointField& vertices, const blockEdgeList& edges, + const blockFaceList& faces, Istream& is ) : - blockDescriptor(vertices, edges, is), - points_(0), - cells_(0), - boundaryPatches_(0) -{} + blockDescriptor(vertices, edges, faces, is) +{ + createPoints(); + createBoundary(); +} Foam::block::block(const blockDescriptor& blockDesc) : - blockDescriptor(blockDesc), - points_(0), - cells_(0), - boundaryPatches_(0) -{} - - -// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // - -Foam::block::~block() -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -const Foam::pointField& Foam::block::points() const + blockDescriptor(blockDesc) { - if (points_.empty()) - { - createPoints(); - } - - return points_; -} - - -const Foam::labelListList& Foam::block::cells() const -{ - if (cells_.empty()) - { - createCells(); - } - - return cells_; -} - - -const Foam::labelListListList& Foam::block::boundaryPatches() const -{ - if (boundaryPatches_.empty()) - { - createBoundary(); - } - - return boundaryPatches_; + createPoints(); + createBoundary(); } diff --git a/src/mesh/blockMesh/block/block.H b/src/mesh/blockMesh/block/block.H index 0f48b975e..7b55c5f59 100644 --- a/src/mesh/blockMesh/block/block.H +++ b/src/mesh/blockMesh/block/block.H @@ -68,25 +68,19 @@ class block // Private data //- List of points - mutable pointField points_; - - //- List of cells - mutable labelListList cells_; + pointField points_; //- Boundary patches - mutable labelListListList boundaryPatches_; + FixedList>, 6> boundaryPatches_; // Private Member Functions //- Creates vertices for cells filling the block - void createPoints() const; - - //- Creates cells for filling the block - void createCells() const; + void createPoints(); //- Creates boundary patch faces for the block - void createBoundary() const; + void createBoundary(); //- Disallow default bitwise copy construct block(const block&); @@ -94,6 +88,7 @@ class block //- Disallow default bitwise assignment void operator=(const block&); + public: // Constructors @@ -103,6 +98,7 @@ public: ( const pointField& vertices, const blockEdgeList& edges, + const blockFaceList& faces, Istream& is ); @@ -122,50 +118,42 @@ public: { const pointField& points_; const blockEdgeList& edges_; + const blockFaceList& faces_; public: - iNew(const pointField& points, const blockEdgeList& edges) + iNew + ( + const pointField& points, + const blockEdgeList& edges, + const blockFaceList& faces + ) : points_(points), - edges_(edges) + edges_(edges), + faces_(faces) {} autoPtr operator()(Istream& is) const { - return autoPtr(new block(points_, edges_, is)); + return autoPtr(new block(points_, edges_, faces_, is)); } }; - //- Destructor - ~block(); - - // Member Functions // Access - //- Return the block definition - inline const blockDescriptor& blockDef() const; - - //- Vertex label offset for a particular i,j,k position - inline label vtxLabel(label i, label j, label k) const; - //- Return the points for filling the block - const pointField& points() const; + inline const pointField& points() const; //- Return the cells for filling the block - const labelListList& cells() const; + List> cells() const; //- Return the boundary patch faces for the block - const labelListListList& boundaryPatches() const; - - - // Edit - - //- Clear geometry (internal points, cells, boundaryPatches) - void clearGeom(); + inline const FixedList>, 6>& + boundaryPatches() const; // Ostream Operator diff --git a/src/mesh/blockMesh/block/blockCreate.C b/src/mesh/blockMesh/block/blockCreate.C index 0aac98f19..8f63d629f 100644 --- a/src/mesh/blockMesh/block/blockCreate.C +++ b/src/mesh/blockMesh/block/blockCreate.C @@ -27,7 +27,22 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::block::createPoints() const +#define w0 w[0][i] +#define w1 w[1][i] +#define w2 w[2][i] +#define w3 w[3][i] + +#define w4 w[4][j] +#define w5 w[5][j] +#define w6 w[6][j] +#define w7 w[7][j] + +#define w8 w[8][k] +#define w9 w[9][k] +#define w10 w[10][k] +#define w11 w[11][k] + +void Foam::block::createPoints() { // Set local variables for mesh specification const label ni = density().x(); @@ -45,168 +60,285 @@ void Foam::block::createPoints() const const point& p011 = blockPoint(7); // List of edge point and weighting factors - List p[12]; + pointField p[12]; scalarList w[12]; - edgesPointsWeights(p, w); + label nCurvedEdges = edgesPointsWeights(p, w); - // - // Generate vertices - // - points_.clear(); points_.setSize(nPoints()); - for (label k = 0; k <= nk; k++) - { - for (label j = 0; j <= nj; j++) - { - for (label i = 0; i <= ni; i++) - { - const label vertexNo = vtxLabel(i, j, k); + points_[pointLabel(0, 0, 0)] = p000; + points_[pointLabel(ni, 0, 0)] = p100; + points_[pointLabel(ni, nj, 0)] = p110; + points_[pointLabel(0, nj, 0)] = p010; + points_[pointLabel(0, 0, nk)] = p001; + points_[pointLabel(ni, 0, nk)] = p101; + points_[pointLabel(ni, nj, nk)] = p111; + points_[pointLabel(0, nj, nk)] = p011; - // Calculate the importance factors for all edges + for (label k=0; k<=nk; k++) + { + for (label j=0; j<=nj; j++) + { + for (label i=0; i<=ni; i++) + { + // Skip block vertices + if (vertex(i, j, k)) continue; + + const label vijk = pointLabel(i, j, k); + + // Calculate the weighting factors for all edges // x-direction - scalar impx1 = - ( - (1 - w[0][i])*(1 - w[4][j])*(1 - w[8][k]) - + w[0][i]*(1 - w[5][j])*(1 - w[9][k]) - ); + scalar wx1 = (1 - w0)*(1 - w4)*(1 - w8) + w0*(1 - w5)*(1 - w9); + scalar wx2 = (1 - w1)*w4*(1 - w11) + w1*w5*(1 - w10); + scalar wx3 = (1 - w2)*w7*w11 + w2*w6*w10; + scalar wx4 = (1 - w3)*(1 - w7)*w8 + w3*(1 - w6)*w9; - scalar impx2 = - ( - (1 - w[1][i])*w[4][j]*(1 - w[11][k]) - + w[1][i]*w[5][j]*(1 - w[10][k]) - ); - - scalar impx3 = - ( - (1 - w[2][i])*w[7][j]*w[11][k] - + w[2][i]*w[6][j]*w[10][k] - ); - - scalar impx4 = - ( - (1 - w[3][i])*(1 - w[7][j])*w[8][k] - + w[3][i]*(1 - w[6][j])*w[9][k] - ); - - const scalar magImpx = impx1 + impx2 + impx3 + impx4; - impx1 /= magImpx; - impx2 /= magImpx; - impx3 /= magImpx; - impx4 /= magImpx; + const scalar sumWx = wx1 + wx2 + wx3 + wx4; + wx1 /= sumWx; + wx2 /= sumWx; + wx3 /= sumWx; + wx4 /= sumWx; // y-direction - scalar impy1 = - ( - (1 - w[4][j])*(1 - w[0][i])*(1 - w[8][k]) - + w[4][j]*(1 - w[1][i])*(1 - w[11][k]) - ); + scalar wy1 = (1 - w4)*(1 - w0)*(1 - w8) + w4*(1 - w1)*(1 - w11); + scalar wy2 = (1 - w5)*w0*(1 - w9) + w5*w1*(1 - w10); + scalar wy3 = (1 - w6)*w3*w9 + w6*w2*w10; + scalar wy4 = (1 - w7)*(1 - w3)*w8 + w7*(1 - w2)*w11; - scalar impy2 = - ( - (1 - w[5][j])*w[0][i]*(1 - w[9][k]) - + w[5][j]*w[1][i]*(1 - w[10][k]) - ); - - scalar impy3 = - ( - (1 - w[6][j])*w[3][i]*w[9][k] - + w[6][j]*w[2][i]*w[10][k] - ); - - scalar impy4 = - ( - (1 - w[7][j])*(1 - w[3][i])*w[8][k] - + w[7][j]*(1 - w[2][i])*w[11][k] - ); - - const scalar magImpy = impy1 + impy2 + impy3 + impy4; - impy1 /= magImpy; - impy2 /= magImpy; - impy3 /= magImpy; - impy4 /= magImpy; + const scalar sumWy = wy1 + wy2 + wy3 + wy4; + wy1 /= sumWy; + wy2 /= sumWy; + wy3 /= sumWy; + wy4 /= sumWy; // z-direction - scalar impz1 = - ( - (1 - w[8][k])*(1 - w[0][i])*(1 - w[4][j]) - + w[8][k]*(1 - w[3][i])*(1 - w[7][j]) - ); + scalar wz1 = (1 - w8)*(1 - w0)*(1 - w4) + w8*(1 - w3)*(1 - w7); + scalar wz2 = (1 - w9)*w0*(1 - w5) + w9*w3*(1 - w6); + scalar wz3 = (1 - w10)*w1*w5 + w10*w2*w6; + scalar wz4 = (1 - w11)*(1 - w1)*w4 + w11*(1 - w2)*w7; - scalar impz2 = - ( - (1 - w[9][k])*w[0][i]*(1 - w[5][j]) - + w[9][k]*w[3][i]*(1 - w[6][j]) - ); - - scalar impz3 = - ( - (1 - w[10][k])*w[1][i]*w[5][j] - + w[10][k]*w[2][i]*w[6][j] - ); - - scalar impz4 = - ( - (1 - w[11][k])*(1 - w[1][i])*w[4][j] - + w[11][k]*(1 - w[2][i])*w[7][j] - ); - - const scalar magImpz = impz1 + impz2 + impz3 + impz4; - impz1 /= magImpz; - impz2 /= magImpz; - impz3 /= magImpz; - impz4 /= magImpz; + const scalar sumWz = wz1 + wz2 + wz3 + wz4; + wz1 /= sumWz; + wz2 /= sumWz; + wz3 /= sumWz; + wz4 /= sumWz; // 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]; + const vector edgex1 = p000 + (p100 - p000)*w0; + const vector edgex2 = p010 + (p110 - p010)*w1; + const vector edgex3 = p011 + (p111 - p011)*w2; + const vector edgex4 = p001 + (p101 - p001)*w3; - 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 edgey1 = p000 + (p010 - p000)*w4; + const vector edgey2 = p100 + (p110 - p100)*w5; + const vector edgey3 = p101 + (p111 - p101)*w6; + const vector edgey4 = p001 + (p011 - p001)*w7; - 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]; + const vector edgez1 = p000 + (p001 - p000)*w8; + const vector edgez2 = p100 + (p101 - p100)*w9; + const vector edgez3 = p110 + (p111 - p110)*w10; + const vector edgez4 = p010 + (p011 - p010)*w11; // Add the contributions - points_[vertexNo] = + points_[vijk] = ( - 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; + wx1*edgex1 + wx2*edgex2 + wx3*edgex3 + wx4*edgex4 + + wy1*edgey1 + wy2*edgey2 + wy3*edgey3 + wy4*edgey4 + + wz1*edgez1 + wz2*edgez2 + wz3*edgez3 + wz4*edgez4 + )/3; - // 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); + // Apply curved-edge correction if block has curved edges + if (nCurvedEdges) + { + // Calculate the correction vectors + const vector corx1 = wx1*(p[0][i] - edgex1); + const vector corx2 = wx2*(p[1][i] - edgex2); + const vector corx3 = wx3*(p[2][i] - edgex3); + const vector corx4 = wx4*(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 cory1 = wy1*(p[4][j] - edgey1); + const vector cory2 = wy2*(p[5][j] - edgey2); + const vector cory3 = wy3*(p[6][j] - edgey3); + const vector cory4 = wy4*(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); + const vector corz1 = wz1*(p[8][k] - edgez1); + const vector corz2 = wz2*(p[9][k] - edgez2); + const vector corz3 = wz3*(p[10][k] - edgez3); + const vector corz4 = wz4*(p[11][k] - edgez4); - points_[vertexNo] += + points_[vijk] += + ( + corx1 + corx2 + corx3 + corx4 + + cory1 + cory2 + cory3 + cory4 + + corz1 + corz2 + corz3 + corz4 + ); + } + } + } + } + + if (!nCurvedFaces()) return; + + // Apply curvature correction to face points + FixedList facePoints(this->facePoints(points_)); + correctFacePoints(facePoints); + + // Loop over points and apply the correction from the from the i-faces + for (label ii=0; ii<=ni; ii++) + { + // Process the points on the curved faces last + label i = (ii + 1)%(ni + 1); + + for (label j=0; j<=nj; j++) + { + for (label k=0; k<=nk; k++) + { + // Skip non-curved faces and edges + if (flatFaceOrEdge(i, j, k)) continue; + + const label vijk = pointLabel(i, j, k); + + scalar wf0 = ( - corx1 + corx2 + corx3 + corx4 - + cory1 + cory2 + cory3 + cory4 - + corz1 + corz2 + corz3 + corz4 + (1 - w0)*(1 - w4)*(1 - w8) + + (1 - w1)*w4*(1 - w11) + + (1 - w2)*w7*w11 + + (1 - w3)*(1 - w7)*w8 + ); + + scalar wf1 = + ( + w0*(1 - w5)*(1 - w9) + + w1*w5*(1 - w10) + + w2*w5*w10 + + w3*(1 - w6)*w9 + ); + + const scalar sumWf = wf0 + wf1; + wf0 /= sumWf; + wf1 /= sumWf; + + points_[vijk] += + ( + wf0 + *( + facePoints[0][facePointLabel(0, j, k)] + - points_[pointLabel(0, j, k)] + ) + + wf1 + *( + facePoints[1][facePointLabel(1, j, k)] + - points_[pointLabel(ni, j, k)] + ) + ); + } + } + } + + // Loop over points and apply the correction from the from the j-faces + for (label jj=0; jj<=nj; jj++) + { + // Process the points on the curved faces last + label j = (jj + 1)%(nj + 1); + + for (label i=0; i<=ni; i++) + { + for (label k=0; k<=nk; k++) + { + // Skip flat faces and edges + if (flatFaceOrEdge(i, j, k)) continue; + + const label vijk = pointLabel(i, j, k); + + scalar wf2 = + ( + (1 - w4)*(1 - w1)*(1 - w8) + + (1 - w5)*w0*(1 - w9) + + (1 - w6)*w3*w9 + + (1 - w7)*(1 - w3)*w8 + ); + + scalar wf3 = + ( + w4*(1 - w1)*(1 - w11) + + w5*w1*(1 - w10) + + w6*w2*w10 + + w7*(1 - w2)*w11 + ); + + const scalar sumWf = wf2 + wf3; + wf2 /= sumWf; + wf3 /= sumWf; + + points_[vijk] += + ( + wf2 + *( + facePoints[2][facePointLabel(2, i, k)] + - points_[pointLabel(i, 0, k)] + ) + + wf3 + *( + facePoints[3][facePointLabel(3, i, k)] + - points_[pointLabel(i, nj, k)] + ) + ); + } + } + } + + // Loop over points and apply the correction from the from the k-faces + for (label kk=0; kk<=nk; kk++) + { + // Process the points on the curved faces last + label k = (kk + 1)%(nk + 1); + + for (label i=0; i<=ni; i++) + { + for (label j=0; j<=nj; j++) + { + // Skip flat faces and edges + if (flatFaceOrEdge(i, j, k)) continue; + + const label vijk = pointLabel(i, j, k); + + scalar wf4 = + ( + (1 - w8)*(1 - w0)*(1 - w4) + + (1 - w9)*w0*(1 - w5) + + (1 - w10)*w1*w5 + + (1 - w11)*(1 - w1)*w4 + ); + + scalar wf5 = + ( + w8*(1 - w3)*(1 - w7) + + w9*w3*(1 - w6) + + w10*w2*w6 + + w11*(1 - w2)*w7 + ); + + const scalar sumWf = wf4 + wf5; + wf4 /= sumWf; + wf5 /= sumWf; + + points_[vijk] += + ( + wf4 + *( + facePoints[4][facePointLabel(4, i, j)] + - points_[pointLabel(i, j, 0)] + ) + + wf5 + *( + facePoints[5][facePointLabel(5, i, j)] + - points_[pointLabel(i, j, nk)] + ) ); } } @@ -214,230 +346,164 @@ void Foam::block::createPoints() const } -void Foam::block::createCells() const +Foam::List> Foam::block::cells() const { const label ni = density().x(); const label nj = density().y(); const label nk = density().z(); - // - // Generate cells - // - cells_.clear(); - cells_.setSize(nCells()); + List> cells(nCells()); - label cellNo = 0; + label celli = 0; - for (label k = 0; k < nk; k++) + for (label k=0; k>, 6>& +Foam::block::boundaryPatches() const { - return *this; + return boundaryPatches_; } diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C index 2d6b308ec..7ccf24d42 100644 --- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C +++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C @@ -102,6 +102,32 @@ void Foam::blockDescriptor::check(const Istream& is) } +void Foam::blockDescriptor::findCurvedFaces() +{ + const faceList blockFaces(blockShape().faces()); + + forAll(blockFaces, blockFacei) + { + forAll(faces_, facei) + { + if + ( + face::sameVertices + ( + faces_[facei].vertices(), + blockFaces[blockFacei] + ) + ) + { + curvedFaces_[blockFacei] = facei; + nCurvedFaces_++; + break; + } + } + } +} + + // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // Foam::blockDescriptor::blockDescriptor @@ -109,6 +135,7 @@ Foam::blockDescriptor::blockDescriptor const cellShape& bshape, const pointField& vertices, const blockEdgeList& edges, + const blockFaceList& faces, const Vector