From e53520ecbfddee15d71885ab483d07f883b2bc53 Mon Sep 17 00:00:00 2001
From: Henry Weller <http://cfd.direct>
Date: Sun, 9 Oct 2016 15:16:41 +0100
Subject: [PATCH] blockMesh::blockDescriptor: Generate the edge points and
 weights on demand

---
 src/mesh/blockMesh/block/blockCreate.C        | 144 ++++++++----------
 .../blockDescriptor/blockDescriptor.C         |  10 --
 .../blockDescriptor/blockDescriptor.H         |  34 ++---
 .../blockDescriptor/blockDescriptorEdges.C    |  97 ++++++------
 .../blockDescriptor/blockDescriptorI.H        |  14 --
 src/mesh/blockMesh/blockMesh/blockMesh.C      |   2 +-
 6 files changed, 131 insertions(+), 170 deletions(-)

diff --git a/src/mesh/blockMesh/block/blockCreate.C b/src/mesh/blockMesh/block/blockCreate.C
index 27def57ecd..0aac98f197 100644
--- a/src/mesh/blockMesh/block/blockCreate.C
+++ b/src/mesh/blockMesh/block/blockCreate.C
@@ -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] +=
                 (
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
index 2433f94489..2d6b308ecd 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C
@@ -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();
 }
 
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
index f2bd8d2af0..0a95cacd07 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.H
@@ -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
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
index 5840450d7b..55895c34a8 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorEdges.C
@@ -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);
 }
 
 
diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
index a4a81d1ffa..37439b1a0b 100644
--- a/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
+++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptorI.H
@@ -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_;
diff --git a/src/mesh/blockMesh/blockMesh/blockMesh.C b/src/mesh/blockMesh/blockMesh/blockMesh.C
index cceed4ebfa..90f75c8014 100644
--- a/src/mesh/blockMesh/blockMesh/blockMesh.C
+++ b/src/mesh/blockMesh/blockMesh/blockMesh.C
@@ -48,7 +48,7 @@ Foam::blockMesh::blockMesh(const IOdictionary& dict, const word& regionName)
 
     if (fastMerge)
     {
-         calcMergeInfoFast();
+        calcMergeInfoFast();
     }
     else
     {