Commit 8b72dbed authored by Mark Olesen's avatar Mark Olesen
Browse files

sampling changes

- cuttingPlane and sampledPatch now use PrimitiveMeshedSurface for their
  storage.
- cuttingPlane always triangulates since the cut faces tend to look quite
  horrible anyhow.
- moved triangulation request out of sampledSurface to the inherited
  classes. It might disappear anyhow.
parent a1d7080d
EXE_INC = \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/surfMesh/lnInclude \
-I$(LIB_SRC)/triSurface/lnInclude \
-I$(LIB_SRC)/lagrangian/basic/lnInclude
LIB_LIBS = \
-lfiniteVolume \
-lmeshTools \
-lsurfMesh \
-ltriSurface
......@@ -45,7 +45,7 @@ void Foam::cuttingPlane::calcCutCells
(
const primitiveMesh& mesh,
const scalarField& dotProducts,
const labelList& cellIdLabels
const UList<label>& cellIdLabels
)
{
const labelListList& cellEdges = mesh.cellEdges();
......@@ -103,11 +103,12 @@ void Foam::cuttingPlane::calcCutCells
// Determine for each edge the intersection point. Calculates
// - cutPoints_ : coordinates of all intersection points
// - edgePoint : per edge -1 or the index into cutPoints_
Foam::labelList Foam::cuttingPlane::intersectEdges
// - edgePoint : per edge -1 or the index into cutPoints
void Foam::cuttingPlane::intersectEdges
(
const primitiveMesh& mesh,
const scalarField& dotProducts
const scalarField& dotProducts,
List<label>& edgePoint
)
{
// Use the dotProducts to find out the cut edges.
......@@ -115,7 +116,7 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
const pointField& points = mesh.points();
// Per edge -1 or the label of the intersection point
labelList edgePoint(edges.size(), -1);
edgePoint.setSize(edges.size());
DynamicList<point> dynCuttingPoints(4*cutCells_.size());
......@@ -129,7 +130,9 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
|| (dotProducts[e[1]] < zeroish && dotProducts[e[0]] > positive)
)
{
// Edge is cut.
// Edge is cut
edgePoint[edgeI] = dynCuttingPoints.size();
const point& p0 = points[e[0]];
const point& p1 = points[e[1]];
......@@ -139,7 +142,7 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
{
dynCuttingPoints.append(p0);
}
else if (alpha > 1.0)
else if (alpha >= 1.0)
{
dynCuttingPoints.append(p1);
}
......@@ -147,15 +150,14 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
{
dynCuttingPoints.append((1-alpha)*p0 + alpha*p1);
}
edgePoint[edgeI] = dynCuttingPoints.size() - 1;
}
else
{
edgePoint[edgeI] = -1;
}
}
dynCuttingPoints.shrink();
cutPoints_.transfer(dynCuttingPoints);
return edgePoint;
this->storedPoints().transfer(dynCuttingPoints);
}
......@@ -164,7 +166,7 @@ Foam::labelList Foam::cuttingPlane::intersectEdges
bool Foam::cuttingPlane::walkCell
(
const primitiveMesh& mesh,
const labelList& edgePoint,
const UList<label>& edgePoint,
const label cellI,
const label startEdgeI,
DynamicList<label>& faceVerts
......@@ -175,6 +177,7 @@ bool Foam::cuttingPlane::walkCell
label nIter = 0;
faceVerts.clear();
do
{
faceVerts.append(edgePoint[edgeI]);
......@@ -255,11 +258,17 @@ bool Foam::cuttingPlane::walkCell
void Foam::cuttingPlane::walkCellCuts
(
const primitiveMesh& mesh,
const labelList& edgePoint
const UList<label>& edgePoint
)
{
cutFaces_.setSize(cutCells_.size());
label cutFaceI = 0;
const pointField& cutPoints = this->points();
// use dynamic lists to handle triangulation and/or missed cuts
DynamicList<face> dynCutFaces(cutCells_.size());
DynamicList<label> dynCutCells(cutCells_.size());
// scratch space for calculating the face vertices
DynamicList<label> faceVerts(10);
forAll(cutCells_, i)
{
......@@ -290,7 +299,6 @@ void Foam::cuttingPlane::walkCellCuts
}
// Walk from starting edge around the circumference of the cell.
DynamicList<label> faceVerts(2*mesh.faces()[cellI].size());
bool okCut = walkCell
(
mesh,
......@@ -302,54 +310,46 @@ void Foam::cuttingPlane::walkCellCuts
if (okCut)
{
faceVerts.shrink();
face cutFace(faceVerts);
face f(faceVerts);
// Orient face.
if ((cutFace.normal(cutPoints_) && normal()) < 0)
// Orient face to point in the same direction as the plane normal
if ((f.normal(cutPoints) && normal()) < 0)
{
cutFace = cutFace.reverseFace();
f = f.reverseFace();
}
cutFaces_[cutFaceI++] = cutFace;
// the cut faces are usually quite ugly, so always triangulate
label nTri = f.triangles(cutPoints, dynCutFaces);
while (nTri--)
{
dynCutCells.append(cellI);
}
}
}
cutFaces_.setSize(cutFaceI);
this->storedFaces().transfer(dynCutFaces);
cutCells_.transfer(dynCutCells);
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
// Construct without cutting
Foam::cuttingPlane::cuttingPlane(const plane& newPlane)
Foam::cuttingPlane::cuttingPlane(const plane& pln)
:
plane(newPlane)
plane(pln)
{}
// Construct from components
Foam::cuttingPlane::cuttingPlane
(
const primitiveMesh& mesh,
const plane& newPlane
)
:
plane(newPlane)
{
reCut(mesh);
}
// Construct from mesh reference and plane, restricted to a list of cells
Foam::cuttingPlane::cuttingPlane
(
const primitiveMesh& mesh,
const plane& newPlane,
const labelList& cellIdLabels
const plane& pln,
const UList<label>& cellIdLabels
)
:
plane(newPlane)
plane(pln)
{
reCut(mesh, cellIdLabels);
}
......@@ -358,101 +358,51 @@ Foam::cuttingPlane::cuttingPlane
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
// reCut mesh with existing planeDesc
void Foam::cuttingPlane::reCut
(
const primitiveMesh& mesh
)
{
cutCells_.clear();
cutPoints_.clear();
cutFaces_.clear();
scalarField dotProducts = (mesh.points() - refPoint()) & normal();
//// Perturb points cuts so edges are cut properly.
//const tmp<scalarField> tdotProducts = stabilise(rawDotProducts, SMALL);
//const scalarField& dotProducts = tdotProducts();
// Determine cells that are (probably) cut.
calcCutCells(mesh, dotProducts);
// Determine cutPoints and return list of edge cuts. (per edge -1 or
// the label of the intersection point (in cutPoints_)
labelList edgePoint(intersectEdges(mesh, dotProducts));
// Do topological walk around cell to find closed loop.
walkCellCuts(mesh, edgePoint);
}
// recut mesh with existing planeDesc
void Foam::cuttingPlane::reCut
(
const primitiveMesh& mesh,
const labelList& cellIdLabels
const UList<label>& cellIdLabels
)
{
MeshStorage::clear();
cutCells_.clear();
cutPoints_.clear();
cutFaces_.clear();
scalarField dotProducts = (mesh.points() - refPoint()) & normal();
// Determine cells that are (probably) cut.
calcCutCells(mesh, dotProducts, cellIdLabels);
// Determine cutPoints and return list of edge cuts. (per edge -1 or
// the label of the intersection point (in cutPoints_)
labelList edgePoint(intersectEdges(mesh, dotProducts));
// Determine cutPoints and return list of edge cuts.
// per edge -1 or the label of the intersection point
labelList edgePoint;
intersectEdges(mesh, dotProducts, edgePoint);
// Do topological walk around cell to find closed loop.
walkCellCuts(mesh, edgePoint);
}
// Return plane used
const Foam::plane& Foam::cuttingPlane::planeDesc() const
{
return static_cast<const plane&>(*this);
}
// Return vectorField of cutting points
const Foam::pointField& Foam::cuttingPlane::points() const
{
return cutPoints_;
}
// Return unallocFaceList of points in cells
const Foam::faceList& Foam::cuttingPlane::faces() const
{
return cutFaces_;
}
// Return labelList of cut cells
const Foam::labelList& Foam::cuttingPlane::cells() const
{
return cutCells_;
}
bool Foam::cuttingPlane::cut()
// remap action on triangulation
void Foam::cuttingPlane::remapFaces
(
const UList<label>& faceMap
)
{
if (cutCells_.size() > 0)
{
return true;
}
else
// recalculate the cells cut
if (&faceMap && faceMap.size())
{
return false;
MeshStorage::remapFaces(faceMap);
List<label> newCutCells(faceMap.size());
forAll(faceMap, faceI)
{
newCutCells[faceI] = cutCells_[faceMap[faceI]];
}
cutCells_.transfer(newCutCells);
}
}
// * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * //
void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
......@@ -465,11 +415,9 @@ void Foam::cuttingPlane::operator=(const cuttingPlane& rhs)
<< abort(FatalError);
}
static_cast<MeshStorage&>(*this) = rhs;
static_cast<plane&>(*this) = rhs;
cutCells_ = rhs.cells();
cutPoints_ = rhs.points();
cutFaces_ = rhs.faces();
cutCells_ = rhs.cutCells();
}
......
......@@ -28,7 +28,8 @@ Class
Description
Constructs plane through mesh.
No attempt at resolving degenerate cases.
No attempt at resolving degenerate cases. Since the cut faces are
usually quite ugly, they will always be triangulated.
Note
When the cutting plane coincides with a mesh face, the cell edge on the
......@@ -45,6 +46,7 @@ SourceFiles
#include "plane.H"
#include "pointField.H"
#include "faceList.H"
#include "MeshedSurface.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......@@ -54,25 +56,22 @@ namespace Foam
class primitiveMesh;
/*---------------------------------------------------------------------------*\
Class cuttingPlane Declaration
Class cuttingPlane Declaration
\*---------------------------------------------------------------------------*/
class cuttingPlane
:
public PrimitiveMeshedSurface<face>,
public plane
{
//- Private typedefs for convenience
typedef PrimitiveMeshedSurface<face> MeshStorage;
// Private data
//- List of cells cut by the plane
labelList cutCells_;
//- Intersection points
pointField cutPoints_;
//- Cut faces in terms of intersection points
faceList cutFaces_;
// Private Member Functions
//- Determine cut cells, possibly restricted to a list of cells
......@@ -80,22 +79,23 @@ class cuttingPlane
(
const primitiveMesh&,
const scalarField& dotProducts,
const labelList& cellIdLabels = labelList::null()
const UList<label>& cellIdLabels = UList<label>::null()
);
//- Determine intersection points (cutPoints_).
labelList intersectEdges
//- Determine intersection points (cutPoints).
void intersectEdges
(
const primitiveMesh&,
const scalarField& dotProducts
const scalarField& dotProducts,
List<label>& edgePoint
);
//- Walk around circumference of cell starting from startEdgeI crossing
//- Walk circumference of cell, starting from startEdgeI crossing
// only cut edges. Record cutPoint labels in faceVerts.
static bool walkCell
(
const primitiveMesh&,
const labelList& edgePoint,
const UList<label>& edgePoint,
const label cellI,
const label startEdgeI,
DynamicList<label>& faceVerts
......@@ -105,7 +105,7 @@ class cuttingPlane
void walkCellCuts
(
const primitiveMesh& mesh,
const labelList& edgePoint
const UList<label>& edgePoint
);
......@@ -116,53 +116,51 @@ protected:
//- Construct plane description without cutting
cuttingPlane(const plane&);
// Protected Member Functions
//- recut mesh with existing planeDesc
void reCut(const primitiveMesh&);
//- recut mesh with existing planeDesc, restricted to a list of cells
void reCut
(
const primitiveMesh&,
const labelList& cellIdLabels
const UList<label>& cellIdLabels = UList<label>::null()
);
//- remap action on triangulation or cleanup
virtual void remapFaces(const UList<label>& faceMap);
public:
// Constructors
//- Construct from components: Mesh reference and plane
cuttingPlane(const primitiveMesh&, const plane&);
//- Construct from mesh reference and plane,
// restricted to a list of cells
// possibly restricted to a list of cells
cuttingPlane
(
const primitiveMesh&,
const plane&,
const labelList& cellIdLabels
const UList<label>& cellIdLabels = UList<label>::null()
);
// Member Functions
//- Return plane used
const plane& planeDesc() const;
//- Return pointField of cutting points
const pointField& points() const;
//- Return faceList of points in cells
const faceList& faces() const;
const plane& planeDesc() const
{
return static_cast<const plane&>(*this);
}
//- Return labelList of cut cells
const labelList& cells() const;
//- Return List of cells cut by the plane
const labelList& cutCells() const
{
return cutCells_;
}
//- Return true or false to question: have any cells been cut?
bool cut();
bool cut() const
{
return (cutCells_.size() != 0);
}
//- Sample the cell field
template<class Type>
......@@ -185,7 +183,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "cuttingPlaneSample.C"
# include "cuttingPlaneTemplates.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -39,7 +39,7 @@ Foam::tmp<Foam::Field<Type> > Foam::cuttingPlane::sample
const Field<Type>& sf
) const
{
return tmp<Field<Type> >(new Field<Type>(sf, cells()));
return tmp<Field<Type> >(new Field<Type>(sf, cutCells()));
}
......
......@@ -44,55 +44,30 @@ namespace Foam
void Foam::sampledPatch::createGeometry()
{
clearGeom();
points_.clear();
faces_.clear();
sampledSurface::clearGeom();
MeshStorage::clear();
patchFaceLabels_.clear();
if (patchIndex() != -1)
if (patchIndex_ != -1)
{
const polyPatch& patch = mesh().boundaryMesh()[patchIndex()];
const faceList& localFaces = patch.localFaces();
const polyPatch& p = mesh().boundaryMesh()[patchIndex_];
this->storedPoints() = p.localPoints();
this->storedFaces() = p.localFaces();
points_ = patch.localPoints();
if (triangulate())
// an identity map
patchFaceLabels_.setSize(faces().size());
forAll(patchFaceLabels_, i)
{
// Count triangles
label nTri = 0;
forAll(localFaces, faceI)
{
const face& f = localFaces[faceI];
nTri += f.nTriangles(patch.localPoints());
}
faces_.setSize(nTri);
patchFaceLabels_.setSize(nTri);
// split and fill mesh face references
nTri = 0;
forAll(localFaces, faceI)
{
const face& f = localFaces[faceI];
label fillIndex = nTri;
f.triangles(patch.localPoints(), nTri, faces_);
while (fillIndex < nTri)
{
patchFaceLabels_[fillIndex++] = faceI;
}
}
patchFaceLabels_[i] = i;
}
else
// triangulate uses remapFaces()
// - this is somewhat less efficient since it recopies the faces
// that we just created, but we probably don't want to do this
// too often anyhow.
if (triangulate_)
{
faces_ = localFaces;
patchFaceLabels_.setSize(faces_.size());
forAll(localFaces, i)
{
patchFaceLabels_[i] = i;
}
MeshStorage::triangulate();
}
}
......@@ -114,11 +89,10 @@ Foam::sampledPatch::sampledPatch
const bool triangulate
)
:
sampledSurface(name, mesh, triangulate),
sampledSurface(name, mesh),
patchName_(patchName),
patchIndex_(mesh.boundaryMesh().findPatchID(patchName_)),
points_(0),
faces_(0),
triangulate_(triangulate),
patchFaceLabels_(0)
{
createGeometry();
......@@ -135,12 +109,9 @@ Foam::sampledPatch::sampledPatch
sampledSurface(name, mesh, dict),
patchName_(dict.lookup("patchName")),
patchIndex_(mesh.boundaryMesh().findPatchID(patchName_)),
points_(0),
faces_(0),
triangulate_(dict.lookupOrDefault("triangulate", false)),
patchFaceLabels_(0)
{
// default: non-triangulated
triangulate() = dict.lookupOrDefault("triangulate", false);
createGeometry();
}
......@@ -162,6 +133,28 @@ void Foam::sampledPatch::correct(const bool meshChanged)
}
// remap action on triangulation
void Foam::sampledPatch::remapFaces
(
const UList<label>& faceMap
)
{