From 8f656fffab81e77f4c8a58362c9ccab60fd8c5c7 Mon Sep 17 00:00:00 2001 From: mattijs Date: Fri, 15 Jan 2010 17:34:26 +0000 Subject: [PATCH] Added functions to searchableSurface to get coordinates and store/retrieve an elementwise field. Added numbering on searhableSurfaceCollection. --- .../distributedTriSurfaceMesh.C | 150 +++----- .../distributedTriSurfaceMesh.H | 14 +- .../searchableSurface/searchableBox.C | 15 + .../searchableSurface/searchableBox.H | 3 + .../searchableSurface/searchableCylinder.C | 8 + .../searchableSurface/searchableCylinder.H | 4 + .../searchableSurface/searchablePlane.H | 10 +- .../searchableSurface/searchablePlate.H | 7 + .../searchableSurface/searchableSphere.H | 7 + .../searchableSurface/searchableSurface.H | 16 + .../searchableSurfaceCollection.C | 355 ++++++++++++------ .../searchableSurfaceCollection.H | 38 ++ .../searchableSurfaceWithGaps.H | 42 +++ .../searchableSurface/triSurfaceMesh.C | 147 ++++++-- .../searchableSurface/triSurfaceMesh.H | 26 +- 15 files changed, 587 insertions(+), 255 deletions(-) diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C index 60ad087b70..40ad36bbc3 100644 --- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C +++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.C @@ -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(*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& info, labelList& values ) const { - const triSurfaceLabelField& fld = lookupObject - ( - 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; } - - // Get query data (= local index of triangle) - // ~~~~~~~~~~~~~~ - - labelList triangleIndex(info.size()); - autoPtr mapPtr - ( - calcLocalQueries - ( - info, - triangleIndex - ) - ); - const mapDistribute& map = mapPtr(); - - - // Do my tests - // ~~~~~~~~~~~ - - values.setSize(triangleIndex.size()); - - forAll(triangleIndex, i) + if (foundObject("values")) { - label triI = triangleIndex[i]; - values[i] = fld[triI]; + const triSurfaceLabelField& fld = lookupObject + ( + "values" + ); + + + // Get query data (= local index of triangle) + // ~~~~~~~~~~~~~~ + + labelList triangleIndex(info.size()); + autoPtr mapPtr + ( + calcLocalQueries + ( + info, + triangleIndex + ) + ); + const mapDistribute& map = mapPtr(); + + + // Do my tests + // ~~~~~~~~~~~ + + values.setSize(triangleIndex.size()); + + forAll(triangleIndex, i) + { + label triI = triangleIndex[i]; + values[i] = fld[triI]; + } + + + // Send back results + // ~~~~~~~~~~~~~~~~~ + + map.distribute + ( + Pstream::nonBlocking, + List(0), + info.size(), + map.constructMap(), // what to send + map.subMap(), // what to receive + values + ); } - - - // Send back results - // ~~~~~~~~~~~~~~~~~ - - map.distribute - ( - Pstream::nonBlocking, - List(0), - info.size(), - map.constructMap(), // what to send - map.subMap(), // what to receive - values - ); } diff --git a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H index dba72d2abf..577d12e4c3 100644 --- a/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H +++ b/src/meshTools/searchableSurface/distributedTriSurfaceMesh.H @@ -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&, @@ -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&, - labelList& values - ) const; + virtual void getField(const List&, labelList&) const; //- Subset the part of surface that is overlapping bounds. static triSurface overlappingSurface diff --git a/src/meshTools/searchableSurface/searchableBox.C b/src/meshTools/searchableSurface/searchableBox.C index a4cc593562..cb854c3772 100644 --- a/src/meshTools/searchableSurface/searchableBox.C +++ b/src/meshTools/searchableSurface/searchableBox.C @@ -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, diff --git a/src/meshTools/searchableSurface/searchableBox.H b/src/meshTools/searchableSurface/searchableBox.H index 9d1f2e8217..77e0ecd98d 100644 --- a/src/meshTools/searchableSurface/searchableBox.H +++ b/src/meshTools/searchableSurface/searchableBox.H @@ -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. diff --git a/src/meshTools/searchableSurface/searchableCylinder.C b/src/meshTools/searchableSurface/searchableCylinder.C index fcaa42194a..0df222ccb0 100644 --- a/src/meshTools/searchableSurface/searchableCylinder.C +++ b/src/meshTools/searchableSurface/searchableCylinder.C @@ -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, diff --git a/src/meshTools/searchableSurface/searchableCylinder.H b/src/meshTools/searchableSurface/searchableCylinder.H index cae0f058db..98b4d91e41 100644 --- a/src/meshTools/searchableSurface/searchableCylinder.H +++ b/src/meshTools/searchableSurface/searchableCylinder.H @@ -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. diff --git a/src/meshTools/searchableSurface/searchablePlane.H b/src/meshTools/searchableSurface/searchablePlane.H index 87777e3209..c0b58cc8ac 100644 --- a/src/meshTools/searchableSurface/searchablePlane.H +++ b/src/meshTools/searchableSurface/searchablePlane.H @@ -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. diff --git a/src/meshTools/searchableSurface/searchablePlate.H b/src/meshTools/searchableSurface/searchablePlate.H index 7bc900f554..b05414bb6c 100644 --- a/src/meshTools/searchableSurface/searchablePlate.H +++ b/src/meshTools/searchableSurface/searchablePlate.H @@ -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. diff --git a/src/meshTools/searchableSurface/searchableSphere.H b/src/meshTools/searchableSurface/searchableSphere.H index 297890b9db..fb4aeb22f8 100644 --- a/src/meshTools/searchableSurface/searchableSphere.H +++ b/src/meshTools/searchableSurface/searchableSphere.H @@ -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. diff --git a/src/meshTools/searchableSurface/searchableSurface.H b/src/meshTools/searchableSurface/searchableSurface.H index 35b9dc9fb5..daf8b76daa 100644 --- a/src/meshTools/searchableSurface/searchableSurface.H +++ b/src/meshTools/searchableSurface/searchableSurface.H @@ -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&, labelList& values) + const + { + values.clear(); + } }; diff --git a/src/meshTools/searchableSurface/searchableSurfaceCollection.C b/src/meshTools/searchableSurface/searchableSurfaceCollection.C index ba1dd32988..9279087c19 100644 --- a/src/meshTools/searchableSurface/searchableSurfaceCollection.C +++ b/src/meshTools/searchableSurface/searchableSurfaceCollection.C @@ -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& info, + List >& 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(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(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 { + // Multiple surfaces. Sort by surface. + + // Per surface the hit + List > surfInfo; + // Per surface the original position + List > infoMap; + sortHits(info, surfInfo, infoMap); + region.setSize(info.size()); region = -1; - // 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. + // Do region tests - pointField samples(info.size()); - forAll(info, pointI) + if (mergeSubRegions_) { - if (info[pointI].hit()) + // Actually no need for surfInfo. Just take region for surface. + forAll(infoMap, surfI) { - samples[pointI] = info[pointI].hitPoint(); - } - else - { - samples[pointI] = vector::zero; - } - } - //scalarField minDistSqr(info.size(), SMALL); - scalarField minDistSqr(info.size(), GREAT); - - labelList nearestSurf; - List nearestInfo; - findNearest - ( - samples, - minDistSqr, - nearestInfo, - nearestSurf - ); - - // Check - { - forAll(info, pointI) - { - 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]; } } } - - forAll(subGeom_, surfI) + else { - // Collect points from my surface - labelList indices(findIndices(nearestSurf, surfI)); - - if (mergeSubRegions_) - { - forAll(indices, i) - { - region[indices[i]] = regionOffset_[surfI]; - } - } - else + forAll(infoMap, surfI) { labelList surfRegion; - subGeom_[surfI].getRegion - ( - List - ( - UIndirectList(info, indices) - ), - surfRegion - ); - forAll(indices, i) + subGeom_[surfI].getRegion(surfInfo[surfI], surfRegion); + + const labelList& map = infoMap[surfI]; + forAll(map, i) { - region[indices[i]] = regionOffset_[surfI] + surfRegion[i]; + region[map[i]] = regionOffset_[surfI] + surfRegion[i]; } } } @@ -494,52 +564,26 @@ void Foam::searchableSurfaceCollection::getNormal } else { + // Multiple surfaces. Sort by surface. + + // Per surface the hit + List > surfInfo; + // Per surface the original position + List > infoMap; + sortHits(info, surfInfo, infoMap); + normal.setSize(info.size()); - // See above - crap retest to find surface point originates from. - pointField samples(info.size()); - forAll(info, pointI) + // Do region tests + forAll(surfInfo, surfI) { - if (info[pointI].hit()) - { - samples[pointI] = info[pointI].hitPoint(); - } - else - { - samples[pointI] = vector::zero; - } - } - //scalarField minDistSqr(info.size(), SMALL); - scalarField minDistSqr(info.size(), GREAT); - - labelList nearestSurf; - List nearestInfo; - findNearest - ( - samples, - minDistSqr, - nearestInfo, - nearestSurf - ); - - - forAll(subGeom_, surfI) - { - // Collect points from my surface - labelList indices(findIndices(nearestSurf, surfI)); - vectorField surfNormal; - subGeom_[surfI].getNormal - ( - List - ( - UIndirectList(info, indices) - ), - surfNormal - ); - forAll(indices, i) + subGeom_[surfI].getNormal(surfInfo[surfI], surfNormal); + + const labelList& map = infoMap[surfI]; + forAll(map, i) { - normal[indices[i]] = surfNormal[i]; + normal[map[i]] = surfNormal[i]; } } } @@ -561,4 +605,99 @@ void Foam::searchableSurfaceCollection::getVolumeType } +void Foam::searchableSurfaceCollection::distribute +( + const List& bbs, + const bool keepNonLocal, + autoPtr& faceMap, + autoPtr& pointMap +) +{ + forAll(subGeom_, surfI) + { + // Note:Tranform the bounding boxes? Something like + // pointField bbPoints = + // cmptDivide + // ( + // transform_[surfI].localPosition + // ( + // bbs[i].points() + // ), + // scale_[surfI] + // ); + // treeBoundBox newBb(bbPoints); + + // Note: what to do with faceMap, pointMap from multiple surfaces? + subGeom_[surfI].distribute + ( + bbs, + keepNonLocal, + faceMap, + pointMap + ); + } +} + + +void Foam::searchableSurfaceCollection::setField(const labelList& values) +{ + forAll(subGeom_, surfI) + { + subGeom_[surfI].setField + ( + static_cast + ( + SubList