diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C index abc0c0ba03..416f06c1bf 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C +++ b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.C @@ -31,7 +31,7 @@ License template Foam::simpleMatrix::simpleMatrix(const label mSize) : - scalarSquareMatrix(mSize), + scalarSquareMatrix(mSize, mSize, pTraits::zero), source_(mSize, pTraits::zero) {} diff --git a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H b/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H index b854355003..e203972ffe 100644 --- a/src/OpenFOAM/matrices/simpleMatrix/simpleMatrix.H +++ b/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& source() { return source_; } + //- Return const-access to the source const Field& source() const { return source_; diff --git a/src/mesh/blockMesh/curvedEdges/BSpline.C b/src/mesh/blockMesh/curvedEdges/BSpline.C index 4e53399719..5aaf693c29 100644 --- a/src/mesh/blockMesh/curvedEdges/BSpline.C +++ b/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 M(NKnots+2, pTraits::zero); + simpleMatrix 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); } diff --git a/src/mesh/blockMesh/curvedEdges/arcEdge.C b/src/mesh/blockMesh/curvedEdges/arcEdge.C index bb7b3b7a8c..afb62a8f8b 100644 --- a/src/mesh/blockMesh/curvedEdges/arcEdge.C +++ b/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_; } diff --git a/src/mesh/blockMesh/curvedEdges/polySplineEdge.C b/src/mesh/blockMesh/curvedEdges/polySplineEdge.C index 841f6f986a..866689d8cf 100644 --- a/src/mesh/blockMesh/curvedEdges/polySplineEdge.C +++ b/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); - - label N = spl.nKnots(); - scalar init = 1.0/(N - 1); - scalar interval = (N - scalar(3))/N; + const label NKnots = spl.nKnots(); + const scalar init = 1.0/(NKnots - 1); + scalar interval = (NKnots - scalar(3.0))/NKnots; 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; } diff --git a/src/mesh/blockMesh/curvedEdges/spline.C b/src/mesh/blockMesh/curvedEdges/spline.C index 871c7f7db4..67b8c79559 100644 --- a/src/mesh/blockMesh/curvedEdges/spline.C +++ b/src/mesh/blockMesh/curvedEdges/spline.C @@ -26,17 +26,9 @@ License #include "spline.H" -// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // +// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * // -Foam::spline::spline(const pointField& knotPoints) -: - knots_(knotPoints) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -Foam::scalar Foam::spline::B(const scalar& tau) const +Foam::scalar Foam::spline::B(const scalar& tau) { 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