ENH: change UnsortedMeshedSurface -> meshedSurface for sampledTriSurfaceMesh

- all sampled surface types now consistently use the same storage,
  which allows some more simplifications in the future.

- before/after comparison of the sampledTriSurfaceMesh tested with
  motorbike passenger helmet (serial and parallel). Use the newly added
  'keepIds' functionality to retain the original ids, and can also
  compare them to the original obj file with "GenerateIds" in paraview.
This commit is contained in:
Mark Olesen
2016-11-29 22:56:08 +01:00
parent 3c41b80b38
commit 8b75035f29
3 changed files with 224 additions and 101 deletions

View File

@ -78,6 +78,32 @@ namespace Foam
}
// * * * * * * * * * * * * * Static Member Functions * * * * * * * * * * * * //
void Foam::sampledTriSurfaceMesh::setZoneMap
(
const surfZoneList& zoneLst,
labelList& zoneIds
)
{
label sz = 0;
forAll(zoneLst, zonei)
{
const surfZone& zn = zoneLst[zonei];
sz += zn.size();
}
zoneIds.setSize(sz);
forAll(zoneLst, zonei)
{
const surfZone& zn = zoneLst[zonei];
// Assign sub-zone Ids
SubList<label>(zoneIds, zn.size(), zn.start()) = zonei;
}
}
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
const Foam::indexedOctree<Foam::treeDataFace>&
@ -146,12 +172,7 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
// Global numbering for cells/faces - only used to uniquely identify local
// elements
globalIndex globalCells
(
(sampleSource_ == cells || sampleSource_ == insideCells)
? mesh().nCells()
: mesh().nFaces()
);
globalIndex globalCells(onBoundary() ? mesh().nFaces() : mesh().nCells());
forAll(nearest, i)
{
@ -270,6 +291,37 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
// From original point to compact points
labelList reversePointMap(s.points().size(), -1);
// Handle region-wise sorting (makes things slightly more complicated)
zoneIds_.setSize(s.size(), -1);
// Better not to use triSurface::sortedZones here,
// since we'll sort ourselves
// Get zone/region sizes used, store under the original region Id
Map<label> zoneSizes;
// Recover region names from the input surface
Map<word> zoneNames;
{
const geometricSurfacePatchList& patches = s.patches();
forAll(patches, patchi)
{
zoneNames.set
(
patchi,
(
patches[patchi].name().empty()
? Foam::name("patch%d", patchi)
: patches[patchi].name()
)
);
zoneSizes.set(patchi, 0);
}
}
{
label newPointi = 0;
label newFacei = 0;
@ -278,57 +330,139 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
{
if (cellOrFaceLabels[facei] != -1)
{
faceMap[newFacei++] = facei;
const triSurface::FaceType& f = s[facei];
const label regionid = f.region();
Map<label>::iterator fnd = zoneSizes.find(regionid);
if (fnd != zoneSizes.end())
{
fnd()++;
}
else
{
// This shouldn't happen
zoneSizes.insert(regionid, 1);
zoneNames.set
(
regionid,
Foam::name("patch%d", regionid)
);
}
// Store new faces compact
faceMap[newFacei] = facei;
zoneIds_[newFacei] = regionid;
++newFacei;
// Renumber face points
forAll(f, fp)
{
if (reversePointMap[f[fp]] == -1)
const label labI = f[fp];
if (reversePointMap[labI] == -1)
{
pointMap[newPointi] = f[fp];
reversePointMap[f[fp]] = newPointi++;
pointMap[newPointi] = labI;
reversePointMap[labI] = newPointi++;
}
}
}
}
// Trim
faceMap.setSize(newFacei);
zoneIds_.setSize(newFacei);
pointMap.setSize(newPointi);
}
// Assign start/size (and name) to the newZones
// re-use the lookup to map (zoneId => zoneI)
surfZoneList zoneLst(zoneSizes.size());
label start = 0;
label zoneI = 0;
forAllIter(Map<label>, zoneSizes, iter)
{
// No negative regionids, so Map<label> sorts properly
const label regionid = iter.key();
word name;
Map<word>::const_iterator fnd = zoneNames.find(regionid);
if (fnd != zoneNames.end())
{
name = fnd();
}
if (name.empty())
{
name = ::Foam::name("patch%d", regionid);
}
zoneLst[zoneI] = surfZone
(
name,
0, // initialize with zero size
start,
zoneI
);
// Adjust start for the next zone and save (zoneId => zoneI) mapping
start += iter();
iter() = zoneI++;
}
// At this stage:
// - faceMap to map the (unsorted) compact to original triangle
// - zoneIds for the next sorting
// - zoneSizes contains region -> count information
// Rebuild the faceMap for the sorted order
labelList sortedFaceMap(faceMap.size());
forAll(zoneIds_, facei)
{
label zonei = zoneIds_[facei];
label sortedFacei = zoneLst[zonei].start() + zoneLst[zonei].size()++;
sortedFaceMap[sortedFacei] = faceMap[facei];
}
// zoneIds are now simply flat values
setZoneMap(zoneLst, zoneIds_);
// Replace the faceMap with the properly sorted face map
faceMap.transfer(sortedFaceMap);
if (keepIds_)
{
originalIds_ = faceMap;
}
// Subset cellOrFaceLabels
// Subset cellOrFaceLabels (for compact faces)
cellOrFaceLabels = UIndirectList<label>(cellOrFaceLabels, faceMap)();
// Store any face per point (without using pointFaces())
labelList pointToFace(pointMap.size());
// Create faces and points for subsetted surface
faceList& faces = this->storedFaces();
faces.setSize(faceMap.size());
faceList& surfFaces = this->storedFaces();
surfFaces.setSize(faceMap.size());
labelList& zoneIds = this->storedZoneIds();
zoneIds.setSize(faceMap.size());
this->storedZones().transfer(zoneLst);
forAll(faceMap, i)
forAll(faceMap, facei)
{
const labelledTri& f = s[faceMap[i]];
triFace newF
const labelledTri& origF = s[faceMap[facei]];
face& f = surfFaces[facei];
f = triFace
(
reversePointMap[f[0]],
reversePointMap[f[1]],
reversePointMap[f[2]]
reversePointMap[origF[0]],
reversePointMap[origF[1]],
reversePointMap[origF[2]]
);
faces[i] = newF.triFaceFace();
zoneIds[i] = f.region(); // preserve zone information
forAll(newF, fp)
forAll(f, fp)
{
pointToFace[newF[fp]] = i;
pointToFace[f[fp]] = facei;
}
}
@ -429,16 +563,16 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
else
{
// if sampleSource_ == cells:
// samplePoints_ : n/a
// sampleElements_ : per surface triangle the cell
// samplePoints_ : n/a
// if sampleSource_ == insideCells:
// samplePoints_ : n/a
// sampleElements_ : -1 or per surface triangle the cell
// else:
// samplePoints_ : n/a
// else:
// sampleElements_ : per surface triangle the boundary face
samplePoints_.clear();
// samplePoints_ : n/a
sampleElements_.transfer(cellOrFaceLabels);
samplePoints_.clear();
}
@ -448,12 +582,17 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
OFstream str(mesh().time().path()/"surfaceToMesh.obj");
Info<< "Dumping correspondence from local surface (points or faces)"
<< " to mesh (cells or faces) to " << str.name() << endl;
label vertI = 0;
const vectorField& centres =
(
onBoundary()
? mesh().faceCentres()
: mesh().cellCentres()
);
if (sampledSurface::interpolate())
{
if (sampleSource_ == cells || sampleSource_ == insideCells)
{
label vertI = 0;
forAll(samplePoints_, pointi)
{
meshTools::writeOBJ(str, points()[pointi]);
@ -462,62 +601,27 @@ bool Foam::sampledTriSurfaceMesh::update(const meshSearch& meshSearcher)
meshTools::writeOBJ(str, samplePoints_[pointi]);
vertI++;
label celli = sampleElements_[pointi];
meshTools::writeOBJ(str, mesh().cellCentres()[celli]);
label elemi = sampleElements_[pointi];
meshTools::writeOBJ(str, centres[elemi]);
vertI++;
str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI
<< nl;
str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI << nl;
}
}
else
{
forAll(samplePoints_, pointi)
{
meshTools::writeOBJ(str, points()[pointi]);
vertI++;
meshTools::writeOBJ(str, samplePoints_[pointi]);
vertI++;
label facei = sampleElements_[pointi];
meshTools::writeOBJ(str, mesh().faceCentres()[facei]);
vertI++;
str << "l " << vertI-2 << ' ' << vertI-1 << ' ' << vertI
<< nl;
}
}
}
else
{
if (sampleSource_ == cells || sampleSource_ == insideCells)
{
label vertI = 0;
forAll(sampleElements_, triI)
{
meshTools::writeOBJ(str, faceCentres()[triI]);
vertI++;
label celli = sampleElements_[triI];
meshTools::writeOBJ(str, mesh().cellCentres()[celli]);
vertI++;
str << "l " << vertI-1 << ' ' << vertI << nl;
}
}
else
{
forAll(sampleElements_, triI)
{
meshTools::writeOBJ(str, faceCentres()[triI]);
vertI++;
label facei = sampleElements_[triI];
meshTools::writeOBJ(str, mesh().faceCentres()[facei]);
label elemi = sampleElements_[triI];
meshTools::writeOBJ(str, centres[elemi]);
vertI++;
str << "l " << vertI-1 << ' ' << vertI << nl;
}
}
}
}
needsUpdate_ = false;
return true;
@ -553,6 +657,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
needsUpdate_(true),
keepIds_(false),
originalIds_(),
zoneIds_(),
sampleElements_(0),
samplePoints_(0)
{}
@ -584,6 +689,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
needsUpdate_(true),
keepIds_(dict.lookupOrDefault<Switch>("keepIds", false)),
originalIds_(),
zoneIds_(),
sampleElements_(0),
samplePoints_(0)
{}
@ -617,6 +723,7 @@ Foam::sampledTriSurfaceMesh::sampledTriSurfaceMesh
needsUpdate_(true),
keepIds_(false),
originalIds_(),
zoneIds_(),
sampleElements_(0),
samplePoints_(0)
{}
@ -646,6 +753,7 @@ bool Foam::sampledTriSurfaceMesh::expire()
sampledSurface::clearGeom();
MeshStorage::clear();
zoneIds_.clear();
originalIds_.clear();
boundaryTreePtr_.clear();

View File

@ -62,6 +62,7 @@ Description
SourceFiles
sampledTriSurfaceMesh.C
sampledTriSurfaceMeshTemplates.C
\*---------------------------------------------------------------------------*/
@ -71,6 +72,7 @@ SourceFiles
#include "sampledSurface.H"
#include "triSurfaceMesh.H"
#include "MeshedSurface.H"
#include "MeshedSurfacesFwd.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -87,7 +89,7 @@ class meshSearch;
class sampledTriSurfaceMesh
:
public sampledSurface,
public UnsortedMeshedSurface<face>
public meshedSurface
{
public:
//- Types of communications
@ -95,13 +97,13 @@ public:
{
cells,
insideCells,
boundaryFaces,
boundaryFaces
};
private:
//- Private typedefs for convenience
typedef UnsortedMeshedSurface<face> MeshStorage;
typedef meshedSurface MeshStorage;
// Private data
@ -127,6 +129,9 @@ private:
//- Search tree for all non-coupled boundary faces
mutable autoPtr<indexedOctree<treeDataFace>> boundaryTreePtr_;
//- For compatibility with the meshSurf interface
labelList zoneIds_;
//- From local surface triangle to mesh cell/face.
labelList sampleElements_;
@ -194,6 +199,10 @@ public:
// Member Functions
//- Set new zoneIds list based on the surfZoneList information
static void setZoneMap(const surfZoneList&, labelList& zoneIds);
//- Does the surface need an update?
virtual bool needsUpdate() const;
@ -226,7 +235,7 @@ public:
//- Const access to per-face zone/region information
virtual const labelList& zoneIds() const
{
return MeshStorage::zoneIds();
return zoneIds_;
}
//- Face area vectors
@ -260,6 +269,12 @@ public:
return originalIds_;
}
//- Sampling boundary values instead of cell values
bool onBoundary() const
{
return sampleSource_ == boundaryFaces;
}
//- Sample field on surface
virtual tmp<scalarField> sample
@ -323,7 +338,7 @@ public:
const interpolation<tensor>&
) const;
//- Write
//- Write information
virtual void print(Ostream&) const;
};

View File

@ -3,7 +3,7 @@
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation
\\/ M anipulation |
\\/ M anipulation | Copyright (C) 2016 OpenCFD Ltd.
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@ -38,7 +38,7 @@ Foam::sampledTriSurfaceMesh::sampleField
tmp<Field<Type>> tvalues(new Field<Type>(sampleElements_.size()));
Field<Type>& values = tvalues.ref();
if (sampleSource_ == cells || sampleSource_ == insideCells)
if (!onBoundary())
{
// Sample cells
@ -52,7 +52,7 @@ Foam::sampledTriSurfaceMesh::sampleField
// Sample boundary faces
const polyBoundaryMesh& pbm = mesh().boundaryMesh();
label nBnd = mesh().nFaces()-mesh().nInternalFaces();
const label nBnd = mesh().nFaces()-mesh().nInternalFaces();
// Create flat boundary field
@ -94,7 +94,7 @@ Foam::sampledTriSurfaceMesh::interpolateField
tmp<Field<Type>> tvalues(new Field<Type>(sampleElements_.size()));
Field<Type>& values = tvalues.ref();
if (sampleSource_ == cells || sampleSource_ == insideCells)
if (!onBoundary())
{
// Sample cells.