new searchableSurfaces

This commit is contained in:
mattijs
2008-09-15 12:14:51 +01:00
parent 907e195c72
commit b77b0e5b39
4 changed files with 1140 additions and 0 deletions

View File

@ -0,0 +1,366 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "searchablePlate.H"
#include "addToRunTimeSelectionTable.H"
#include "SortableList.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchablePlate, 0);
addToRunTimeSelectionTable(searchableSurface, searchablePlate, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::direction Foam::searchablePlate::calcNormal(const point& span)
{
direction normalDir = 3;
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (span[dir] < 0)
{
FatalErrorIn("searchablePlate::calcNormal()")
<< "Span should have two positive and one zero entry. Now:"
<< span << exit(FatalError);
}
else if (span[dir] < VSMALL)
{
if (normalDir == 3)
{
normalDir = dir;
}
else
{
// Multiple zero entries. Flag and exit.
normalDir = 3;
break;
}
}
}
if (normalDir == 3)
{
FatalErrorIn("searchablePlate::calcNormal()")
<< "Span should have one and only zero entry. Now:" << span
<< exit(FatalError);
}
return normalDir;
}
// Returns miss or hit with face (always 0)
Foam::pointIndexHit Foam::searchablePlate::findNearest
(
const point& sample,
const scalar nearestDistSqr
) const
{
// For every component direction can be
// left of min, right of max or inbetween.
// - outside points: project first one x plane (either min().x()
// or max().x()), then onto y plane and finally z. You should be left
// with intersection point
// - inside point: find nearest side (compare to mid point). Project onto
// that.
// Project point on plane.
pointIndexHit info(true, sample, 0);
info.rawPoint()[normalDir_] = origin_[normalDir_];
// Clip to edges if outside
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (dir != normalDir_)
{
if (info.rawPoint()[dir] < origin_[dir])
{
info.rawPoint()[dir] = origin_[dir];
}
else if (info.rawPoint()[dir] > origin_[dir]+span_[dir])
{
info.rawPoint()[dir] = origin_[dir]+span_[dir];
}
}
}
// Check if outside. Optimisation: could do some checks on distance already
// on components above
if (magSqr(info.rawPoint() - sample) > nearestDistSqr)
{
info.setMiss();
info.setIndex(-1);
}
return info;
}
Foam::pointIndexHit Foam::searchablePlate::findLine
(
const point& start,
const point& end
) const
{
pointIndexHit info
(
true,
vector::zero,
0
);
const vector dir(end-start);
if (mag(dir[normalDir_]) < VSMALL)
{
info.setMiss();
info.setIndex(-1);
}
else
{
scalar t = (origin_[normalDir_]-start[normalDir_]) / dir[normalDir_];
if (t < 0 || t > 1)
{
info.setMiss();
info.setIndex(-1);
}
else
{
info.rawPoint() = start+t*dir;
info.rawPoint()[normalDir_] = origin_[normalDir_];
// Clip to edges
for (direction dir = 0; dir < vector::nComponents; dir++)
{
if (dir != normalDir_)
{
if (info.rawPoint()[dir] < origin_[dir])
{
info.setMiss();
info.setIndex(-1);
break;
}
else if (info.rawPoint()[dir] > origin_[dir]+span_[dir])
{
info.setMiss();
info.setIndex(-1);
break;
}
}
}
}
}
// Debug
if (info.hit())
{
treeBoundBox bb(origin_, origin_+span_);
bb.min()[normalDir_] -= 1E-6;
bb.max()[normalDir_] += 1E-6;
if (!bb.contains(info.hitPoint()))
{
FatalErrorIn("searchablePlate::findLine(..)")
<< "bb:" << bb << endl
<< "origin_:" << origin_ << endl
<< "span_:" << span_ << endl
<< "normalDir_:" << normalDir_ << endl
<< "hitPoint:" << info.hitPoint()
<< abort(FatalError);
}
}
return info;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchablePlate::searchablePlate
(
const IOobject& io,
const point& origin,
const vector& span
)
:
searchableSurface(io),
origin_(origin),
span_(span),
normalDir_(calcNormal(span_))
{}
Foam::searchablePlate::searchablePlate
(
const IOobject& io,
const dictionary& dict
)
:
searchableSurface(io),
origin_(dict.lookup("origin")),
span_(dict.lookup("span")),
normalDir_(calcNormal(span_))
{}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchablePlate::~searchablePlate()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
const Foam::wordList& Foam::searchablePlate::regions() const
{
if (regions_.size() == 0)
{
regions_.setSize(1);
regions_[0] = "region0";
}
return regions_;
}
void Foam::searchablePlate::findNearest
(
const pointField& samples,
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const
{
info.setSize(samples.size());
forAll(samples, i)
{
info[i] = findNearest(samples[i], nearestDistSqr[i]);
}
}
void Foam::searchablePlate::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
info.setSize(start.size());
forAll(start, i)
{
info[i] = findLine(start[i], end[i]);
}
}
void Foam::searchablePlate::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
findLine(start, end, info);
}
void Foam::searchablePlate::findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
List<pointIndexHit> nearestInfo;
findLine(start, end, nearestInfo);
info.setSize(start.size());
forAll(info, pointI)
{
if (nearestInfo[pointI].hit())
{
info[pointI].setSize(1);
info[pointI][0] = nearestInfo[pointI];
}
else
{
info[pointI].clear();
}
}
}
void Foam::searchablePlate::getRegion
(
const List<pointIndexHit>& info,
labelList& region
) const
{
region.setSize(info.size());
region = 0;
}
void Foam::searchablePlate::getNormal
(
const List<pointIndexHit>& info,
vectorField& normal
) const
{
normal.setSize(info.size());
normal = vector::zero;
forAll(normal, i)
{
normal[i][normalDir_] = 1.0;
}
}
void Foam::searchablePlate::getVolumeType
(
const pointField& points,
List<volumeType>& volType
) const
{
FatalErrorIn
(
"searchableCollection::getVolumeType(const pointField&"
", List<volumeType>&) const"
) << "Volume type not supported for plate."
<< exit(FatalError);
}
// ************************************************************************* //

View File

@ -0,0 +1,215 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
Class
Foam::searchablePlate
Description
Searching on finite plate. Plate has to be aligned with coordinate
axes!
SourceFiles
searchablePlate.C
\*---------------------------------------------------------------------------*/
#ifndef searchablePlate_H
#define searchablePlate_H
#include "searchableSurface.H"
#include "treeBoundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchablePlate Declaration
\*---------------------------------------------------------------------------*/
class searchablePlate
:
public searchableSurface
{
private:
// Private Member Data
const point origin_;
const vector span_;
//- Coordinate direction which is normal
const direction normalDir_;
mutable wordList regions_;
// Private Member Functions
//- Calculate normal direction from span
static direction calcNormal(const point&);
pointIndexHit findNearest
(
const point& sample,
const scalar nearestDistSqr
) const;
pointIndexHit findLine
(
const point& start,
const point& end
) const;
//- Disallow default bitwise copy construct
searchablePlate(const searchablePlate&);
//- Disallow default bitwise assignment
void operator=(const searchablePlate&);
public:
//- Runtime type information
TypeName("searchablePlate");
// Constructors
//- Construct from components
searchablePlate
(
const IOobject& io,
const point& origin,
const point& span
);
//- Construct from dictionary (used by searchableSurface)
searchablePlate
(
const IOobject& io,
const dictionary& dict
);
// Destructor
virtual ~searchablePlate();
// Member Functions
virtual const wordList& regions() const;
//- Whether supports volume type below
virtual bool hasVolumeType() const
{
return false;
}
//- Range of local indices that can be returned.
virtual label size() const
{
return 1;
}
// Multiple point queries.
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>&
) const;
virtual void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
virtual void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const;
//- From a set of points and indices get the region
virtual void getRegion
(
const List<pointIndexHit>&,
labelList& region
) const;
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>&,
vectorField& normal
) const;
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual void getVolumeType
(
const pointField&,
List<volumeType>&
) const;
// regIOobject implementation
bool writeData(Ostream&) const
{
notImplemented("searchablePlate::writeData(Ostream&) const");
return false;
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //

View File

@ -0,0 +1,321 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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 "searchableSurfaceWithGaps.H"
#include "addToRunTimeSelectionTable.H"
#include "SortableList.H"
#include "Time.H"
#include "ListOps.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(searchableSurfaceWithGaps, 0);
addToRunTimeSelectionTable(searchableSurface, searchableSurfaceWithGaps, dict);
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::Pair<Foam::vector> Foam::searchableSurfaceWithGaps::offsetVecs
(
const point& start,
const point& end
) const
{
Pair<vector> offsets(vector::zero, vector::zero);
vector n(end-start);
scalar magN = mag(n);
if (magN > SMALL)
{
n /= magN;
// Do first offset vector. Is the coordinate axes with the smallest
// component along the vector n.
scalar minMag = GREAT;
direction minCmpt = 0;
for (direction cmpt = 0; cmpt < vector::nComponents; cmpt++)
{
if (mag(n[cmpt]) < minMag)
{
minMag = mag(n[cmpt]);
minCmpt = cmpt;
}
}
offsets[0][minCmpt] = 1.0;
// Orthogonalise
offsets[0] -= n[minCmpt]*n;
// Scale
offsets[0] *= gap_/mag(offsets[0]);
// Do second offset vector perp to original edge and first offset vector
offsets[1] = n ^ offsets[0];
offsets[1] *= gap_/mag(offsets[1]);
}
return offsets;
}
void Foam::searchableSurfaceWithGaps::offsetVecs
(
const pointField& start,
const pointField& end,
pointField& offset0,
pointField& offset1
) const
{
offset0.setSize(start.size());
offset1.setSize(start.size());
forAll(start, i)
{
const Pair<vector> offsets(offsetVecs(start[i], end[i]));
offset0[i] = offsets[0];
offset1[i] = offsets[1];
}
}
Foam::label Foam::searchableSurfaceWithGaps::countMisses
(
const List<pointIndexHit>& info,
labelList& missMap
)
{
label nMiss = 0;
forAll(info, i)
{
if (!info[i].hit())
{
nMiss++;
}
}
missMap.setSize(nMiss);
nMiss = 0;
forAll(info, i)
{
if (!info[i].hit())
{
missMap[nMiss++] = i;
}
}
return nMiss;
}
// Anything not a hit in both counts as a hit
Foam::label Foam::searchableSurfaceWithGaps::countMisses
(
const List<pointIndexHit>& plusInfo,
const List<pointIndexHit>& minInfo,
labelList& missMap
)
{
label nMiss = 0;
forAll(plusInfo, i)
{
if (!plusInfo[i].hit() || !minInfo[i].hit())
{
nMiss++;
}
}
missMap.setSize(nMiss);
nMiss = 0;
forAll(plusInfo, i)
{
if (!plusInfo[i].hit() || !minInfo[i].hit())
{
missMap[nMiss++] = i;
}
}
return nMiss;
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableSurfaceWithGaps::searchableSurfaceWithGaps
(
const IOobject& io,
const dictionary& dict
)
:
searchableSurface(io),
gap_(readScalar(dict.lookup("gap"))),
subGeom_(1)
{
const word subGeomName(dict.lookup("surface"));
const searchableSurface& s =
io.db().lookupObject<searchableSurface>(subGeomName);
subGeom_.set(0, &const_cast<searchableSurface&>(s));
}
// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
Foam::searchableSurfaceWithGaps::~searchableSurfaceWithGaps()
{}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
void Foam::searchableSurfaceWithGaps::findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
surface().findLine(start, end, info);
// Count number of misses. Determine map
labelList compactMap;
label nMiss = countMisses(info, compactMap);
if (returnReduce(nMiss, sumOp<label>()) > 0)
{
pointField compactStart(start, compactMap);
pointField compactEnd(end, compactMap);
// Calculate offset vector
pointField offset0, offset1;
offsetVecs
(
compactStart,
compactEnd,
offset0,
offset1
);
// test in pairs: only if both perturbations hit something
// do we accept the hit.
List<pointIndexHit> plusInfo;
surface().findLine(start+offset0, end+offset0, plusInfo);
List<pointIndexHit> minInfo;
surface().findLine(start-offset0, end-offset0, minInfo);
// Extract any hits
forAll(plusInfo, i)
{
if (plusInfo[i].hit() && minInfo[i].hit())
{
info[compactMap[i]] = plusInfo[i].hitPoint()-offset0[i];
}
}
labelList plusMissMap;
nMiss = countMisses(plusInfo, minInfo, plusMissMap);
if (returnReduce(nMiss, sumOp<label>()) > 0)
{
// Extract (inplace possible because of order)
forAll(plusMissMap, i)
{
label mapI = plusMissMap[i];
compactStart[i] = compactStart[mapI];
compactEnd[i] = compactEnd[mapI];
compactMap[i] = compactMap[mapI];
}
compactStart.setSize(plusMissMap.size());
compactEnd.setSize(plusMissMap.size());
compactMap.setSize(plusMissMap.size());
surface().findLine(start+offset1, end+offset1, plusInfo);
surface().findLine(start-offset1, end-offset1, minInfo);
// Extract any hits
forAll(plusInfo, i)
{
if (plusInfo[i].hit() && minInfo[i].hit())
{
info[compactMap[i]] = plusInfo[i].hitPoint()-offset1[i];
}
}
}
}
}
void Foam::searchableSurfaceWithGaps::findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>& info
) const
{
// To be done ...
findLine(start, end, info);
}
void Foam::searchableSurfaceWithGaps::findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >& info
) const
{
// To be done. Assume for now only one intersection.
List<pointIndexHit> nearestInfo;
findLine(start, end, nearestInfo);
info.setSize(start.size());
forAll(info, pointI)
{
if (nearestInfo[pointI].hit())
{
info[pointI].setSize(1);
info[pointI][0] = nearestInfo[pointI];
}
else
{
info[pointI].clear();
}
}
}
// ************************************************************************* //

View File

@ -0,0 +1,238 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 1991-2008 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
Class
Foam::searchableSurfaceWithGaps
Description
searchableSurface using multiple slightly shifted underlying surfaces
to make sure pierces don't go through gaps.
SourceFiles
searchableSurfaceWithGaps.C
\*---------------------------------------------------------------------------*/
#ifndef searchableSurfaceWithGaps_H
#define searchableSurfaceWithGaps_H
#include "searchableSurface.H"
#include "UPtrList.H"
#include "Pair.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
// Forward declaration of classes
/*---------------------------------------------------------------------------*\
Class searchableSurfaceWithGaps Declaration
\*---------------------------------------------------------------------------*/
class searchableSurfaceWithGaps
:
public searchableSurface
{
private:
// Private Member Data
//- gap size in meter
const scalar gap_;
//- Underlying geometry (size 1)
UPtrList<searchableSurface> subGeom_;
// Private Member Functions
Pair<vector> offsetVecs(const point&, const point&) const;
void offsetVecs
(
const pointField& start,
const pointField& end,
pointField& offset0,
pointField& offset1
) const;
static label countMisses
(
const List<pointIndexHit>& info,
labelList& missMap
);
static label countMisses
(
const List<pointIndexHit>& plusInfo,
const List<pointIndexHit>& minInfo,
labelList& missMap
);
//- Disallow default bitwise copy construct
searchableSurfaceWithGaps(const searchableSurfaceWithGaps&);
//- Disallow default bitwise assignment
void operator=(const searchableSurfaceWithGaps&);
public:
//- Runtime type information
TypeName("searchableSurfaceWithGaps");
// Constructors
//- Construct from dictionary (used by searchableSurface)
searchableSurfaceWithGaps
(
const IOobject& io,
const dictionary& dict
);
// Destructor
virtual ~searchableSurfaceWithGaps();
// Member Functions
const searchableSurface& surface() const
{
return subGeom_[0];
}
virtual const wordList& regions() const
{
return surface().regions();
}
//- Whether supports volume type below
virtual bool hasVolumeType() const
{
return surface().hasVolumeType();
}
//- Range of local indices that can be returned.
virtual label size() const
{
return surface().size();
}
// Multiple point queries.
virtual void findNearest
(
const pointField& sample,
const scalarField& nearestDistSqr,
List<pointIndexHit>& info
) const
{
surface().findNearest
(
sample,
nearestDistSqr,
info
);
}
virtual void findLine
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
virtual void findLineAny
(
const pointField& start,
const pointField& end,
List<pointIndexHit>&
) const;
//- Get all intersections in order from start to end.
virtual void findLineAll
(
const pointField& start,
const pointField& end,
List<List<pointIndexHit> >&
) const;
//- From a set of points and indices get the region
virtual void getRegion
(
const List<pointIndexHit>& info,
labelList& region
) const
{
surface().getRegion(info, region);
}
//- From a set of points and indices get the normal
virtual void getNormal
(
const List<pointIndexHit>& info,
vectorField& normal
) const
{
surface().getNormal(info, normal);
}
//- Determine type (inside/outside/mixed) for point. unknown if
// cannot be determined (e.g. non-manifold surface)
virtual void getVolumeType
(
const pointField& samples,
List<volumeType>& info
) const
{
surface().getVolumeType(samples, info);
}
// regIOobject implementation
bool writeData(Ostream& os) const
{
return surface().writeData(os);
}
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //