zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

Merge branch 'master' of ssh://noisy/home/noisy3/OpenFOAM/OpenFOAM-dev

Browse Source
This commit is contained in:
henry
2009-04-03 17:46:32 +01:00
parent c09f9589b7 b2ecadd9f5
commit 150ad8c97a
6 changed files with 222 additions and 80 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
applications/utilities/mesh/advanced/refineHexMesh/refineHexMesh.C
View File

@ -182,6 +182,7 @@ int main(int argc, char *argv[])
if (overwrite)
{
mesh.setInstance(oldInstance);
meshCutter.setInstance(oldInstance);
}
Info<< "Writing mesh to " << runTime.timeName() << endl;

6
src/ODE/ODESolvers/SIBS/SIBS.C
View File

@ -71,8 +71,8 @@ void Foam::SIBS::solve
const ODE& ode,
scalar& x,
scalarField& y,
scalarField& dydx,
const scalar eps,
scalarField& dydx,
const scalar eps,
const scalarField& yScale,
const scalar hTry,
scalar& hDid,
@ -96,7 +96,7 @@ void Foam::SIBS::solve
{
for (register label k=0; k<iq; k++)
{
alpha_[k][iq] =
alpha_[k][iq] =
pow(eps1, (a_[k + 1] - a_[iq + 1])
/((a_[iq + 1] - a_[0] + 1.0)*(2*k + 3)));
}

51
src/ODE/ODESolvers/SIBS/SIBS.H
View File

@ -29,8 +29,8 @@ Description
Foam::SIBS
SourceFiles
SIBSCK.C
SIBSQS.C
SIMPR.C
polyExtrapolate.C
\*---------------------------------------------------------------------------*/
@ -61,7 +61,7 @@ class SIBS
mutable scalarField a_;
mutable scalarSquareMatrix alpha_;
mutable scalarSquareMatrix d_p_;
mutable scalarRectangularMatrix d_p_;
mutable scalarField x_p_;
mutable scalarField err_;
@ -75,30 +75,31 @@ class SIBS
mutable scalar epsOld_, xNew_;
void SIMPR
(
const ODE& ode,
const scalar xStart,
const scalarField& y,
const scalarField& dydx,
const scalarField& dfdx,
const scalarSquareMatrix& dfdy,
const scalar deltaX,
const label nSteps,
scalarField& yEnd
) const;
// Private member functions
void SIMPR
(
const ODE& ode,
const scalar xStart,
const scalarField& y,
const scalarField& dydx,
const scalarField& dfdx,
const scalarSquareMatrix& dfdy,
const scalar deltaX,
const label nSteps,
scalarField& yEnd
) const;
void polyExtrapolate
(
const label iest,
const scalar xest,
const scalarField& yest,
scalarField& yz,
scalarField& dy,
scalarField& x_p,
scalarSquareMatrix& d_p
) const;
void polyExtrapolate
(
const label iest,
const scalar xest,
const scalarField& yest,
scalarField& yz,
scalarField& dy,
scalarField& x_p,
scalarRectangularMatrix& d_p
) const;
public:

2
src/ODE/ODESolvers/SIBS/polyExtrapolate.C
View File

@ -36,7 +36,7 @@ void Foam::SIBS::polyExtrapolate
scalarField& yz,
scalarField& dy,
scalarField& x,
scalarSquareMatrix& d
scalarRectangularMatrix& d
) const
{
label n = yz.size();

240
src/meshTools/searchableSurface/searchableCylinder.C
View File

@ -105,6 +105,13 @@ Foam::pointIndexHit Foam::searchableCylinder::findNearest
}
Foam::scalar Foam::searchableCylinder::radius2(const point& pt) const
{
const vector x = (pt-point1_) ^ unitDir_;
return x&x;
}
// From http://www.gamedev.net/community/forums/topic.asp?topic_id=467789 -
// intersection of cylinder with ray
void Foam::searchableCylinder::findLineAll
@ -118,13 +125,91 @@ void Foam::searchableCylinder::findLineAll
near.setMiss();
far.setMiss();
// Line as P = start + t*V
const vector V(end-start);
vector point1Start(start-point1_);
vector point2Start(start-point2_);
vector point1End(end-point1_);
//Pout<< "point1:" << point1_ << " point2:" << point2_
// << " start:" << start << " end:" << end << endl;
// Quick rejection of complete vector outside endcaps
scalar s1 = point1Start&unitDir_;
scalar s2 = point1End&unitDir_;
const vector x = (start-point1_) ^ unitDir_;
if ((s1 < 0 && s2 < 0) || (s1 > magDir_ && s2 > magDir_))
{
return;
}
// Line as P = start+t*V where V is unit vector and t=[0..mag(end-start)]
vector V(end-start);
scalar magV = mag(V);
if (magV < ROOTVSMALL)
{
return;
}
V /= magV;
// We now get the nearest intersections to start. This can either be
// the intersection with the end plane or with the cylinder side.
// Get the two points (expressed in t) on the end planes. This is to
// clip any cylinder intersection against.
scalar tPoint1;
scalar tPoint2;
// Maintain the two intersections with the endcaps
scalar tNear = VGREAT;
scalar tFar = VGREAT;
{
scalar s = (V&unitDir_);
if (mag(s) > VSMALL)
{
tPoint1 = -s1/s;
tPoint2 = -(point2Start&unitDir_)/s;
if (tPoint2 < tPoint1)
{
Swap(tPoint1, tPoint2);
}
if (tPoint1 > magV || tPoint2 < 0)
{
return;
}
// See if the points on the endcaps are actually inside the cylinder
if (tPoint1 >= 0 && tPoint1 <= magV)
{
if (radius2(start+tPoint1*V) <= sqr(radius_))
{
tNear = tPoint1;
}
}
if (tPoint2 >= 0 && tPoint2 <= magV)
{
if (radius2(start+tPoint2*V) <= sqr(radius_))
{
// Check if already have a near hit from point1
if (tNear <= 1)
{
tFar = tPoint2;
}
else
{
tNear = tPoint2;
}
}
}
}
else
{
// Vector perpendicular to cylinder. Check for outside already done
// above so just set tpoint to allow all.
tPoint1 = -VGREAT;
tPoint2 = VGREAT;
}
}
const vector x = point1Start ^ unitDir_;
const vector y = V ^ unitDir_;
const scalar d = sqr(radius_);
@ -135,11 +220,12 @@ void Foam::searchableCylinder::findLineAll
const scalar disc = b*b-4*a*c;
//Pout<< "a:" << a << " b:" << b << " c:" << c << " disc:" << disc
// << endl;
scalar t1 = -VGREAT;
scalar t2 = VGREAT;
if (disc < 0)
{
// Fully outside
return;
}
else if (disc < ROOTVSMALL)
@ -147,12 +233,40 @@ void Foam::searchableCylinder::findLineAll
// Single solution
if (mag(a) > ROOTVSMALL)
{
scalar t = -b/(2*a);
if (t >= 0 && t <= 1)
t1 = -b/(2*a);
//Pout<< "single solution t:" << t1
// << " for start:" << start << " end:" << end
// << " c:" << c << endl;
if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
{
near.setPoint(start + t*V);
near.setHit();
near.setIndex(0);
// valid. Insert sorted.
if (t1 < tNear)
{
tFar = tNear;
tNear = t1;
}
else if (t1 < tFar)
{
tFar = t1;
}
}
else
{
return;
}
}
else
{
// Aligned with axis. Check if outside radius
//Pout<< "small discriminant:" << disc
// << " for start:" << start << " end:" << end
// << " magV:" << magV
// << " c:" << c << endl;
if (c > 0)
{
return;
}
}
}
@ -162,41 +276,79 @@ void Foam::searchableCylinder::findLineAll
{
scalar sqrtDisc = sqrt(disc);
scalar t1 = (-b + sqrtDisc)/2*a;
scalar t2 = (-b - sqrtDisc)/2*a;
if (t1 < t2)
t1 = (-b - sqrtDisc)/(2*a);
t2 = (-b + sqrtDisc)/(2*a);
if (t2 < t1)
{
if (t1 >= 0 && t1 <= 1)
Swap(t1, t2);
}
if (t1 >= 0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
{
// valid. Insert sorted.
if (t1 < tNear)
{
near.setPoint(start + t1*V);
near.setHit();
near.setIndex(0);
tFar = tNear;
tNear = t1;
}
if (t2 >= 0 && t2 <= 1)
else if (t1 < tFar)
{
far.setPoint(start + t2*V);
far.setHit();
far.setIndex(0);
tFar = t1;
}
}
else
if (t2 >= 0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
{
if (t2 >= 0 && t2 <= 1)
// valid. Insert sorted.
if (t2 < tNear)
{
near.setPoint(start + t2*V);
near.setHit();
near.setIndex(0);
tFar = tNear;
tNear = t2;
}
if (t1 >= 0 && t1 <= 1)
else if (t2 < tFar)
{
far.setPoint(start + t1*V);
far.setHit();
far.setIndex(0);
tFar = t2;
}
}
//Pout<< "two solutions t1:" << t1 << " t2:" << t2
// << " for start:" << start << " end:" << end
// << " magV:" << magV
// << " c:" << c << endl;
}
else
{
// Aligned with axis. Check if outside radius
//Pout<< "large discriminant:" << disc
// << " small a:" << a
// << " for start:" << start << " end:" << end
// << " magV:" << magV
// << " c:" << c << endl;
if (c > 0)
{
return;
}
}
}
// Check tNear, tFar
if (tNear >= 0 && tNear <= magV)
{
near.setPoint(start+tNear*V);
near.setHit();
near.setIndex(0);
if (tFar <= magV)
{
far.setPoint(start+tFar*V);
far.setHit();
far.setIndex(0);
}
}
else if (tFar >= 0 && tFar <= magV)
{
near.setPoint(start+tFar*V);
near.setHit();
near.setIndex(0);
}
}
@ -216,13 +368,7 @@ Foam::searchableCylinder::searchableCylinder
magDir_(mag(point2_-point1_)),
unitDir_((point2_-point1_)/magDir_),
radius_(radius)
{
Pout<< "point1_:" << point1_ << endl;
Pout<< "point2_:" << point2_ << endl;
Pout<< "magDir_:" << magDir_ << endl;
Pout<< "unitDir_:" << unitDir_ << endl;
Pout<< "radius_:" << radius_ << endl;
}
{}
Foam::searchableCylinder::searchableCylinder
@ -237,13 +383,7 @@ Foam::searchableCylinder::searchableCylinder
magDir_(mag(point2_-point1_)),
unitDir_((point2_-point1_)/magDir_),
radius_(readScalar(dict.lookup("radius")))
{
Pout<< "point1_:" << point1_ << endl;
Pout<< "point2_:" << point2_ << endl;
Pout<< "magDir_:" << magDir_ << endl;
Pout<< "unitDir_:" << unitDir_ << endl;
Pout<< "radius_:" << radius_ << endl;
}
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
@ -290,11 +430,10 @@ void Foam::searchableCylinder::findLine
{
info.setSize(start.size());
pointIndexHit b;
forAll(start, i)
{
// Pick nearest intersection. If none intersected take second one.
pointIndexHit b;
findLineAll(start[i], end[i], info[i], b);
if (!info[i].hit() && b.hit())
{
@ -313,11 +452,10 @@ void Foam::searchableCylinder::findLineAny
{
info.setSize(start.size());
pointIndexHit b;
forAll(start, i)
{
// Discard far intersection
pointIndexHit b;
findLineAll(start[i], end[i], info[i], b);
if (!info[i].hit() && b.hit())
{

2
src/meshTools/searchableSurface/searchableCylinder.H
View File

@ -86,6 +86,8 @@ private:
const scalar nearestDistSqr
) const;
scalar radius2(const point& pt) const;
//- Find intersection with cylinder
void findLineAll
(