Commit d34aee63 authored by Andrew Heather's avatar Andrew Heather
Browse files

Merge branch 'olesenm'

parents 3c54d945 7a321282
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2009 OpenCFD Ltd.
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
OpenFOAM is free software; you can redistribute it and/or modify it
under the terms of the GNU General Public License as published by the
Free Software Foundation; either version 2 of the License, or (at your
option) any later version.
OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
for more details.
You should have received a copy of the GNU General Public License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
\*---------------------------------------------------------------------------*/
#include "argList.H"
#include "vector.H"
#include "IFstream.H"
#include "BSpline.H"
using namespace Foam;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:
int main(int argc, char *argv[])
{
argList::noParallel();
argList::validArgs.insert("file .. fileN");
argList args(argc, argv, false, true);
forAll(args.additionalArgs(), argI)
{
const string& srcFile = args.additionalArgs()[argI];
Info<< nl << "reading " << srcFile << nl;
IFstream ifs(srcFile);
List<pointField> splinePointFields(ifs);
forAll(splinePointFields, splineI)
{
Info<<"convert " << splinePointFields[splineI] << " to bspline" << endl;
BSpline spl(splinePointFields[splineI], vector::zero, vector::zero);
Info<< "1/2 = " << spl.position(0.5) << endl;
}
}
return 0;
}
// ************************************************************************* //
BSplineTest.C
EXE = $(FOAM_USER_APPBIN)/BSplineTest
EXE_INC = \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/mesh/blockMesh/lnInclude
EXE_LIBS = -lblockMesh
/*--------------------------------*- C++ -*----------------------------------*\
| ========= | |
| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| \\ / O peration | Version: 1.6 |
| \\ / A nd | Web: http://www.OpenFOAM.org |
| \\/ M anipulation | |
\*---------------------------------------------------------------------------*/
(
// Upper body longitudinal splines.
(
(-0.22685 -0.01125166 0) // 7
(-0.21685 -0.01340204 0)
(-0.20685 -0.01529684 0)
(-0.19685 -0.01694748 0)
(-0.18685 -0.01836538 0)
(-0.17685 -0.01956197 0)
(-0.16685 -0.02054868 0)
(-0.15685 -0.02133693 0)
(-0.14685 -0.02193816 0)
(-0.13685 -0.02236377 0)
(-0.12685 -0.02262521 0)
(-0.11685 -0.02273389 0) // 2
)
(
(-0.22685 0 0.01125166) // 8
(-0.21685 0 0.01340204)
(-0.20685 0 0.01529684)
(-0.19685 0 0.01694748)
(-0.18685 0 0.01836538)
(-0.17685 0 0.01956197)
(-0.16685 0 0.02054868)
(-0.15685 0 0.02133693)
(-0.14685 0 0.02193816)
(-0.13685 0 0.02236377)
(-0.12685 0 0.02262521)
(-0.11685 0 0.02273389) // 3
)
(
(-0.22685 0.01125166 0) // 9
(-0.21685 0.01340204 0)
(-0.20685 0.01529684 0)
(-0.19685 0.01694748 0)
(-0.18685 0.01836538 0)
(-0.17685 0.01956197 0)
(-0.16685 0.02054868 0)
(-0.15685 0.02133693 0)
(-0.14685 0.02193816 0)
(-0.13685 0.02236377 0)
(-0.12685 0.02262521 0)
(-0.11685 0.02273389 0) // 4
)
(
(-0.22685 0 -0.01125166) // 6
(-0.21685 0 -0.01340204)
(-0.20685 0 -0.01529684)
(-0.19685 0 -0.01694748)
(-0.18685 0 -0.01836538)
(-0.17685 0 -0.01956197)
(-0.16685 0 -0.02054868)
(-0.15685 0 -0.02133693)
(-0.14685 0 -0.02193816)
(-0.13685 0 -0.02236377)
(-0.12685 0 -0.02262521)
(-0.11685 0 -0.02273389) // 1
)
);
......@@ -22,7 +22,7 @@ License
along with OpenFOAM; if not, write to the Free Software Foundation,
Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
Namespace
InNamespace
Foam
Description
......@@ -35,8 +35,6 @@ Description
#include "mathematicalConstants.H"
using namespace Foam::constant::mathematical;
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
......@@ -47,13 +45,13 @@ namespace Foam
//- Conversion from degrees to radians
inline scalar degToRad(const scalar deg)
{
return (deg*pi/180.0);
return (deg * constant::mathematical::pi/180.0);
}
//- Conversion from radians to degrees
inline scalar radToDeg(const scalar rad)
{
return (rad*180.0/pi);
return (rad * 180.0/constant::mathematical::pi);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -78,8 +78,7 @@ public:
// Note: this does not initialise the coefficients or the source.
simpleMatrix(const label);
//- Construct given size and initial values for the
// coefficients and source
//- Construct given size and initial values for coefficients and source
simpleMatrix(const label, const scalar, const Type&);
//- Construct from components
......
......@@ -55,9 +55,9 @@ namespace Foam
class ifEqEqOp
{
public:
void operator()(label x, const label y) const
void operator()(label& x, const label y) const
{
x = (x==y) ? x : value;
x = (x == y) ? x : value;
}
};
}
......
......@@ -38,74 +38,70 @@ Foam::pointField Foam::BSpline::findKnots
const vector& sndend
)
{
label newnKnots(allknots.size() + 2);
label NKnots(allknots.size());
pointField newknots(newnKnots);
const label NKnots = allknots.size();
// set up 1/6 and 2/3 which are the matrix elements throughout most
// of the matrix
register scalar oneSixth = 1.0/6.0;
register scalar twoThird = 2.0/3.0;
register const scalar oneSixth = 1.0/6.0;
register const scalar twoThird = 2.0/3.0;
simpleMatrix<vector> M(newnKnots, 0, vector::zero);
simpleMatrix<vector> M(NKnots+2, 0, vector::zero);
// set up the matrix
M[0][0] = -0.5*scalar(NKnots - 1);
M[0][2] = 0.5*scalar(NKnots - 1);
for (register label i=1; i<newnKnots-1; i++)
for (register label i = 1; i <= NKnots; i++)
{
M[i][i-1] = oneSixth;
M[i][i] = twoThird;
M[i][i+1] = oneSixth;
}
M[newnKnots - 1][newnKnots - 3] = -0.5*scalar(NKnots - 1);
M[newnKnots - 1][newnKnots - 1] = 0.5*scalar(NKnots - 1);
M[NKnots+1][NKnots-1] = -0.5*scalar(NKnots - 1);
M[NKnots+1][NKnots+1] = 0.5*scalar(NKnots - 1);
// set up the vector
for (label i=1; i<=NKnots; i++)
for (register label i = 1; i <= NKnots; i++)
{
M.source()[i] = allknots[i-1];
}
// set the gradients at the two ends
// set the gradients at the ends:
if (mag(fstend)<1e-8)
if (mag(fstend) < 1e-8)
{
// set to the default : forward differences on the end knots
// default : forward differences on the end knots
M.source()[0] = allknots[1] - allknots[0];
M.source()[0] /= mag(M.source()[0]);
}
else
{
// use the gradient vector provided
M.source()[0] = fstend/mag(fstend);
}
if (mag(sndend)<1e-8)
{
// default : forward differences on the end knots
M.source()[NKnots+1] = M.source()[NKnots-1] - M.source()[NKnots];
M.source()[NKnots+1] /= mag(M.source()[NKnots+1]);
}
else
{
// set to the gradient vectors provided
M.source()[0] = fstend/mag(fstend);
// use the gradient vector provided
M.source()[NKnots+1] = sndend/mag(sndend);
}
// invert the equation to find the control knots
newknots = M.solve();
return newknots;
// invert the equation to find the control knots
return M.solve();
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::BSpline::BSpline(const pointField& Knots)
:
spline(findKnots(Knots))
{}
Foam::BSpline::BSpline
(
const pointField& Knots,
......
......@@ -56,8 +56,8 @@ class BSpline
pointField findKnots
(
const pointField&,
const vector& fstend = vector::zero,
const vector& sndend = vector::zero
const vector& fstend,
const vector& sndend
);
//- Disallow default bitwise copy construct
......@@ -71,15 +71,12 @@ public:
// Constructors
//- Construct from components
BSpline(const pointField& knots);
//- Construct from components
BSpline
(
const pointField& knots,
const vector& fstend,
const vector& sndend
const vector& fstend = vector::zero,
const vector& sndend = vector::zero
);
......
......@@ -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_;
}
......
......@@ -36,20 +36,6 @@ namespace Foam
addToRunTimeSelectionTable(curvedEdge, polySplineEdge, Istream);
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
//! @cond fileScope
inline label nsize(const label otherKnotsSize, const label nBetweenKnots)
{
return otherKnotsSize*(1 + nBetweenKnots) + nBetweenKnots + 2;
}
//! @endcond fileScope
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
// intervening : returns a list of the points making up the polyLineEdge
......@@ -68,18 +54,25 @@ Foam::pointField Foam::polySplineEdge::intervening
const vector& sndend
)
{
BSpline spl(knotlist(points_, start_, end_, otherknots), fstend, sndend);
label nSize(nsize(otherknots.size(), nBetweenKnots));
BSpline spl
(
knotlist(points_, start_, end_, otherknots),
fstend,
sndend
);
pointField ans(nSize);
const label nSize
(
otherknots.size() * (1 + nBetweenKnots) + nBetweenKnots + 2
);
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 +128,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;
}
......
......@@ -26,17 +26,9 @@ License
#include "spline.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::spline::spline(const pointField& knotPoints)
:
knots_(knotPoints)
{}
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
// * * * * * * * * * * * * * * * 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<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;
......
......@@ -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&);
......
......@@ -63,11 +63,8 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
const scalar dt
) const
{
scalar pf, cf, pr, cr;
label lRef, rRef;
label nSpecie = this->model_.nSpecie();
simpleMatrix<scalar> RR(nSpecie);
const label nSpecie = this->model_.nSpecie();
simpleMatrix<scalar> RR(nSpecie, 0, 0);
for (label i=0; i<nSpecie; i++)
{
......@@ -83,6 +80,9 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
{
const Reaction<ThermoType>& R = this->model_.reactions()[i];
scalar pf, cf, pr, cr;
label lRef, rRef;
scalar omegai = this->model_.omega
(
R, c, T, p, pf, cf, lRef, pr, cr, rRef
......@@ -132,7 +132,7 @@ Foam::scalar Foam::EulerImplicit<CompType, ThermoType>::solve
// estimate the next time step
scalar tMin = GREAT;
label nEqns = this->model_.nEqns();
const label nEqns = this->model_.nEqns();
scalarField c1(nEqns, 0.0);
for (label i=0; i<nSpecie; i++)
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment