Commit 8f656fff authored by mattijs's avatar mattijs
Browse files

Added functions to searchableSurface to get coordinates and store/retrieve an elementwise field.

Added numbering on searhableSurfaceCollection.
parent 4d4276fd
......@@ -912,41 +912,6 @@ Foam::distributedTriSurfaceMesh::independentlyDistributedBbs
}
void Foam::distributedTriSurfaceMesh::calcBounds
(
boundBox& bb,
label& nPoints
) const
{
// Unfortunately nPoints constructs meshPoints() so do compact version
// ourselves
PackedBoolList pointIsUsed(points().size());
nPoints = 0;
bb.min() = point(VGREAT, VGREAT, VGREAT);
bb.max() = point(-VGREAT, -VGREAT, -VGREAT);
const triSurface& s = static_cast<const triSurface&>(*this);
forAll(s, triI)
{
const labelledTri& f = s[triI];
forAll(f, fp)
{
label pointI = f[fp];
if (pointIsUsed.set(pointI, 1u))
{
bb.min() = ::Foam::min(bb.min(), points()[pointI]);
bb.max() = ::Foam::max(bb.max(), points()[pointI]);
nPoints++;
}
}
}
}
// Does any part of triangle overlap bb.
bool Foam::distributedTriSurfaceMesh::overlaps
(
......@@ -1935,20 +1900,7 @@ void Foam::distributedTriSurfaceMesh::getNormal
{
if (!Pstream::parRun())
{
normal.setSize(info.size());
forAll(info, i)
{
if (info[i].hit())
{
normal[i] = faceNormals()[info[i].index()];
}
else
{
// Set to what?
normal[i] = vector::zero;
}
}
triSurfaceMesh::getNormal(info, normal);
return;
}
......@@ -2006,70 +1958,64 @@ void Foam::distributedTriSurfaceMesh::getNormal
void Foam::distributedTriSurfaceMesh::getField
(
const word& fieldName,
const List<pointIndexHit>& info,
labelList& values
) const
{
const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
(
fieldName
);
if (!Pstream::parRun())
{
values.setSize(info.size());
forAll(info, i)
{
if (info[i].hit())
{
values[i] = fld[info[i].index()];
}
}
triSurfaceMesh::getField(info, values);
return;
}
if (foundObject<triSurfaceLabelField>("values"))
{
const triSurfaceLabelField& fld = lookupObject<triSurfaceLabelField>
(
"values"
);
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
labelList triangleIndex(info.size());
autoPtr<mapDistribute> mapPtr
(
calcLocalQueries
// Get query data (= local index of triangle)
// ~~~~~~~~~~~~~~
labelList triangleIndex(info.size());
autoPtr<mapDistribute> mapPtr
(
info,
triangleIndex
)
);
const mapDistribute& map = mapPtr();
calcLocalQueries
(
info,
triangleIndex
)
);
const mapDistribute& map = mapPtr();
// Do my tests
// ~~~~~~~~~~~
// Do my tests
// ~~~~~~~~~~~
values.setSize(triangleIndex.size());
values.setSize(triangleIndex.size());
forAll(triangleIndex, i)
{
label triI = triangleIndex[i];
values[i] = fld[triI];
}
forAll(triangleIndex, i)
{
label triI = triangleIndex[i];
values[i] = fld[triI];
}
// Send back results
// ~~~~~~~~~~~~~~~~~
// Send back results
// ~~~~~~~~~~~~~~~~~
map.distribute
(
Pstream::nonBlocking,
List<labelPair>(0),
info.size(),
map.constructMap(), // what to send
map.subMap(), // what to receive
values
);
map.distribute
(
Pstream::nonBlocking,
List<labelPair>(0),
info.size(),
map.constructMap(), // what to send
map.subMap(), // what to receive
values
);
}
}
......
......@@ -217,9 +217,6 @@ private:
const triSurface&
);
//- Calculate surface bounding box
void calcBounds(boundBox& bb, label& nPoints) const;
//- Does any part of triangle overlap bb.
static bool overlaps
(
......@@ -418,7 +415,7 @@ public:
// Should really be split into a routine to determine decomposition
// and one that does actual distribution but determining
// decomposition with duplicate triangle merging requires
// same amoun as work as actual distribution.
// same amount as work as actual distribution.
virtual void distribute
(
const List<treeBoundBox>&,
......@@ -430,14 +427,9 @@ public:
// Other
//- Specific to triSurfaceMesh: from a set of hits (points and
//- WIP. From a set of hits (points and
// indices) get the specified field. Misses do not get set.
virtual void getField
(
const word& fieldName,
const List<pointIndexHit>&,
labelList& values
) const;
virtual void getField(const List<pointIndexHit>&, labelList&) const;
//- Subset the part of surface that is overlapping bounds.
static triSurface overlappingSurface
......
......@@ -229,6 +229,21 @@ const Foam::wordList& Foam::searchableBox::regions() const
}
Foam::pointField Foam::searchableBox::coordinates() const
{
pointField ctrs(6);
const pointField pts = treeBoundBox::points();
const faceList& fcs = treeBoundBox::faces;
forAll(fcs, i)
{
ctrs[i] = fcs[i].centre(pts);
}
return ctrs;
}
Foam::pointIndexHit Foam::searchableBox::findNearest
(
const point& sample,
......
......@@ -127,6 +127,9 @@ public:
return 6;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const;
// Single point queries.
......
......@@ -40,6 +40,14 @@ addToRunTimeSelectionTable(searchableSurface, searchableCylinder, dict);
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
Foam::pointField Foam::searchableCylinder::coordinates() const
{
pointField ctrs(1, 0.5*(point1_ + point2_));
return ctrs;
}
Foam::pointIndexHit Foam::searchableCylinder::findNearest
(
const point& sample,
......
......@@ -150,6 +150,10 @@ public:
return 1;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const;
// Multiple point queries.
......
......@@ -26,7 +26,7 @@ Class
Foam::searchablePlane
Description
Searching on plane. See plane.H
Searching on (infinite) plane. See plane.H
SourceFiles
searchablePlane.C
......@@ -122,6 +122,14 @@ public:
return 1;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const
{
//notImplemented("searchablePlane::coordinates()")
return pointField(1, refPoint());
}
// Multiple point queries.
......
......@@ -145,6 +145,13 @@ public:
return 1;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const
{
return pointField(1, origin_);
}
// Multiple point queries.
......
......@@ -133,6 +133,13 @@ public:
return 1;
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const
{
return pointField(1, centre_);
}
// Multiple point queries.
......
......@@ -190,6 +190,10 @@ public:
return size();
}
//- Get representative set of element coordinates
// Usually the element centres (should be of length size()).
virtual pointField coordinates() const = 0;
// Single point queries.
......@@ -319,6 +323,18 @@ public:
)
{}
//- WIP. Store element-wise field.
virtual void setField(const labelList& values)
{}
//- WIP. From a set of hits (points and
// indices) get the specified field. Misses do not get set. Return
// empty field if not supported.
virtual void getField(const List<pointIndexHit>&, labelList& values)
const
{
values.clear();
}
};
......
......@@ -101,7 +101,11 @@ void Foam::searchableSurfaceCollection::findNearest
minDistSqr[pointI] = distSqr;
nearestInfo[pointI].setPoint(globalPt);
nearestInfo[pointI].setHit();
nearestInfo[pointI].setIndex(hitInfo[pointI].index());
nearestInfo[pointI].setIndex
(
hitInfo[pointI].index()
+ indexOffset_[surfI]
);
nearestSurf[pointI] = surfI;
}
}
......@@ -110,6 +114,62 @@ void Foam::searchableSurfaceCollection::findNearest
}
// Sort hits into per-surface bins. Misses are rejected. Maintains map back
// to position
void Foam::searchableSurfaceCollection::sortHits
(
const List<pointIndexHit>& info,
List<List<pointIndexHit> >& surfInfo,
labelListList& infoMap
) const
{
// Count hits per surface.
labelList nHits(subGeom_.size(), 0);
forAll(info, pointI)
{
if (info[pointI].hit())
{
label index = info[pointI].index();
label surfI = findLower(indexOffset_, index+1);
nHits[surfI]++;
}
}
// Per surface the hit
surfInfo.setSize(subGeom_.size());
// Per surface the original position
infoMap.setSize(subGeom_.size());
forAll(surfInfo, surfI)
{
surfInfo[surfI].setSize(nHits[surfI]);
infoMap[surfI].setSize(nHits[surfI]);
}
nHits = 0;
forAll(info, pointI)
{
if (info[pointI].hit())
{
label index = info[pointI].index();
label surfI = findLower(indexOffset_, index+1);
// Store for correct surface and adapt indices back to local
// ones
label localI = nHits[surfI]++;
surfInfo[surfI][localI] = pointIndexHit
(
info[pointI].hit(),
info[pointI].rawPoint(),
index-indexOffset_[surfI]
);
infoMap[surfI][localI] = pointI;
}
}
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
Foam::searchableSurfaceCollection::searchableSurfaceCollection
......@@ -123,11 +183,13 @@ Foam::searchableSurfaceCollection::searchableSurfaceCollection
scale_(dict.size()),
transform_(dict.size()),
subGeom_(dict.size()),
mergeSubRegions_(dict.lookup("mergeSubRegions"))
mergeSubRegions_(dict.lookup("mergeSubRegions")),
indexOffset_(dict.size()+1)
{
Info<< "SearchableCollection : " << name() << endl;
label surfI = 0;
label startIndex = 0;
forAllConstIter(dictionary, dict, iter)
{
if (dict.isDict(iter().keyword()))
......@@ -153,8 +215,24 @@ Foam::searchableSurfaceCollection::searchableSurfaceCollection
const searchableSurface& s =
io.db().lookupObject<searchableSurface>(subGeomName);
// I don't know yet how to handle the globalSize combined with
// regionOffset. Would cause non-consecutive indices locally
// if all indices offset by globalSize() of the local region...
if (s.size() != s.globalSize())
{
FatalErrorIn
(
"searchableSurfaceCollection::searchableSurfaceCollection"
"(const IOobject&, const dictionary&)"
) << "Cannot use a distributed surface in a collection."
<< exit(FatalError);
}
subGeom_.set(surfI, &const_cast<searchableSurface&>(s));
indexOffset_[surfI] = startIndex;
startIndex += subGeom_[surfI].size();
Info<< " instance : " << instance_[surfI] << endl;
Info<< " surface : " << s.name() << endl;
Info<< " scale : " << scale_[surfI] << endl;
......@@ -163,10 +241,13 @@ Foam::searchableSurfaceCollection::searchableSurfaceCollection
surfI++;
}
}
indexOffset_[surfI] = startIndex;
instance_.setSize(surfI);
scale_.setSize(surfI);
transform_.setSize(surfI);
subGeom_.setSize(surfI);
indexOffset_.setSize(surfI+1);
}
......@@ -212,12 +293,36 @@ const Foam::wordList& Foam::searchableSurfaceCollection::regions() const
Foam::label Foam::searchableSurfaceCollection::size() const
{
label n = 0;
return indexOffset_[indexOffset_.size()-1];
}
Foam::pointField Foam::searchableSurfaceCollection::coordinates() const
{
// Get overall size
pointField coords(size());
// Append individual coordinates
label coordI = 0;
forAll(subGeom_, surfI)
{
n += subGeom_[surfI].size();
const pointField subCoords = subGeom_[surfI].coordinates();
forAll(subCoords, i)
{
coords[coordI++] = transform_[surfI].globalPosition
(
cmptMultiply
(
subCoords[i],
scale_[surfI]
)
);
}
}
return n;
return coords;
}
......@@ -296,6 +401,11 @@ void Foam::searchableSurfaceCollection::findLine
);
info[pointI] = hitInfo[pointI];
info[pointI].rawPoint() = nearest[pointI];
info[pointI].setIndex
(
hitInfo[pointI].index()
+ indexOffset_[surfI]
);
}
}
}
......@@ -397,82 +507,42 @@ void Foam::searchableSurfaceCollection::getRegion
}
else
{
region.setSize(info.size());
region = -1;
// Multiple surfaces. Sort by surface.
// Which region did point come from. Retest for now to see which
// surface it originates from - crap solution! Should use global indices
// in index inside pointIndexHit to do this better.
// Per surface the hit
List<List<pointIndexHit> > surfInfo;
// Per surface the original position
List<List<label> > infoMap;
sortHits(info, surfInfo, infoMap);
pointField samples(info.size());
forAll(info, pointI)
{
if (info[pointI].hit())
{
samples[pointI] = info[pointI].hitPoint();
}
else
{
samples[pointI] = vector::zero;
}
}
//scalarField minDistSqr(info.size(), SMALL);
scalarField minDistSqr(info.size(), GREAT);
region.setSize(info.size());
region = -1;
labelList nearestSurf;
List<pointIndexHit> nearestInfo;
findNearest
(
samples,
minDistSqr,
nearestInfo,
nearestSurf
);
// Do region tests
// Check
if (mergeSubRegions_)
{
forAll(info, pointI)
// Actually no need for surfInfo. Just take region for surface.
forAll(infoMap, surfI)
{
if (info[pointI].hit() && nearestSurf[pointI] == -1)
const labelList& map = infoMap[surfI];
forAll(map, i)
{
FatalErrorIn
(
"searchableSurfaceCollection::getRegion(..)"
) << "pointI:" << pointI
<< " sample:" << samples[pointI]
<< " nearest:" << nearestInfo[pointI]
<< " nearestsurf:" << nearestSurf[pointI]
<< abort(FatalError);
region[map[i]] = regionOffset_[surfI];
}
}