simpleMatrix initialize with coefficients with 0.0

- the element of least surprise

src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C

@@ -31,7 +31,7 @@ License
 template<class Type>
 Foam::simpleMatrix<Type>::simpleMatrix(const label mSize)
 :
-    scalarSquareMatrix(mSize),
+    scalarSquareMatrix(mSize, mSize, pTraits<scalar>::zero),
     source_(mSize, pTraits<Type>::zero)
 {}
 

src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H

@@ -26,7 +26,7 @@ Class
     Foam::simpleMatrix
 
 Description
-    Foam::simpleMatrix
+    A simple square matrix solver with scalar coefficients.
 
 SourceFiles
     simpleMatrix.C
@@ -74,10 +74,10 @@ public:
 
     // Constructors
 
-        //- Construct given size
+        //- Construct given size, initializing matrix coefficients to zero
         simpleMatrix(const label);
 
-        //- Construct given size and initial value
+        //- Construct given size, and initial value for matrix coefficients
         simpleMatrix(const label, const scalar&);
 
         //- Construct from components
@@ -94,11 +94,13 @@ public:
 
         // Access
 
+            //- Return access to the source
             Field<Type>& source()
             {
                 return source_;
             }
 
+            //- Return const-access to the source
             const Field<Type>& source() const
             {
                 return source_;

src/mesh/blockMesh/curvedEdges/BSpline.C

@@ -46,7 +46,7 @@ Foam::pointField Foam::BSpline::findKnots
     register const scalar oneSixth = 1.0/6.0;
     register const scalar twoThird = 2.0/3.0;
 
-    simpleMatrix<vector> M(NKnots+2, pTraits<scalar>::zero);
+    simpleMatrix<vector> M(NKnots+2);
 
     // set up the matrix
     M[0][0] = -0.5*scalar(NKnots - 1);
@@ -68,9 +68,9 @@ Foam::pointField Foam::BSpline::findKnots
         M.source()[i] = allknots[i-1];
     }
 
+    // set the gradients at the ends:
 
-    // set the gradients at the two ends
+    if (mag(fstend) < 1e-8)
-    if (mag(fstend)<1e-8)
     {
         // default : forward differences on the end knots
         M.source()[0] = allknots[1] - allknots[0];
@@ -78,7 +78,7 @@ Foam::pointField Foam::BSpline::findKnots
     }
     else
     {
-        // set to the gradient vector provided
+        // use the gradient vector provided
         M.source()[0] = fstend/mag(fstend);
     }
 
@@ -90,7 +90,7 @@ Foam::pointField Foam::BSpline::findKnots
     }
     else
     {
-        // set to the gradient vector provided
+        // use the gradient vector provided
         M.source()[NKnots+1] = sndend/mag(sndend);
     }
 

src/mesh/blockMesh/curvedEdges/arcEdge.C

@@ -75,7 +75,7 @@ Foam::cylindricalCS Foam::arcEdge::calcAngle()
     angle_ = radToDeg(acos(tmp));
 
     // check if the vectors define an exterior or an interior arcEdge
-    if (((r1  ^ r2)&(r1 ^ r3)) < 0.0)
+    if (((r1 ^ r2)&(r1 ^ r3)) < 0.0)
     {
         angle_ = 360.0 - angle_;
     }

src/mesh/blockMesh/curvedEdges/polySplineEdge.C

@@ -68,18 +68,22 @@ Foam::pointField Foam::polySplineEdge::intervening
     const vector& sndend
 )
 {
-    BSpline spl(knotlist(points_, start_, end_, otherknots), fstend, sndend);
+    BSpline spl
+    (
+        knotlist(points_, start_, end_, otherknots),
+        fstend,
+        sndend
+    );
 
-    label nSize(nsize(otherknots.size(), nBetweenKnots));
+    const label nSize(nsize(otherknots.size(), nBetweenKnots));
 
-    pointField ans(nSize);
+    const label NKnots = spl.nKnots();
-
+    const scalar init = 1.0/(NKnots - 1);
-    label N = spl.nKnots();
+    scalar interval = (NKnots - scalar(3.0))/NKnots;
-    scalar init = 1.0/(N - 1);
-    scalar interval = (N - scalar(3))/N;
     interval /= otherknots.size() + 1;
     interval /= nBetweenKnots + 1;
 
+    pointField ans(nSize);
     ans[0] = points_[start_];
 
     register scalar index(init);
@@ -135,17 +139,8 @@ Foam::polySplineEdge::polySplineEdge
     vector fstend(is);
     vector sndend(is);
 
-    controlPoints_.setSize(nsize(otherKnots_.size(), nInterKnots));
-    // why does this need to be here (to avoid a crash)?
-    // 'intervening' uses BSpline to solve the new points
-    // it seems to be going badly there
-    distances_.setSize(controlPoints_.size());
-
     controlPoints_ = intervening(otherKnots_, nInterKnots, fstend, sndend);
     calcDistances();
-
-    // Info<< "polyLine[" << start_ << " " << end_
-    //     << "] controlPoints " << controlPoints_ << endl;
 }
 
 

src/mesh/blockMesh/curvedEdges/spline.C

@@ -26,17 +26,9 @@ License
 
 #include "spline.H"
 
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
 
-Foam::spline::spline(const pointField& knotPoints)
+Foam::scalar Foam::spline::B(const scalar& tau)
-:
-    knots_(knotPoints)
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-Foam::scalar Foam::spline::B(const scalar& tau) const
 {
     if (tau <= -2.0 || tau >= 2.0)
     {
@@ -44,7 +36,7 @@ Foam::scalar Foam::spline::B(const scalar& tau) const
     }
     else if (tau <= -1.0)
     {
-        return pow((2.0 + tau),3.0)/6.0;
+        return pow((2.0 + tau), 3.0)/6.0;
     }
     else if (tau <= 0.0)
     {
@@ -56,7 +48,7 @@ Foam::scalar Foam::spline::B(const scalar& tau) const
     }
     else if (tau <= 2.0)
     {
-        return pow((2.0 - tau),3.0)/6.0;
+        return pow((2.0 - tau), 3.0)/6.0;
     }
     else
     {
@@ -70,13 +62,23 @@ Foam::scalar Foam::spline::B(const scalar& tau) const
 }
 
 
-Foam::vector Foam::spline::position(const scalar mu1) const
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::spline::spline(const pointField& knotPoints)
+:
+    knots_(knotPoints)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::vector Foam::spline::position(const scalar mu) const
 {
     vector loc(vector::zero);
 
-    for (register label i=0; i<knots_.size(); i++)
+    for (register label i=0; i < knots_.size(); i++)
     {
-        loc += B((knots_.size() - 1)*mu1 - i)*knots_[i];
+        loc += B((knots_.size() - 1)*mu - i)*knots_[i];
     }
 
     return loc;

src/mesh/blockMesh/curvedEdges/spline.H

@@ -27,6 +27,7 @@ Class
 
 Description
     Define a basic spline on nKnots knots.
+
     The spline does not go anywhere near these knots
     (will act as a base type for various splines that will have real uses)
 
@@ -59,7 +60,7 @@ class spline
     // Private Member Functions
 
         //- Blending function for constructing spline
-        scalar B(const scalar&) const;
+        static scalar B(const scalar&);
 
         //- Disallow default bitwise copy construct
         spline(const spline&);