diff --git a/applications/utilities/postProcessing/sampling/sample/sampleDict b/applications/utilities/postProcessing/sampling/sample/sampleDict
index 92b1c02e35b983ca4e52f377cdc6bcf2ae5fec76..88c86bf482eda4b44f3d96bca85e4ed6468cebe5 100644
--- a/applications/utilities/postProcessing/sampling/sample/sampleDict
+++ b/applications/utilities/postProcessing/sampling/sample/sampleDict
@@ -31,18 +31,66 @@ setFormat raw;
// dx : DX scalar or vector format
// vtk : VTK ascii format
// raw : x y z value format for use with e.g. gnuplot 'splot'.
+// boundaryData: written in a form that can be used with timeVaryingMapped
+// boundary conditions (move to constant/boundaryData)
+// starcd : Star-CD format (geometry in .cel/.vrt, data in .usr)
+// nastran : Nastran (.nas) format
//
// Note:
-// other formats such as obj, stl, etc can also be written (by proxy)
-// but without any values!
+// - other formats such as obj, stl, etc can also be written (by proxy)
+// but without any values!
+// - boundaryData format: typical use:
+// - use a surfaces functionObject to sample a patch:
+// type surfaces;
+// surfaceFormat boundaryData;
+// fields ( p );
+// surfaces
+// (
+// inlet
+// {
+// type patch;
+// patches (inpletPatch);
+// interpolate false;
+// }
+// );
+// - the boundaryData writer will write postProcessing/surfaces/inlet/
+// - move this to constant/boundaryData/inlet of the destination case
+// - use a timeVaryingMappedFixedValue bc to read&interpolate and use
+// as fixedValue.
+// - boundaryData: 'interpolate false' writes face centres, 'interpolate true'
+// writes points. For 2D case the face centres might only be a single column
+// so cannot be used in timeVaryingMapped with standard mapping (since
+// uses triangulation to do interpolation).
+
surfaceFormat vtk;
-// optionally define extra controls for the output formats
+// Optionally define extra controls for the output formats
formatOptions
{
ensight
{
+ // ascii/binary format
format ascii;
+ //collateTimes true; // write single file containing multiple timesteps
+ // (only for static surfaces)
+ }
+ nastran
+ {
+ // From OpenFOAM field name to Nastran field name
+ fields
+ (
+ (U PLOAD4)
+ (p PLOAD2)
+ );
+ // Optional scale
+ scale 1.0;
+ // Optional format
+ // short/long/free format
+ // short: scalar field width 8. default.
+ // long : scalar field width 16
+ // free : comma separated, field width according to writePrecision
+ // setting
+ format short; // short/long/free
}
}
@@ -143,9 +191,19 @@ sets
type patchSeed;
axis xyz;
patches (".*Wall.*");
- // Number of points to seed. Divided amongst all processors according
- // to fraction of patches they hold.
+
+ // Subset patch faces by:
+
+ // 1. Number of points to seed. Divided amongst all processors
+ // according to fraction of patches they hold.
maxPoints 100;
+
+ // 2. Specified set of locations. This selects for every the specified
+ // point the nearest patch face. (in addition the number of points
+ // is also truncated by the maxPoints setting)
+ // The difference with patchCloud is that this selects patch
+ // face centres, not an arbitrary location on the face.
+ points ((0.049 0.099 0.005)(0.051 0.054 0.005));
}
);
@@ -249,8 +307,10 @@ surfaces
isoValue 0.5;
interpolate true;
- //zone ABC; // Optional: zone only
- //exposedPatchName fixedWalls; // Optional: zone only
+ // zone ABC; // Optional: zone only
+ // exposedPatchName fixedWalls; // Optional: zone only
+
+ // bounds (1 1 1)(2 2 2); // Optional: limit extent
// regularise false; // Optional: do not simplify
// mergeTol 1e-10; // Optional: fraction of mesh bounding box
@@ -265,6 +325,9 @@ surfaces
isoValue 0.5;
interpolate false;
regularise false; // do not simplify
+
+ // bounds (1 1 1)(2 2 2); // Optional: limit extent
+
// mergeTol 1e-10; // Optional: fraction of mesh bounding box
// to merge points (default=1e-6)
}
@@ -281,8 +344,10 @@ surfaces
}
interpolate true;
- //zone ABC; // Optional: zone only
- //exposedPatchName fixedWalls; // Optional: zone only
+ // zone ABC; // Optional: zone only
+ // exposedPatchName fixedWalls; // Optional: zone only
+
+ // bounds (1 1 1)(2 2 2); // Optional: limit extent
// regularise false; // Optional: do not simplify
// mergeTol 1e-10; // Optional: fraction of mesh bounding box
@@ -305,6 +370,8 @@ surfaces
// of isoSurfaceCell
interpolate false;
regularise false; // Optional: do not simplify
+ // bounds (1 1 1)(2 2 2); // Optional: limit extent
+
// mergeTol 1e-10; // Optional: fraction of mesh bounding box
// to merge points (default=1e-6)
}
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C
new file mode 100644
index 0000000000000000000000000000000000000000..8085ccd8f0ca7093040c19ca2e9c3e444fd8c4ad
--- /dev/null
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C
@@ -0,0 +1,242 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | Copyright (C) 2015 OpenFOAM Foundation
+ \\/ 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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "triangle.H"
+#include "triPoints.H"
+#include "plane.H"
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+template<class Point, class PointRef>
+template<class AboveOp, class BelowOp>
+inline void Foam::triangle<Point, PointRef>::triSliceWithPlane
+(
+ const plane& pl,
+ const triPoints& tri,
+ AboveOp& aboveOp,
+ BelowOp& belowOp
+)
+{
+ // distance to cutting plane
+ FixedList<scalar, 3> d;
+
+ // determine how many of the points are above the cutting plane
+ label nPos = 0;
+ label posI = -1;
+ label negI = -1;
+ forAll(tri, i)
+ {
+ d[i] = ((tri[i] - pl.refPoint()) & pl.normal());
+
+ if (d[i] > 0)
+ {
+ nPos++;
+ posI = i;
+ }
+ else
+ {
+ negI = i;
+ }
+ }
+
+ if (nPos == 3)
+ {
+ aboveOp(tri);
+ }
+ else if (nPos == 2)
+ {
+ // point under the plane
+ label i0 = negI;
+
+ // indices of remaining points
+ label i1 = d.fcIndex(i0);
+ label i2 = d.fcIndex(i1);
+
+ // determine the two intersection points
+ point p01 = planeIntersection(d, tri, i0, i1);
+ point p02 = planeIntersection(d, tri, i0, i2);
+
+ aboveOp(triPoints(tri[i1], tri[i2], p02));
+ aboveOp(triPoints(tri[i1], p02, p01));
+ belowOp(triPoints(tri[i0], p01, p02));
+ }
+ else if (nPos == 1)
+ {
+ // point above the plane
+ label i0 = posI;
+
+ // indices of remaining points
+ label i1 = d.fcIndex(i0);
+ label i2 = d.fcIndex(i1);
+
+ // determine the two intersection points
+ point p01 = planeIntersection(d, tri, i0, i1);
+ point p02 = planeIntersection(d, tri, i0, i2);
+
+ belowOp(triPoints(tri[i1], tri[i2], p02));
+ belowOp(triPoints(tri[i1], p02, p01));
+ aboveOp(triPoints(tri[i0], p01, p02));
+ }
+ else
+ {
+ // All below
+ belowOp(tri);
+ }
+}
+
+
+template<class Point, class PointRef>
+template<class AboveOp, class BelowOp>
+inline void Foam::triangle<Point, PointRef>::sliceWithPlane
+(
+ const plane& pl,
+ AboveOp& aboveOp,
+ BelowOp& belowOp
+) const
+{
+ triSliceWithPlane(pl, triPoints(a_, b_, c_), aboveOp, belowOp);
+}
+
+
+template<class Point, class PointRef>
+template<class InsideOp, class OutsideOp>
+inline void Foam::triangle<Point, PointRef>::triangleOverlap
+(
+ const vector& n,
+ const triangle<Point, PointRef>& tgt,
+ InsideOp& insideOp,
+ OutsideOp& outsideOp
+) const
+{
+ // There are two possibilities with this algorithm - we either cut
+ // the outside triangles with all the edges or not (and keep them
+ // as disconnected triangles). In the first case
+ // we cannot do any evaluation short cut so we've chosen not to re-cut
+ // the outside triangles.
+
+
+ triIntersectionList insideTrisA;
+ label nInsideA = 0;
+ storeOp insideOpA(insideTrisA, nInsideA);
+
+ triIntersectionList outsideTrisA;
+ label nOutsideA = 0;
+ storeOp outsideOpA(outsideTrisA, nOutsideA);
+
+
+ const triPoints thisTri(a_, b_, c_);
+
+
+ // Cut original triangle with tgt edge 0.
+ // From *this to insideTrisA, outsideTrisA.
+ {
+ scalar s = Foam::mag(tgt.b() - tgt.a());
+ const plane pl0(tgt.a(), tgt.b(), tgt.b() + s*n);
+ triSliceWithPlane(pl0, thisTri, insideOpA, outsideOpA);
+ }
+
+ // Shortcut if nothing cut
+
+ if (insideOpA.nTris_ == 0)
+ {
+ outsideOp(thisTri);
+ return;
+ }
+
+ if (outsideOpA.nTris_ == 0)
+ {
+ insideOp(thisTri);
+ return;
+ }
+
+
+ // Cut all triangles with edge 1.
+ // From insideTrisA to insideTrisB, outsideTrisA
+
+ triIntersectionList insideTrisB;
+ label nInsideB = 0;
+ storeOp insideOpB(insideTrisB, nInsideB);
+
+ //triIntersectionList outsideTrisB;
+ //label nOutsideB = 0;
+ //storeOp outsideOpB(outsideTrisB, nOutsideB);
+
+ {
+ scalar s = Foam::mag(tgt.c() - tgt.b());
+ const plane pl0(tgt.b(), tgt.c(), tgt.c() + s*n);
+
+ for (label i = 0; i < insideOpA.nTris_; i++)
+ {
+ const triPoints& tri = insideOpA.tris_[i];
+ triSliceWithPlane(pl0, tri, insideOpB, outsideOpA);
+ }
+
+ //// Recut outside triangles (not necessary if only interested in
+ //// intersection properties)
+ //for (label i = 0; i < outsideOpA.nTris_; i++)
+ //{
+ // const triPoints& tri = outsideOpA.tris_[i];
+ // triSliceWithPlane(pl0, tri, outsideOpB, outsideOpB);
+ //}
+ }
+
+
+ // Cut all triangles with edge 2.
+ // From insideTrisB to insideTrisA, outsideTrisA
+ {
+ scalar s = Foam::mag(tgt.a() - tgt.c());
+ const plane pl0(tgt.c(), tgt.a(), tgt.a() + s*n);
+
+ insideOpA.nTris_ = 0;
+ //outsideOpA.nTris_ = 0;
+ for (label i = 0; i < insideOpB.nTris_; i++)
+ {
+ const triPoints& tri = insideOpB.tris_[i];
+ triSliceWithPlane(pl0, tri, insideOpA, outsideOpA);
+ }
+
+ //// Recut outside triangles (not necessary if only interested in
+ //// intersection properties)
+ //for (label i = 0; i < outsideOpB.nTris_; i++)
+ //{
+ // const triPoints& tri = outsideOpB.tris_[i];
+ // triSliceWithPlane(pl0, tri, outsideOpA, outsideOpA);
+ //}
+ }
+
+ // Transfer from A to argument
+ for (label i = 0; i < insideOpA.nTris_; i++)
+ {
+ insideOp(insideOpA.tris_[i]);
+ }
+
+ for (label i = 0; i < outsideOpA.nTris_; i++)
+ {
+ outsideOp(outsideOpA.tris_[i]);
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
index 6e8030c797658bb16108b4d32bad9766d3adaf7f..de477d20e05277b6c1f85660643e43d6ecac312a 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -52,6 +52,8 @@ namespace Foam
class Istream;
class Ostream;
+class triPoints;
+class plane;
// Forward declaration of friend functions and operators
@@ -71,6 +73,8 @@ inline Ostream& operator<<
const triangle<Point, PointRef>&
);
+typedef triangle<point, const point&> triPointRef;
+
/*---------------------------------------------------------------------------*\
Class triangle Declaration
@@ -79,27 +83,98 @@ inline Ostream& operator<<
template<class Point, class PointRef>
class triangle
{
+public:
+
+ // Public typedefs
+
+ //- Storage type for triangles originating from intersecting triangle
+ // with another triangle
+ typedef FixedList<triPoints, 27> triIntersectionList;
+
+ //- Return types for classify
+ enum proxType
+ {
+ NONE,
+ POINT, // Close to point
+ EDGE // Close to edge
+ };
+
+
+ // Public classes
+
+ // Classes for use in sliceWithPlane. What to do with decomposition
+ // of triangle.
+
+ //- Dummy
+ class dummyOp
+ {
+ public:
+ inline void operator()(const triPoints&);
+ };
+
+ //- Sum resulting areas
+ class sumAreaOp
+ {
+ public:
+ scalar area_;
+
+ inline sumAreaOp();
+
+ inline void operator()(const triPoints&);
+ };
+
+ //- Store resulting tris
+ class storeOp
+ {
+ public:
+ triIntersectionList& tris_;
+ label& nTris_;
+
+ inline storeOp(triIntersectionList&, label&);
+
+ inline void operator()(const triPoints&);
+ };
+
+
+private:
+
// Private data
PointRef a_, b_, c_;
-public:
+ // Private Member Functions
- //- Return types for classify
- enum proxType
- {
- NONE,
- POINT, // Close to point
- EDGE // Close to edge
- };
+ //- Helper: calculate intersection point
+ inline static point planeIntersection
+ (
+ const FixedList<scalar, 3>& d,
+ const triPoints& t,
+ const label negI,
+ const label posI
+ );
+
+ //- Helper: slice triangle with plane
+ template<class AboveOp, class BelowOp>
+ inline static void triSliceWithPlane
+ (
+ const plane& pl,
+ const triPoints& tri,
+ AboveOp& aboveOp,
+ BelowOp& belowOp
+ );
+public:
+
// Constructors
//- Construct from three points
inline triangle(const Point& a, const Point& b, const Point& c);
+ //- Construct from three points
+ inline triangle(const FixedList<Point, 3>&);
+
//- Construct from three points in the list of points
// The indices could be from triFace etc.
inline triangle
@@ -237,6 +312,27 @@ public:
pointHit& edgePoint
) const;
+ //- Decompose triangle into triangles above and below plane
+ template<class AboveOp, class BelowOp>
+ inline void sliceWithPlane
+ (
+ const plane& pl,
+ AboveOp& aboveOp,
+ BelowOp& belowOp
+ ) const;
+
+ //- Decompose triangle into triangles inside and outside
+ // (with respect to user provided normal) other
+ // triangle.
+ template<class InsideOp, class OutsideOp>
+ inline void triangleOverlap
+ (
+ const vector& n,
+ const triangle<Point, PointRef>& tri,
+ InsideOp& insideOp,
+ OutsideOp& outsideOp
+ ) const;
+
// IOstream operators
@@ -264,6 +360,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+#ifdef NoRepository
+# include "triangle.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
#endif
// ************************************************************************* //
diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
index fbc5a14286b8d5572dcfadff2d1c0280d4ebcde1..e3f8993e2eac9c4c13515c45defb3c810f6e2eaa 100644
--- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
+++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -25,6 +25,7 @@ License
#include "IOstreams.H"
#include "pointHit.H"
+#include "triPoints.H"
#include "mathematicalConstants.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
@@ -43,6 +44,18 @@ inline Foam::triangle<Point, PointRef>::triangle
{}
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::triangle
+(
+ const FixedList<Point, 3>& tri
+)
+:
+ a_(tri[0]),
+ b_(tri[1]),
+ c_(tri[2])
+{}
+
+
template<class Point, class PointRef>
inline Foam::triangle<Point, PointRef>::triangle
(
@@ -804,6 +817,66 @@ inline Foam::pointHit Foam::triangle<Point, PointRef>::nearestPoint
}
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::dummyOp::operator()
+(
+ const triPoints&
+)
+{}
+
+
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::sumAreaOp::sumAreaOp()
+:
+ area_(0.0)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::sumAreaOp::operator()
+(
+ const triPoints& tri
+)
+{
+ area_ += triangle<Point, const Point&>(tri).mag();
+}
+
+
+template<class Point, class PointRef>
+inline Foam::triangle<Point, PointRef>::storeOp::storeOp
+(
+ triIntersectionList& tris,
+ label& nTris
+)
+:
+ tris_(tris),
+ nTris_(nTris)
+{}
+
+
+template<class Point, class PointRef>
+inline void Foam::triangle<Point, PointRef>::storeOp::operator()
+(
+ const triPoints& tri
+)
+{
+ tris_[nTris_++] = tri;
+}
+
+
+template<class Point, class PointRef>
+inline Foam::point Foam::triangle<Point, PointRef>::planeIntersection
+(
+ const FixedList<scalar, 3>& d,
+ const triPoints& t,
+ const label negI,
+ const label posI
+)
+{
+ return (d[posI]*t[negI] - d[negI]*t[posI])/(-d[negI] + d[posI]);
+}
+
+
// * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class Point, class PointRef>
diff --git a/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H b/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H
index 9fea06cf2c758c665a84d42c6195b34ab0086265..b832f7e2bede3c8bf524a951b241d3c794f7d049 100644
--- a/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H
+++ b/src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H
@@ -33,8 +33,8 @@ SeeAlso
\*---------------------------------------------------------------------------*/
-#ifndef plane_H
-#define plane_H
+#ifndef planeExtrusion_H
+#define planeExtrusion_H
#include "linearNormal.H"
diff --git a/src/sampling/sampledSet/patchSeed/patchSeedSet.C b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
index 3560fe42046778b9e36d572a8c59bdd432d9ce36..9cdcc7d230c6b8e2e2117a58265fb55be6e7e436 100644
--- a/src/sampling/sampledSet/patchSeed/patchSeedSet.C
+++ b/src/sampling/sampledSet/patchSeed/patchSeedSet.C
@@ -30,9 +30,8 @@ License
#include "treeDataFace.H"
#include "Time.H"
#include "meshTools.H"
-//#include "Random.H"
-// For 'facePoint' helper function only
#include "mappedPatchBase.H"
+#include "indirectPrimitivePatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -85,19 +84,145 @@ void Foam::patchSeedSet::calcSamples
}
- label totalSize = returnReduce(sz, sumOp<label>());
+ if (!rndGenPtr_.valid())
+ {
+ rndGenPtr_.reset(new Random(0));
+ }
+ Random& rndGen = rndGenPtr_();
+
+
+ if (selectedLocations_.size())
+ {
+ DynamicList<label> newPatchFaces(patchFaces.size());
+
+ // Find the nearest patch face
+ {
+ // 1. All processors find nearest local patch face for all
+ // selectedLocations
+
+ // All the info for nearest. Construct to miss
+ List<mappedPatchBase::nearInfo> nearest(selectedLocations_.size());
+
+ const indirectPrimitivePatch pp
+ (
+ IndirectList<face>(mesh().faces(), patchFaces),
+ mesh().points()
+ );
+
+ treeBoundBox patchBb
+ (
+ treeBoundBox(pp.points(), pp.meshPoints()).extend
+ (
+ rndGen,
+ 1e-4
+ )
+ );
+ patchBb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+ patchBb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
+
+ indexedOctree<treeDataFace> boundaryTree
+ (
+ treeDataFace // all information needed to search faces
+ (
+ false, // do not cache bb
+ mesh(),
+ patchFaces // boundary faces only
+ ),
+ patchBb, // overall search domain
+ 8, // maxLevel
+ 10, // leafsize
+ 3.0 // duplicity
+ );
+
+ // Get some global dimension so all points are equally likely
+ // to be found
+ const scalar globalDistSqr
+ (
+ //magSqr
+ //(
+ // boundBox
+ // (
+ // pp.points(),
+ // pp.meshPoints(),
+ // true
+ // ).span()
+ //)
+ GREAT
+ );
+
+ forAll(selectedLocations_, sampleI)
+ {
+ const point& sample = selectedLocations_[sampleI];
+
+ pointIndexHit& nearInfo = nearest[sampleI].first();
+ nearInfo = boundaryTree.findNearest
+ (
+ sample,
+ globalDistSqr
+ );
+
+ if (!nearInfo.hit())
+ {
+ nearest[sampleI].second().first() = Foam::sqr(GREAT);
+ nearest[sampleI].second().second() =
+ Pstream::myProcNo();
+ }
+ else
+ {
+ point fc(pp[nearInfo.index()].centre(pp.points()));
+ nearInfo.setPoint(fc);
+ nearest[sampleI].second().first() = magSqr(fc-sample);
+ nearest[sampleI].second().second() =
+ Pstream::myProcNo();
+ }
+ }
+
+
+ // 2. Reduce on master. Select nearest processor.
+
+ // Find nearest. Combine on master.
+ Pstream::listCombineGather(nearest, mappedPatchBase::nearestEqOp());
+ Pstream::listCombineScatter(nearest);
+
+
+ // 3. Pick up my local faces that have won
+
+ forAll(nearest, sampleI)
+ {
+ if (nearest[sampleI].first().hit())
+ {
+ label procI = nearest[sampleI].second().second();
+ label index = nearest[sampleI].first().index();
+
+ if (procI == Pstream::myProcNo())
+ {
+ newPatchFaces.append(pp.addressing()[index]);
+ }
+ }
+ }
+ }
+
+ if (debug)
+ {
+ Pout<< "Found " << newPatchFaces.size()
+ << " out of " << selectedLocations_.size()
+ << " on local processor" << endl;
+ }
+
+ patchFaces.transfer(newPatchFaces);
+ }
// Shuffle and truncate if in random mode
+ label totalSize = returnReduce(patchFaces.size(), sumOp<label>());
+
if (maxPoints_ < totalSize)
{
// Check what fraction of maxPoints_ I need to generate locally.
- label myMaxPoints = label(scalar(sz)/totalSize*maxPoints_);
+ label myMaxPoints =
+ label(scalar(patchFaces.size())/totalSize*maxPoints_);
- rndGenPtr_.reset(new Random(123456));
- Random& rndGen = rndGenPtr_();
-
- labelList subset = identity(sz);
+ labelList subset = identity(patchFaces.size());
for (label iter = 0; iter < 4; iter++)
{
forAll(subset, i)
@@ -115,7 +240,7 @@ void Foam::patchSeedSet::calcSamples
if (debug)
{
Pout<< "In random mode : selected " << patchFaces.size()
- << " faces out of " << sz << endl;
+ << " faces out of " << patchFaces.size() << endl;
}
}
@@ -135,6 +260,9 @@ void Foam::patchSeedSet::calcSamples
forAll(patchFaces, i)
{
label faceI = patchFaces[i];
+
+ // Slightly shift point in since on warped face face-diagonal
+ // decomposition might be outside cell for face-centre decomposition!
pointIndexHit info = mappedPatchBase::facePoint
(
mesh(),
@@ -217,9 +345,15 @@ Foam::patchSeedSet::patchSeedSet
wordReList(dict.lookup("patches"))
)
),
- //searchDist_(readScalar(dict.lookup("maxDistance"))),
- //offsetDist_(readScalar(dict.lookup("offsetDist"))),
- maxPoints_(readLabel(dict.lookup("maxPoints")))
+ maxPoints_(readLabel(dict.lookup("maxPoints"))),
+ selectedLocations_
+ (
+ dict.lookupOrDefault<pointField>
+ (
+ "points",
+ pointField(0)
+ )
+ )
{
genSamples();
diff --git a/src/sampling/sampledSet/patchSeed/patchSeedSet.H b/src/sampling/sampledSet/patchSeed/patchSeedSet.H
index a6e795ad2be1a9d641efb2cf9418b757792fcfc4..c2143a209852416839b8b6b7b4525dd986601e1f 100644
--- a/src/sampling/sampledSet/patchSeed/patchSeedSet.H
+++ b/src/sampling/sampledSet/patchSeed/patchSeedSet.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -58,18 +58,15 @@ class patchSeedSet
//- Patches to sample
const labelHashSet patchSet_;
-// //- Maximum distance to look for nearest
-// const scalar searchDist_;
-//
-// //- Offset distance
-// const scalar offsetDist_;
-
- //- Maximum number of patch faces to seed
+ //- Maximum number of patch faces to seed (if in random subset mode)
const label maxPoints_;
//- Random number generator (if maxPoints < num patch faces)
autoPtr<Random> rndGenPtr_;
+ //- Patch faces to seed selected based on nearness to supplied points
+ const pointField selectedLocations_;
+
// Private Member Functions
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
index 2dafe3e98ad268312e11ca50d68c84a2669b83ec..fd2927ca719e7225fb8e55cfc1c5e72c640745fc 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -278,7 +278,8 @@ void Foam::distanceSurface::createGeometry()
cellDistance,
pointDistance_,
distance_,
- regularise_
+ regularise_,
+ bounds_
)
);
}
@@ -291,7 +292,8 @@ void Foam::distanceSurface::createGeometry()
cellDistance,
pointDistance_,
distance_,
- regularise_
+ regularise_,
+ bounds_
)
);
}
@@ -336,6 +338,7 @@ Foam::distanceSurface::distanceSurface
cell_(dict.lookupOrDefault("cell", true)),
regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", false)),
+ bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
zoneKey_(keyType::null),
needsUpdate_(true),
isoSurfCellPtr_(NULL),
@@ -364,7 +367,8 @@ Foam::distanceSurface::distanceSurface
const bool signedDistance,
const bool cell,
const Switch regularise,
- const Switch average
+ const Switch average,
+ const boundBox& bounds
)
:
sampledSurface(name, mesh, interpolate),
@@ -390,6 +394,7 @@ Foam::distanceSurface::distanceSurface
cell_(cell),
regularise_(regularise),
average_(average),
+ bounds_(bounds),
zoneKey_(keyType::null),
needsUpdate_(true),
isoSurfCellPtr_(NULL),
diff --git a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
index 73cd555e05e35a04dee8c5854745d8ba36fcbfed..aa3f63251144fa3a46ba4140c416070d87b0f3df 100644
--- a/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
+++ b/src/sampling/sampledSurface/distanceSurface/distanceSurface.H
@@ -75,6 +75,9 @@ class distanceSurface
//- Whether to recalculate cell values as average of point values
const Switch average_;
+ //- Optional bounding box to trim triangles against
+ const boundBox bounds_;
+
//- If restricted to zones, name of this zone or a regular expression
keyType zoneKey_;
@@ -144,7 +147,8 @@ public:
const bool signedDistance,
const bool cell,
const Switch regularise,
- const Switch average
+ const Switch average,
+ const boundBox& bounds = boundBox::greatBox
);
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.C b/src/sampling/sampledSurface/isoSurface/isoSurface.C
index e3105b35a56410fe8c6c1ee73425ac45874fcb95..9dbe88440d49a2e43bdbcff0f96d2f24e34685bb 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -32,12 +32,36 @@ License
#include "surfaceFields.H"
#include "OFstream.H"
#include "meshTools.H"
+#include "triSurfaceSearch.H"
+#include "surfaceIntersection.H"
+#include "intersectedSurface.H"
+#include "searchableBox.H"
+#include "triSurfaceMesh.H"
+#include "triPoints.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
namespace Foam
{
defineTypeNameAndDebug(isoSurface, 0);
+
+// Helper class for slicing triangles
+class storeOp
+{
+public:
+ DynamicList<triPoints>& tris_;
+
+ inline storeOp(DynamicList<triPoints>& tris)
+ :
+ tris_(tris)
+ {}
+
+ inline void operator()(const triPoints& tri)
+ {
+ tris_.append(tri);
+ }
+};
+
}
@@ -473,43 +497,6 @@ void Foam::isoSurface::calcCutTypes
}
-// Return the two common points between two triangles
-Foam::labelPair Foam::isoSurface::findCommonPoints
-(
- const labelledTri& tri0,
- const labelledTri& tri1
-)
-{
- labelPair common(-1, -1);
-
- label fp0 = 0;
- label fp1 = findIndex(tri1, tri0[fp0]);
-
- if (fp1 == -1)
- {
- fp0 = 1;
- fp1 = findIndex(tri1, tri0[fp0]);
- }
-
- if (fp1 != -1)
- {
- // So tri0[fp0] is tri1[fp1]
-
- // Find next common point
- label fp0p1 = tri0.fcIndex(fp0);
- label fp1p1 = tri1.fcIndex(fp1);
- label fp1m1 = tri1.rcIndex(fp1);
-
- if (tri0[fp0p1] == tri1[fp1p1] || tri0[fp0p1] == tri1[fp1m1])
- {
- common[0] = tri0[fp0];
- common[1] = tri0[fp0p1];
- }
- }
- return common;
-}
-
-
// Caculate centre of surface.
Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
{
@@ -523,73 +510,6 @@ Foam::point Foam::isoSurface::calcCentre(const triSurface& s)
}
-// Replace surface (localPoints, localTris) with single point. Returns
-// point. Destructs arguments.
-Foam::pointIndexHit Foam::isoSurface::collapseSurface
-(
- pointField& localPoints,
- DynamicList<labelledTri, 64>& localTris
-)
-{
- pointIndexHit info(false, vector::zero, localTris.size());
-
- if (localTris.size() == 1)
- {
- const labelledTri& tri = localTris[0];
- info.setPoint(tri.centre(localPoints));
- info.setHit();
- }
- else if (localTris.size() == 2)
- {
- // Check if the two triangles share an edge.
- const labelledTri& tri0 = localTris[0];
- const labelledTri& tri1 = localTris[0];
-
- labelPair shared = findCommonPoints(tri0, tri1);
-
- if (shared[0] != -1)
- {
- info.setPoint
- (
- 0.5
- * (
- tri0.centre(localPoints)
- + tri1.centre(localPoints)
- )
- );
- info.setHit();
- }
- }
- else if (localTris.size())
- {
- // Check if single region. Rare situation.
- triSurface surf
- (
- localTris,
- geometricSurfacePatchList(0),
- localPoints,
- true
- );
- localTris.clearStorage();
-
- labelList faceZone;
- label nZones = surf.markZones
- (
- boolList(surf.nEdges(), false),
- faceZone
- );
-
- if (nZones == 1)
- {
- info.setPoint(calcCentre(surf));
- info.setHit();
- }
- }
-
- return info;
-}
-
-
// Determine per cell centre whether all the intersections get collapsed
// to a single point
void Foam::isoSurface::calcSnappedCc
@@ -1133,505 +1053,288 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
}
-// Does face use valid vertices?
-bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
+void Foam::isoSurface::trimToPlanes
+(
+ const PtrList<plane>& planes,
+ const triPointRef& tri,
+ DynamicList<point>& newTriPoints
+)
{
- // Simple check on indices ok.
-
- const labelledTri& f = surf[faceI];
+ // Buffer for generated triangles
+ DynamicList<triPoints> insideTrisA;
+ storeOp insideOpA(insideTrisA);
- if
- (
- (f[0] < 0) || (f[0] >= surf.points().size())
- || (f[1] < 0) || (f[1] >= surf.points().size())
- || (f[2] < 0) || (f[2] >= surf.points().size())
- )
- {
- WarningIn("validTri(const triSurface&, const label)")
- << "triangle " << faceI << " vertices " << f
- << " uses point indices outside point range 0.."
- << surf.points().size()-1 << endl;
+ // Buffer for generated triangles
+ DynamicList<triPoints> insideTrisB;
+ storeOp insideOpB(insideTrisB);
- return false;
- }
+ triPointRef::dummyOp dop;
- if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
- {
- WarningIn("validTri(const triSurface&, const label)")
- << "triangle " << faceI
- << " uses non-unique vertices " << f
- << endl;
- return false;
- }
+ // Store starting triangle in insideTrisA
+ insideOpA(triPoints(tri.a(), tri.b(), tri.c()));
- // duplicate triangle check
- const labelList& fFaces = surf.faceFaces()[faceI];
+ bool useA = true;
- // Check if faceNeighbours use same points as this face.
- // Note: discards normal information - sides of baffle are merged.
- forAll(fFaces, i)
+ forAll(planes, faceI)
{
- label nbrFaceI = fFaces[i];
+ const plane& pl = planes[faceI];
- if (nbrFaceI <= faceI)
+ if (useA)
{
- // lower numbered faces already checked
- continue;
+ insideOpB.tris_.clear();
+ forAll(insideOpA.tris_, i)
+ {
+ const triPoints& tri = insideOpA.tris_[i];
+ triPointRef(tri).sliceWithPlane(pl, insideOpB, dop);
+ }
+ }
+ else
+ {
+ insideOpA.tris_.clear();
+ forAll(insideOpB.tris_, i)
+ {
+ const triPoints& tri = insideOpB.tris_[i];
+ triPointRef(tri).sliceWithPlane(pl, insideOpA, dop);
+ }
}
+ useA = !useA;
+ }
- const labelledTri& nbrF = surf[nbrFaceI];
- if
- (
- ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
- && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
- && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
- )
+ // Transfer
+ if (useA)
+ {
+ forAll(insideOpA.tris_, i)
{
- WarningIn("validTri(const triSurface&, const label)")
- << "triangle " << faceI << " vertices " << f
- << " fc:" << f.centre(surf.points())
- << " has the same vertices as triangle " << nbrFaceI
- << " vertices " << nbrF
- << " fc:" << nbrF.centre(surf.points())
- << endl;
-
- return false;
+ const triPoints& tri = insideOpA.tris_[i];
+ newTriPoints.append(tri[0]);
+ newTriPoints.append(tri[1]);
+ newTriPoints.append(tri[2]);
+ }
+ }
+ else
+ {
+ forAll(insideOpB.tris_, i)
+ {
+ const triPoints& tri = insideOpB.tris_[i];
+ newTriPoints.append(tri[0]);
+ newTriPoints.append(tri[1]);
+ newTriPoints.append(tri[2]);
}
}
- return true;
}
-void Foam::isoSurface::calcAddressing
+void Foam::isoSurface::trimToBox
(
- const triSurface& surf,
- List<FixedList<label, 3> >& faceEdges,
- labelList& edgeFace0,
- labelList& edgeFace1,
- Map<labelList>& edgeFacesRest
-) const
+ const treeBoundBox& bb,
+ DynamicList<point>& triPoints,
+ DynamicList<label>& triMap // trimmed to original tris
+)
{
- const pointField& points = surf.points();
-
- pointField edgeCentres(3*surf.size());
- label edgeI = 0;
- forAll(surf, triI)
- {
- const labelledTri& tri = surf[triI];
- edgeCentres[edgeI++] = 0.5*(points[tri[0]]+points[tri[1]]);
- edgeCentres[edgeI++] = 0.5*(points[tri[1]]+points[tri[2]]);
- edgeCentres[edgeI++] = 0.5*(points[tri[2]]+points[tri[0]]);
- }
-
- pointField mergedCentres;
- labelList oldToMerged;
- bool hasMerged = mergePoints
- (
- edgeCentres,
- mergeDistance_,
- false,
- oldToMerged,
- mergedCentres
- );
-
if (debug)
{
- Pout<< "isoSurface : detected "
- << mergedCentres.size()
- << " geometric edges on " << surf.size() << " triangles." << endl;
- }
-
- if (!hasMerged)
- {
- return;
+ Pout<< "isoSurface : trimming to " << bb << endl;
}
-
- // Determine faceEdges
- faceEdges.setSize(surf.size());
- edgeI = 0;
- forAll(surf, triI)
+ // Generate inwards pointing planes
+ PtrList<plane> planes(6);
+ const pointField pts(bb.treeBoundBox::points());
+ forAll(treeBoundBox::faces, faceI)
{
- faceEdges[triI][0] = oldToMerged[edgeI++];
- faceEdges[triI][1] = oldToMerged[edgeI++];
- faceEdges[triI][2] = oldToMerged[edgeI++];
+ const face& f = treeBoundBox::faces[faceI];
+ const vector& n = treeBoundBox::faceNormals[faceI];
+ planes.set(faceI, new plane(pts[f[0]], -n));
}
+ label nTris = triPoints.size()/3;
- // Determine edgeFaces
- edgeFace0.setSize(mergedCentres.size());
- edgeFace0 = -1;
- edgeFace1.setSize(mergedCentres.size());
- edgeFace1 = -1;
- edgeFacesRest.clear();
-
- // Overflow edge faces for geometric shared edges that turned
- // out to be different anyway.
- EdgeMap<labelList> extraEdgeFaces(mergedCentres.size()/100);
+ DynamicList<point> newTriPoints(triPoints.size()/16);
+ triMap.setCapacity(nTris/16);
- forAll(oldToMerged, oldEdgeI)
+ label vertI = 0;
+ for (label triI = 0; triI < nTris; triI++)
{
- label triI = oldEdgeI / 3;
- label edgeI = oldToMerged[oldEdgeI];
+ const point& p0 = triPoints[vertI++];
+ const point& p1 = triPoints[vertI++];
+ const point& p2 = triPoints[vertI++];
+
+ label oldNPoints = newTriPoints.size();
+ trimToPlanes
+ (
+ planes,
+ triPointRef(p0, p1, p2),
+ newTriPoints
+ );
- if (edgeFace0[edgeI] == -1)
+ label nCells = (newTriPoints.size()-oldNPoints)/3;
+ for (label i = 0; i < nCells; i++)
{
- // First triangle for edge
- edgeFace0[edgeI] = triI;
+ triMap.append(triI);
}
- else
- {
- //- Check that the two triangles actually topologically
- // share an edge
- const labelledTri& prevTri = surf[edgeFace0[edgeI]];
- const labelledTri& tri = surf[triI];
+ }
- label fp = oldEdgeI % 3;
+ if (debug)
+ {
+ Pout<< "isoSurface : trimmed from " << nTris
+ << " down to " << triMap.size()
+ << " triangles." << endl;
+ }
- edge e(tri[fp], tri[tri.fcIndex(fp)]);
+ triPoints.transfer(newTriPoints);
+}
- label prevTriIndex = -1;
- forAll(prevTri, i)
- {
- if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
- {
- prevTriIndex = i;
- break;
- }
- }
+void Foam::isoSurface::trimToBox
+(
+ const treeBoundBox& bb,
+ DynamicList<point>& triPoints, // new points
+ DynamicList<label>& triMap, // map from (new) triangle to original
+ labelList& triPointMap, // map from (new) point to original
+ labelList& interpolatedPoints, // labels of newly introduced points
+ List<FixedList<label, 3> >& interpolatedOldPoints,// and their interpolation
+ List<FixedList<scalar, 3> >& interpolationWeights
+)
+{
+ const List<point> oldTriPoints(triPoints);
- if (prevTriIndex == -1)
- {
- // Different edge. Store for later.
- EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
+ // Trim triPoints, return map
+ trimToBox(bb, triPoints, triMap);
- if (iter != extraEdgeFaces.end())
- {
- labelList& eFaces = iter();
- label sz = eFaces.size();
- eFaces.setSize(sz+1);
- eFaces[sz] = triI;
- }
- else
- {
- extraEdgeFaces.insert(e, labelList(1, triI));
- }
- }
- else if (edgeFace1[edgeI] == -1)
- {
- edgeFace1[edgeI] = triI;
- }
- else
- {
- //WarningIn("orientSurface(triSurface&)")
- // << "Edge " << edgeI << " with centre "
- // << mergedCentres[edgeI]
- // << " used by more than two triangles: "
- // << edgeFace0[edgeI] << ", "
- // << edgeFace1[edgeI] << " and " << triI << endl;
- Map<labelList>::iterator iter = edgeFacesRest.find(edgeI);
-
- if (iter != edgeFacesRest.end())
- {
- labelList& eFaces = iter();
- label sz = eFaces.size();
- eFaces.setSize(sz+1);
- eFaces[sz] = triI;
- }
- else
- {
- edgeFacesRest.insert(edgeI, labelList(1, triI));
- }
- }
- }
- }
- // Add extraEdgeFaces
- edgeI = edgeFace0.size();
+ // Find point correspondence:
+ // - one-to-one for preserved original points (triPointMap)
+ // - interpolation for newly introduced points
+ // (interpolatedOldPoints)
+ label sz = oldTriPoints.size()/100;
+ DynamicList<label> dynInterpolatedPoints(sz);
+ DynamicList<FixedList<label, 3> > dynInterpolatedOldPoints(sz);
+ DynamicList<FixedList<scalar, 3> > dynInterpolationWeights(sz);
- edgeFace0.setSize(edgeI + extraEdgeFaces.size());
- edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
- forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
+ triPointMap.setSize(triPoints.size());
+ forAll(triMap, triI)
{
- const labelList& eFaces = iter();
+ label oldTriI = triMap[triI];
- // The current edge will become edgeI. Replace all occurrences in
- // faceEdges
- forAll(eFaces, i)
+ // Find point correspondence. Assumes coordinate is bit-copy.
+ for (label i = 0; i < 3; i++)
{
- label triI = eFaces[i];
- const labelledTri& tri = surf[triI];
+ label pointI = 3*triI+i;
+ const point& pt = triPoints[pointI];
- FixedList<label, 3>& fEdges = faceEdges[triI];
- forAll(tri, fp)
+ // Compare to old-triangle's points
+ label matchPointI = -1;
+ for (label j = 0; j < 3; j++)
{
- edge e(tri[fp], tri[tri.fcIndex(fp)]);
-
- if (e == iter.key())
+ label oldPointI = 3*oldTriI+j;
+ if (pt == oldTriPoints[oldPointI])
{
- fEdges[fp] = edgeI;
+ matchPointI = oldPointI;
break;
}
}
- }
-
- // Add face to edgeFaces
-
- edgeFace0[edgeI] = eFaces[0];
-
- if (eFaces.size() >= 2)
- {
- edgeFace1[edgeI] = eFaces[1];
+ triPointMap[pointI] = matchPointI;
- if (eFaces.size() > 2)
+ // If new point: calculate and store interpolation
+ if (matchPointI == -1)
{
- edgeFacesRest.insert
- (
- edgeI,
- SubList<label>(eFaces, eFaces.size()-2, 2)
- );
+ dynInterpolatedPoints.append(pointI);
+
+ FixedList<label, 3> oldPoints;
+ oldPoints[0] = 3*oldTriI;
+ oldPoints[1] = 3*oldTriI+1;
+ oldPoints[2] = 3*oldTriI+2;
+ dynInterpolatedOldPoints.append(oldPoints);
+
+ triPointRef tri(oldTriPoints, oldPoints);
+ scalarList bary;
+ tri.barycentric(pt, bary);
+ FixedList<scalar, 3> weights;
+ weights[0] = bary[0];
+ weights[1] = bary[1];
+ weights[2] = bary[2];
+ dynInterpolationWeights.append(weights);
}
}
-
- edgeI++;
}
-}
-
-
-void Foam::isoSurface::walkOrientation
-(
- const triSurface& surf,
- const List<FixedList<label, 3> >& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const label seedTriI,
- labelList& flipState
-)
-{
- // Do walk for consistent orientation.
- DynamicList<label> changedFaces(surf.size());
-
- changedFaces.append(seedTriI);
-
- while (changedFaces.size())
- {
- DynamicList<label> newChangedFaces(changedFaces.size());
- forAll(changedFaces, i)
- {
- label triI = changedFaces[i];
- const labelledTri& tri = surf[triI];
- const FixedList<label, 3>& fEdges = faceEdges[triI];
-
- forAll(fEdges, fp)
- {
- label edgeI = fEdges[fp];
-
- // my points:
- label p0 = tri[fp];
- label p1 = tri[tri.fcIndex(fp)];
-
- label nbrI =
- (
- edgeFace0[edgeI] != triI
- ? edgeFace0[edgeI]
- : edgeFace1[edgeI]
- );
-
- if (nbrI != -1 && flipState[nbrI] == -1)
- {
- const labelledTri& nbrTri = surf[nbrI];
-
- // nbr points
- label nbrFp = findIndex(nbrTri, p0);
-
- if (nbrFp == -1)
- {
- FatalErrorIn("isoSurface::walkOrientation(..)")
- << "triI:" << triI
- << " tri:" << tri
- << " p0:" << p0
- << " p1:" << p1
- << " fEdges:" << fEdges
- << " edgeI:" << edgeI
- << " edgeFace0:" << edgeFace0[edgeI]
- << " edgeFace1:" << edgeFace1[edgeI]
- << " nbrI:" << nbrI
- << " nbrTri:" << nbrTri
- << abort(FatalError);
- }
-
-
- label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-
- bool sameOrientation = (p1 == nbrP1);
-
- if (flipState[triI] == 0)
- {
- flipState[nbrI] = (sameOrientation ? 0 : 1);
- }
- else
- {
- flipState[nbrI] = (sameOrientation ? 1 : 0);
- }
- newChangedFaces.append(nbrI);
- }
- }
- }
-
- changedFaces.transfer(newChangedFaces);
- }
+ interpolatedPoints.transfer(dynInterpolatedPoints);
+ interpolatedOldPoints.transfer(dynInterpolatedOldPoints);
+ interpolationWeights.transfer(dynInterpolationWeights);
}
-void Foam::isoSurface::orientSurface
-(
- triSurface& surf,
- const List<FixedList<label, 3> >& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const Map<labelList>& edgeFacesRest
-)
+// Does face use valid vertices?
+bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
{
- // -1 : unvisited
- // 0 : leave as is
- // 1 : flip
- labelList flipState(surf.size(), -1);
-
- label seedTriI = 0;
-
- while (true)
- {
- // Find first unvisited triangle
- for
- (
- ;
- seedTriI < surf.size() && flipState[seedTriI] != -1;
- seedTriI++
- )
- {}
-
- if (seedTriI == surf.size())
- {
- break;
- }
-
- // Note: Determine orientation of seedTriI?
- // for now assume it is ok
- flipState[seedTriI] = 0;
+ // Simple check on indices ok.
- walkOrientation
- (
- surf,
- faceEdges,
- edgeFace0,
- edgeFace1,
- seedTriI,
- flipState
- );
- }
+ const labelledTri& f = surf[faceI];
- // Do actual flipping
- surf.clearOut();
- forAll(surf, triI)
+ if
+ (
+ (f[0] < 0) || (f[0] >= surf.points().size())
+ || (f[1] < 0) || (f[1] >= surf.points().size())
+ || (f[2] < 0) || (f[2] >= surf.points().size())
+ )
{
- if (flipState[triI] == 1)
- {
- labelledTri tri(surf[triI]);
-
- surf[triI][0] = tri[0];
- surf[triI][1] = tri[2];
- surf[triI][2] = tri[1];
- }
- else if (flipState[triI] == -1)
- {
- FatalErrorIn
- (
- "isoSurface::orientSurface(triSurface&, const label)"
- ) << "problem" << abort(FatalError);
- }
- }
-}
-
+ WarningIn("validTri(const triSurface&, const label)")
+ << "triangle " << faceI << " vertices " << f
+ << " uses point indices outside point range 0.."
+ << surf.points().size()-1 << endl;
-// Checks if triangle is connected through edgeI only.
-bool Foam::isoSurface::danglingTriangle
-(
- const FixedList<label, 3>& fEdges,
- const labelList& edgeFace1
-)
-{
- label nOpen = 0;
- forAll(fEdges, i)
- {
- if (edgeFace1[fEdges[i]] == -1)
- {
- nOpen++;
- }
+ return false;
}
- if (nOpen == 1 || nOpen == 2 || nOpen == 3)
- {
- return true;
- }
- else
+ if ((f[0] == f[1]) || (f[0] == f[2]) || (f[1] == f[2]))
{
+ WarningIn("validTri(const triSurface&, const label)")
+ << "triangle " << faceI
+ << " uses non-unique vertices " << f
+ << endl;
return false;
}
-}
+ // duplicate triangle check
-// Mark triangles to keep. Returns number of dangling triangles.
-Foam::label Foam::isoSurface::markDanglingTriangles
-(
- const List<FixedList<label, 3> >& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const Map<labelList>& edgeFacesRest,
- boolList& keepTriangles
-)
-{
- keepTriangles.setSize(faceEdges.size());
- keepTriangles = true;
-
- label nDangling = 0;
+ const labelList& fFaces = surf.faceFaces()[faceI];
- // Remove any dangling triangles
- forAllConstIter(Map<labelList>, edgeFacesRest, iter)
+ // Check if faceNeighbours use same points as this face.
+ // Note: discards normal information - sides of baffle are merged.
+ forAll(fFaces, i)
{
- // These are all the non-manifold edges. Filter out all triangles
- // with only one connected edge (= this edge)
-
- label edgeI = iter.key();
- const labelList& otherEdgeFaces = iter();
+ label nbrFaceI = fFaces[i];
- // Remove all dangling triangles
- if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
- {
- keepTriangles[edgeFace0[edgeI]] = false;
- nDangling++;
- }
- if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
+ if (nbrFaceI <= faceI)
{
- keepTriangles[edgeFace1[edgeI]] = false;
- nDangling++;
+ // lower numbered faces already checked
+ continue;
}
- forAll(otherEdgeFaces, i)
+
+ const labelledTri& nbrF = surf[nbrFaceI];
+
+ if
+ (
+ ((f[0] == nbrF[0]) || (f[0] == nbrF[1]) || (f[0] == nbrF[2]))
+ && ((f[1] == nbrF[0]) || (f[1] == nbrF[1]) || (f[1] == nbrF[2]))
+ && ((f[2] == nbrF[0]) || (f[2] == nbrF[1]) || (f[2] == nbrF[2]))
+ )
{
- label triI = otherEdgeFaces[i];
- if (danglingTriangle(faceEdges[triI], edgeFace1))
- {
- keepTriangles[triI] = false;
- nDangling++;
- }
+ WarningIn("validTri(const triSurface&, const label)")
+ << "triangle " << faceI << " vertices " << f
+ << " fc:" << f.centre(surf.points())
+ << " has the same vertices as triangle " << nbrFaceI
+ << " vertices " << nbrF
+ << " fc:" << nbrF.centre(surf.points())
+ << endl;
+
+ return false;
}
}
- return nDangling;
+ return true;
}
@@ -1715,6 +1418,7 @@ Foam::isoSurface::isoSurface
const scalarField& pVals,
const scalar iso,
const bool regularise,
+ const boundBox& bounds,
const scalar mergeTol
)
:
@@ -1722,6 +1426,7 @@ Foam::isoSurface::isoSurface
pVals_(pVals),
iso_(iso),
regularise_(regularise),
+ bounds_(bounds),
mergeDistance_(mergeTol*mesh_.bounds().mag())
{
if (debug)
@@ -1957,138 +1662,115 @@ Foam::isoSurface::isoSurface
}
-
- DynamicList<point> triPoints(nCutCells_);
- DynamicList<label> triMeshCells(nCutCells_);
-
- generateTriPoints
- (
- cValsPtr_(),
- pVals_,
-
- meshC,
- mesh_.points(),
-
- snappedPoints,
- snappedCc,
- snappedPoint,
-
- triPoints,
- triMeshCells
- );
-
- if (debug)
{
- Pout<< "isoSurface : generated " << triMeshCells.size()
- << " unmerged triangles from " << triPoints.size()
- << " unmerged points." << endl;
- }
+ DynamicList<point> triPoints(3*nCutCells_);
+ DynamicList<label> triMeshCells(nCutCells_);
-
- // Merge points and compact out non-valid triangles
- labelList triMap; // merged to unmerged triangle
- triSurface::operator=
- (
- stitchTriPoints
+ generateTriPoints
(
- true, // check for duplicate tris
- triPoints,
- triPointMergeMap_, // unmerged to merged point
- triMap
- )
- );
+ cValsPtr_(),
+ pVals_,
- if (debug)
- {
- Pout<< "isoSurface : generated " << triMap.size()
- << " merged triangles." << endl;
- }
+ meshC,
+ mesh_.points(),
- meshCells_.setSize(triMap.size());
- forAll(triMap, i)
- {
- meshCells_[i] = triMeshCells[triMap[i]];
- }
+ snappedPoints,
+ snappedCc,
+ snappedPoint,
- if (debug)
- {
- Pout<< "isoSurface : checking " << size()
- << " triangles for validity." << endl;
+ triPoints, // 3 points of the triangle
+ triMeshCells // per triangle the originating cell
+ );
- forAll(*this, triI)
+ if (debug)
{
- // Copied from surfaceCheck
- validTri(*this, triI);
+ Pout<< "isoSurface : generated " << triMeshCells.size()
+ << " unmerged triangles from " << triPoints.size()
+ << " unmerged points." << endl;
}
- }
-
-
- if (false)
- {
- List<FixedList<label, 3> > faceEdges;
- labelList edgeFace0, edgeFace1;
- Map<labelList> edgeFacesRest;
+ label nOldPoints = triPoints.size();
- while (true)
+ // Trimmed to original triangle
+ DynamicList<label> trimTriMap;
+ // Trimmed to original point
+ labelList trimTriPointMap;
+ if (bounds_ != boundBox::greatBox)
{
- // Calculate addressing
- calcAddressing
+ trimToBox
(
- *this,
- faceEdges,
- edgeFace0,
- edgeFace1,
- edgeFacesRest
+ treeBoundBox(bounds_),
+ triPoints, // new points
+ trimTriMap, // map from (new) triangle to original
+ trimTriPointMap, // map from (new) point to original
+ interpolatedPoints_, // labels of newly introduced points
+ interpolatedOldPoints_, // and their interpolation
+ interpolationWeights_
);
+ triMeshCells = labelField(triMeshCells, trimTriMap);
+ }
- // See if any dangling triangles
- boolList keepTriangles;
- label nDangling = markDanglingTriangles
+
+ // Merge points and compact out non-valid triangles
+ labelList triMap; // merged to unmerged triangle
+ triSurface::operator=
+ (
+ stitchTriPoints
(
- faceEdges,
- edgeFace0,
- edgeFace1,
- edgeFacesRest,
- keepTriangles
- );
+ true, // check for duplicate tris
+ triPoints,
+ triPointMergeMap_, // unmerged to merged point
+ triMap
+ )
+ );
- if (debug)
- {
- Pout<< "isoSurface : detected " << nDangling
- << " dangling triangles." << endl;
- }
+ if (debug)
+ {
+ Pout<< "isoSurface : generated " << triMap.size()
+ << " merged triangles." << endl;
+ }
- if (nDangling == 0)
- {
- break;
- }
- // Create face map (new to old)
- labelList subsetTriMap(findIndices(keepTriangles, true));
+ if (bounds_ != boundBox::greatBox)
+ {
+ // Adjust interpolatedPoints_
+ inplaceRenumber(triPointMergeMap_, interpolatedPoints_);
- labelList subsetPointMap;
- labelList reversePointMap;
- triSurface::operator=
- (
- subsetMesh
- (
- *this,
- subsetTriMap,
- reversePointMap,
- subsetPointMap
- )
- );
- meshCells_ = labelField(meshCells_, subsetTriMap);
- inplaceRenumber(reversePointMap, triPointMergeMap_);
+ // Adjust triPointMergeMap_
+ labelList newTriPointMergeMap(nOldPoints, -1);
+ forAll(trimTriPointMap, trimPointI)
+ {
+ label oldPointI = trimTriPointMap[trimPointI];
+ if (oldPointI >= 0)
+ {
+ label pointI = triPointMergeMap_[trimPointI];
+ if (pointI >= 0)
+ {
+ newTriPointMergeMap[oldPointI] = pointI;
+ }
+ }
+ }
+ triPointMergeMap_.transfer(newTriPointMergeMap);
}
- orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
+ meshCells_.setSize(triMap.size());
+ forAll(triMap, i)
+ {
+ meshCells_[i] = triMeshCells[triMap[i]];
+ }
}
-
if (debug)
{
+ Pout<< "isoSurface : checking " << size()
+ << " triangles for validity." << endl;
+
+ forAll(*this, triI)
+ {
+ // Copied from surfaceCheck
+ validTri(*this, triI);
+ }
+
fileName stlFile = mesh_.time().path() + ".stl";
Pout<< "Dumping surface to " << stlFile << endl;
triSurface::write(stlFile);
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurface.H b/src/sampling/sampledSurface/isoSurface/isoSurface.H
index 9f6dc7d76e0318abe34e3455c9ffeb9bd0de77b1..c74586a2ce1f4f1ee07aa9f59261582beef3339b 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurface.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -76,6 +76,8 @@ namespace Foam
{
class fvMesh;
+class plane;
+class treeBoundBox;
/*---------------------------------------------------------------------------*\
Class isoSurface Declaration
@@ -116,10 +118,12 @@ class isoSurface
//- Regularise?
const Switch regularise_;
+ //- Optional bounds
+ const boundBox bounds_;
+
//- When to merge points
const scalar mergeDistance_;
-
//- Whether face might be cut
List<cellCutType> faceCutType_;
@@ -135,6 +139,15 @@ class isoSurface
//- For every unmerged triangle point the point in the triSurface
labelList triPointMergeMap_;
+ //- triSurface points that have weighted interpolation
+ DynamicList<label> interpolatedPoints_;
+
+ //- corresponding original, unmerged points
+ DynamicList<FixedList<label, 3> > interpolatedOldPoints_;
+
+ //- corresponding weights
+ DynamicList<FixedList<scalar, 3> > interpolationWeights_;
+
// Private Member Functions
@@ -193,20 +206,8 @@ class isoSurface
const scalarField& pVals
);
- static labelPair findCommonPoints
- (
- const labelledTri&,
- const labelledTri&
- );
-
static point calcCentre(const triSurface&);
- static pointIndexHit collapseSurface
- (
- pointField& localPoints,
- DynamicList<labelledTri, 64>& localTris
- );
-
//- Determine per cc whether all near cuts can be snapped to single
// point.
void calcSnappedCc
@@ -323,6 +324,17 @@ class isoSurface
DynamicList<label>& triMeshCells
) const;
+ template<class Type>
+ static tmp<Field<Type> > interpolate
+ (
+ const label nPoints,
+ const labelList& triPointMergeMap,
+ const labelList& interpolatedPoints,
+ const List<FixedList<label, 3> >& interpolatedOldPoints,
+ const List<FixedList<scalar, 3> >& interpolationWeights,
+ const DynamicList<Type>& unmergedValues
+ );
+
triSurface stitchTriPoints
(
const bool checkDuplicates,
@@ -331,56 +343,37 @@ class isoSurface
labelList& triMap // merged to unmerged triangle
) const;
- //- Check single triangle for (topological) validity
- static bool validTri(const triSurface&, const label);
-
- //- Determine edge-face addressing
- void calcAddressing
- (
- const triSurface& surf,
- List<FixedList<label, 3> >& faceEdges,
- labelList& edgeFace0,
- labelList& edgeFace1,
- Map<labelList>& edgeFacesRest
- ) const;
-
- //- Determine orientation
- static void walkOrientation
+ //- Trim triangle to planes
+ static void trimToPlanes
(
- const triSurface& surf,
- const List<FixedList<label, 3> >& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const label seedTriI,
- labelList& flipState
+ const PtrList<plane>& planes,
+ const triPointRef& tri,
+ DynamicList<point>& newTriPoints
);
- //- Orient surface
- static void orientSurface
+ //- Trim all triangles to box
+ static void trimToBox
(
- triSurface&,
- const List<FixedList<label, 3> >& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const Map<labelList>& edgeFacesRest
+ const treeBoundBox& bb,
+ DynamicList<point>& triPoints,
+ DynamicList<label>& triMeshCells
);
- //- Is triangle (given by 3 edges) not fully connected?
- static bool danglingTriangle
+ //- Trim all triangles to box. Determine interpolation
+ // for existing and new points
+ static void trimToBox
(
- const FixedList<label, 3>& fEdges,
- const labelList& edgeFace1
+ const treeBoundBox& bb,
+ DynamicList<point>& triPoints,
+ DynamicList<label>& triMap,
+ labelList& triPointMap,
+ labelList& interpolatedPoints,
+ List<FixedList<label, 3> >& interpolatedOldPoints,
+ List<FixedList<scalar, 3> >& interpolationWeights
);
- //- Mark all non-fully connected triangles
- static label markDanglingTriangles
- (
- const List<FixedList<label, 3> >& faceEdges,
- const labelList& edgeFace0,
- const labelList& edgeFace1,
- const Map<labelList>& edgeFacesRest,
- boolList& keepTriangles
- );
+ //- Check single triangle for (topological) validity
+ static bool validTri(const triSurface&, const label);
static triSurface subsetMesh
(
@@ -392,6 +385,10 @@ class isoSurface
public:
+ //- Declare friendship with isoSurfaceCell to share some functionality
+ friend class isoSurfaceCell;
+
+
//- Runtime type information
TypeName("isoSurface");
@@ -407,24 +404,19 @@ public:
const scalarField& pointIsoVals,
const scalar iso,
const bool regularise,
+ const boundBox& bounds = boundBox::greatBox,
const scalar mergeTol = 1e-6 // fraction of bounding box
);
// Member Functions
- //- For every face original cell in mesh
+ //- For every triangle the original cell in mesh
const labelList& meshCells() const
{
return meshCells_;
}
- //- For every unmerged triangle point the point in the triSurface
- const labelList& triPointMergeMap() const
- {
- return triPointMergeMap_;
- }
-
//- Interpolates cCoords,pCoords. Uses the references to the original
// fields used to create the iso surface.
template<class Type>
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
index 8123541ebee8bbabcf84c9586a0bc523b5d90c6a..f0f140cfd4a791c96b59b2a6eb6334caf51cfc9e 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -30,6 +30,9 @@ License
#include "tetMatcher.H"
#include "syncTools.H"
#include "addToRunTimeSelectionTable.H"
+#include "Time.H"
+#include "triPoints.H"
+#include "isoSurface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -1222,144 +1225,6 @@ void Foam::isoSurfaceCell::calcAddressing
}
-//void Foam::isoSurfaceCell::walkOrientation
-//(
-// const triSurface& surf,
-// const List<FixedList<label, 3> >& faceEdges,
-// const labelList& edgeFace0,
-// const labelList& edgeFace1,
-// const label seedTriI,
-// labelList& flipState
-//)
-//{
-// // Do walk for consistent orientation.
-// DynamicList<label> changedFaces(surf.size());
-//
-// changedFaces.append(seedTriI);
-//
-// while (changedFaces.size())
-// {
-// DynamicList<label> newChangedFaces(changedFaces.size());
-//
-// forAll(changedFaces, i)
-// {
-// label triI = changedFaces[i];
-// const labelledTri& tri = surf[triI];
-// const FixedList<label, 3>& fEdges = faceEdges[triI];
-//
-// forAll(fEdges, fp)
-// {
-// label edgeI = fEdges[fp];
-//
-// // my points:
-// label p0 = tri[fp];
-// label p1 = tri[tri.fcIndex(fp)];
-//
-// label nbrI =
-// (
-// edgeFace0[edgeI] != triI
-// ? edgeFace0[edgeI]
-// : edgeFace1[edgeI]
-// );
-//
-// if (nbrI != -1 && flipState[nbrI] == -1)
-// {
-// const labelledTri& nbrTri = surf[nbrI];
-//
-// // nbr points
-// label nbrFp = findIndex(nbrTri, p0);
-// label nbrP1 = nbrTri[nbrTri.rcIndex(nbrFp)];
-//
-// bool sameOrientation = (p1 == nbrP1);
-//
-// if (flipState[triI] == 0)
-// {
-// flipState[nbrI] = (sameOrientation ? 0 : 1);
-// }
-// else
-// {
-// flipState[nbrI] = (sameOrientation ? 1 : 0);
-// }
-// newChangedFaces.append(nbrI);
-// }
-// }
-// }
-//
-// changedFaces.transfer(newChangedFaces);
-// }
-//}
-//
-//
-//void Foam::isoSurfaceCell::orientSurface
-//(
-// triSurface& surf,
-// const List<FixedList<label, 3> >& faceEdges,
-// const labelList& edgeFace0,
-// const labelList& edgeFace1,
-// const Map<labelList>& edgeFacesRest
-//)
-//{
-// // -1 : unvisited
-// // 0 : leave as is
-// // 1 : flip
-// labelList flipState(surf.size(), -1);
-//
-// label seedTriI = 0;
-//
-// while (true)
-// {
-// // Find first unvisited triangle
-// for
-// (
-// ;
-// seedTriI < surf.size() && flipState[seedTriI] != -1;
-// seedTriI++
-// )
-// {}
-//
-// if (seedTriI == surf.size())
-// {
-// break;
-// }
-//
-// // Note: Determine orientation of seedTriI?
-// // for now assume it is ok
-// flipState[seedTriI] = 0;
-//
-// walkOrientation
-// (
-// surf,
-// faceEdges,
-// edgeFace0,
-// edgeFace1,
-// seedTriI,
-// flipState
-// );
-// }
-//
-// // Do actual flipping
-// surf.clearOut();
-// forAll(surf, triI)
-// {
-// if (flipState[triI] == 1)
-// {
-// labelledTri tri(surf[triI]);
-//
-// surf[triI][0] = tri[0];
-// surf[triI][1] = tri[2];
-// surf[triI][2] = tri[1];
-// }
-// else if (flipState[triI] == -1)
-// {
-// FatalErrorIn
-// (
-// "isoSurfaceCell::orientSurface(triSurface&, const label)"
-// ) << "problem" << abort(FatalError);
-// }
-// }
-//}
-
-
// Checks if triangle is connected through edgeI only.
bool Foam::isoSurfaceCell::danglingTriangle
(
@@ -1517,6 +1382,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
const scalarField& pVals,
const scalar iso,
const bool regularise,
+ const boundBox& bounds,
const scalar mergeTol
)
:
@@ -1524,6 +1390,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
cVals_(cVals),
pVals_(pVals),
iso_(iso),
+ bounds_(bounds),
mergeDistance_(mergeTol*mesh.bounds().mag())
{
if (debug)
@@ -1607,57 +1474,105 @@ Foam::isoSurfaceCell::isoSurfaceCell
}
+ {
+ DynamicList<point> triPoints(nCutCells_);
+ DynamicList<label> triMeshCells(nCutCells_);
- DynamicList<point> triPoints(nCutCells_);
- DynamicList<label> triMeshCells(nCutCells_);
+ generateTriPoints
+ (
+ cVals,
+ pVals,
- generateTriPoints
- (
- cVals,
- pVals,
+ mesh_.cellCentres(),
+ mesh_.points(),
- mesh_.cellCentres(),
- mesh_.points(),
+ snappedPoints,
+ snappedCc,
+ snappedPoint,
- snappedPoints,
- snappedCc,
- snappedPoint,
+ triPoints,
+ triMeshCells
+ );
- triPoints,
- triMeshCells
- );
+ if (debug)
+ {
+ Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
+ << " unmerged triangles." << endl;
+ }
- if (debug)
- {
- Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
- << " unmerged triangles." << endl;
- }
- // Merge points and compact out non-valid triangles
- labelList triMap; // merged to unmerged triangle
- triSurface::operator=
- (
- stitchTriPoints
+ label nOldPoints = triPoints.size();
+
+ // Trimmed to original triangle
+ DynamicList<label> trimTriMap;
+ // Trimmed to original point
+ labelList trimTriPointMap;
+ if (bounds_ != boundBox::greatBox)
+ {
+ isoSurface::trimToBox
+ (
+ treeBoundBox(bounds_),
+ triPoints, // new points
+ trimTriMap, // map from (new) triangle to original
+ trimTriPointMap, // map from (new) point to original
+ interpolatedPoints_, // labels of newly introduced points
+ interpolatedOldPoints_, // and their interpolation
+ interpolationWeights_
+ );
+ triMeshCells = labelField(triMeshCells, trimTriMap);
+ }
+
+
+
+ // Merge points and compact out non-valid triangles
+ labelList triMap;
+ triSurface::operator=
(
- regularise, // check for duplicate tris
- triPoints,
- triPointMergeMap_, // unmerged to merged point
- triMap
- )
- );
+ stitchTriPoints
+ (
+ regularise, // check for duplicate tris
+ triPoints,
+ triPointMergeMap_, // unmerged to merged point
+ triMap // merged to unmerged triangle
+ )
+ );
- if (debug)
- {
- Pout<< "isoSurfaceCell : generated " << triMap.size()
- << " merged triangles." << endl;
- }
+ if (debug)
+ {
+ Pout<< "isoSurfaceCell : generated " << triMap.size()
+ << " merged triangles." << endl;
+ }
- meshCells_.setSize(triMap.size());
- forAll(triMap, i)
- {
- meshCells_[i] = triMeshCells[triMap[i]];
+ if (bounds_ != boundBox::greatBox)
+ {
+ // Adjust interpolatedPoints_
+ inplaceRenumber(triPointMergeMap_, interpolatedPoints_);
+
+ // Adjust triPointMergeMap_
+ labelList newTriPointMergeMap(nOldPoints, -1);
+ forAll(trimTriPointMap, trimPointI)
+ {
+ label oldPointI = trimTriPointMap[trimPointI];
+ if (oldPointI >= 0)
+ {
+ label pointI = triPointMergeMap_[trimPointI];
+ if (pointI >= 0)
+ {
+ newTriPointMergeMap[oldPointI] = pointI;
+ }
+ }
+ }
+ triPointMergeMap_.transfer(newTriPointMergeMap);
+ }
+
+ meshCells_.setSize(triMap.size());
+ forAll(triMap, i)
+ {
+ meshCells_[i] = triMeshCells[triMap[i]];
+ }
}
+
if (debug)
{
Pout<< "isoSurfaceCell : checking " << size()
@@ -1730,8 +1645,6 @@ Foam::isoSurfaceCell::isoSurfaceCell
meshCells_ = labelField(meshCells_, subsetTriMap);
inplaceRenumber(reversePointMap, triPointMergeMap_);
}
-
- //orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
}
}
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
index 7e12efc7bbead6bdac2d7bc23561734530dfe3ca..910f3fe46ed7bdeef88d173fa1c64cceb705331c 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -48,6 +48,7 @@ SourceFiles
#include "labelPair.H"
#include "pointIndexHit.H"
#include "PackedBoolList.H"
+#include "boundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -55,6 +56,7 @@ namespace Foam
{
class polyMesh;
+class plane;
/*---------------------------------------------------------------------------*\
Class isoSurfaceCell Declaration
@@ -91,6 +93,9 @@ class isoSurfaceCell
//- isoSurfaceCell value
const scalar iso_;
+ //- Optional bounds
+ const boundBox bounds_;
+
//- When to merge points
const scalar mergeDistance_;
@@ -106,6 +111,15 @@ class isoSurfaceCell
//- For every unmerged triangle point the point in the triSurface
labelList triPointMergeMap_;
+ //- triSurface points that have weighted interpolation
+ DynamicList<label> interpolatedPoints_;
+
+ //- corresponding original, unmerged points
+ DynamicList<FixedList<label, 3> > interpolatedOldPoints_;
+
+ //- corresponding weights
+ DynamicList<FixedList<scalar, 3> > interpolationWeights_;
+
// Private Member Functions
@@ -214,19 +228,15 @@ class isoSurfaceCell
void generateTriPoints
(
const DynamicList<Type>& snapped,
- const scalar isoVal0,
const scalar s0,
const Type& p0,
const label p0Index,
- const scalar isoVal1,
const scalar s1,
const Type& p1,
const label p1Index,
- const scalar isoVal2,
const scalar s2,
const Type& p2,
const label p2Index,
- const scalar isoVal3,
const scalar s3,
const Type& p3,
const label p3Index,
@@ -271,27 +281,6 @@ class isoSurfaceCell
Map<labelList>& edgeFacesRest
) const;
- ////- Determine orientation
- //static void walkOrientation
- //(
- // const triSurface& surf,
- // const List<FixedList<label, 3> >& faceEdges,
- // const labelList& edgeFace0,
- // const labelList& edgeFace1,
- // const label seedTriI,
- // labelList& flipState
- //);
-
- ////- Orient surface
- //static void orientSurface
- //(
- // triSurface&,
- // const List<FixedList<label, 3> >& faceEdges,
- // const labelList& edgeFace0,
- // const labelList& edgeFace1,
- // const Map<labelList>& edgeFacesRest
- //);
-
//- Is triangle (given by 3 edges) not fully connected?
static bool danglingTriangle
(
@@ -317,8 +306,6 @@ class isoSurfaceCell
labelList& newToOldPoints
);
- //- Combine all triangles inside a cell into a minimal triangulation
- void combineCellTriangles();
public:
@@ -336,6 +323,7 @@ public:
const scalarField& pointValues,
const scalar iso,
const bool regularise,
+ const boundBox& bounds = boundBox::greatBox,
const scalar mergeTol = 1e-6 // fraction of bounding box
);
@@ -348,24 +336,6 @@ public:
return meshCells_;
}
- //- For every unmerged triangle point the point in the triSurface
- const labelList triPointMergeMap() const
- {
- return triPointMergeMap_;
- }
-
-
- //- Interpolates cCoords,pCoords. Takes the original fields
- // used to create the iso surface.
- template<class Type>
- tmp<Field<Type> > interpolate
- (
- const scalarField& cVals,
- const scalarField& pVals,
- const Field<Type>& cCoords,
- const Field<Type>& pCoords
- ) const;
-
//- Interpolates cCoords,pCoords.
template<class Type>
tmp<Field<Type> > interpolate
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
index c78aafeee50525ed763dab6863083546a160084e..8fdee4831b82a2b908addba634e66485077267fe 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
@@ -26,6 +26,7 @@ License
#include "isoSurfaceCell.H"
#include "polyMesh.H"
#include "tetMatcher.H"
+#include "isoSurface.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@@ -76,22 +77,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
const DynamicList<Type>& snapped,
- const scalar isoVal0,
const scalar s0,
const Type& p0,
const label p0Index,
- const scalar isoVal1,
const scalar s1,
const Type& p1,
const label p1Index,
- const scalar isoVal2,
const scalar s2,
const Type& p2,
const label p2Index,
- const scalar isoVal3,
const scalar s3,
const Type& p3,
const label p3Index,
@@ -124,160 +121,196 @@ void Foam::isoSurfaceCell::generateTriPoints
case 0x0F:
break;
- case 0x0E:
case 0x01:
+ case 0x0E:
{
- // 0 is common point. Orient such that normal points in positive
- // gradient direction
- if (isoVal0 >= isoVal1)
- {
- pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
- pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- }
- else
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+ );
+ if (triIndex == 0x0E)
{
- pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
- pts.append(generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index));
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x0D:
case 0x02:
+ case 0x0D:
{
- // 1 is common point
- if (isoVal1 >= isoVal0)
- {
- pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- }
- else
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+ );
+
+ if (triIndex == 0x0D)
{
- pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x0C:
case 0x03:
+ case 0x0C:
{
- Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
- Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
-
- if (isoVal0 >= isoVal3)
- {
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- pts.append(s02);
- pts.append(s13);
- pts.append(s13);
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- pts.append(s02);
- }
- else
+ Type p0p2 =
+ generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
+ Type p1p3 =
+ generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
+
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+ );
+ pts.append(p1p3);
+ pts.append(p0p2);
+
+ pts.append(p1p3);
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+ );
+ pts.append(p0p2);
+
+ if (triIndex == 0x0C)
{
- pts.append(s02);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- pts.append(s13);
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- pts.append(s13);
- pts.append(s02);
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-5], pts[sz-4]);
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x0B:
case 0x04:
+ case 0x0B:
{
- // 2 is common point
- if (isoVal2 >= isoVal0)
+ pts.append
+ (
+ generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)
+ );
+
+ if (triIndex == 0x0B)
{
- pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
- pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
- }
- else
- {
- pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index));
- pts.append(generatePoint(snapped,s2,p2,p2Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index));
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x0A:
case 0x05:
+ case 0x0A:
{
- Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
- Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
- if (isoVal3 >= isoVal0)
+ Type p0p1 =
+ generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+ Type p2p3 =
+ generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+ pts.append(p0p1);
+ pts.append(p2p3);
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)
+ );
+
+ pts.append(p0p1);
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)
+ );
+ pts.append(p2p3);
+
+ if (triIndex == 0x0A)
{
- pts.append(s01);
- pts.append(s23);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- pts.append(s01);
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- pts.append(s23);
- }
- else
- {
- pts.append(s23);
- pts.append(s01);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index));
- pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
- pts.append(s01);
- pts.append(s23);
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-5], pts[sz-4]);
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x09:
case 0x06:
+ case 0x09:
{
- Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
- Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
-
- if (isoVal3 >= isoVal1)
+ Type p0p1 =
+ generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
+ Type p2p3 =
+ generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
+
+ pts.append(p0p1);
+ pts.append
+ (
+ generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)
+ );
+ pts.append(p2p3);
+
+ pts.append(p0p1);
+ pts.append(p2p3);
+ pts.append
+ (
+ generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)
+ );
+
+ if (triIndex == 0x09)
{
- pts.append(s01);
- pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
- pts.append(s23);
- pts.append(s01);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
- pts.append(s23);
- }
- else
- {
- pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index));
- pts.append(s01);
- pts.append(s23);
- pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index));
- pts.append(s01);
- pts.append(s23);
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-5], pts[sz-4]);
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
- case 0x07:
case 0x08:
+ case 0x07:
{
- // 3 is common point
- if (isoVal3 >= isoVal0)
+ pts.append
+ (
+ generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)
+ );
+ pts.append
+ (
+ generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)
+ );
+ if (triIndex == 0x07)
{
- pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
- pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
- }
- else
- {
- pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index));
- pts.append(generatePoint(snapped,s3,p3,p3Index,s0,p0,p0Index));
- pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index));
+ // Flip normals
+ label sz = pts.size();
+ Swap(pts[sz-2], pts[sz-1]);
}
}
break;
@@ -341,22 +374,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
snappedPoints,
- pVals_[f0[1]],
pVals[f0[1]],
pCoords[f0[1]],
snappedPoint[f0[1]],
- pVals_[f0[0]],
pVals[f0[0]],
pCoords[f0[0]],
snappedPoint[f0[0]],
- pVals_[f0[2]],
pVals[f0[2]],
pCoords[f0[2]],
snappedPoint[f0[2]],
- pVals_[oppositeI],
pVals[oppositeI],
pCoords[oppositeI],
snappedPoint[oppositeI],
@@ -370,22 +399,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
snappedPoints,
- pVals_[f0[0]],
pVals[f0[0]],
pCoords[f0[0]],
snappedPoint[f0[0]],
- pVals_[f0[1]],
pVals[f0[1]],
pCoords[f0[1]],
snappedPoint[f0[1]],
- pVals_[f0[2]],
pVals[f0[2]],
pCoords[f0[2]],
snappedPoint[f0[2]],
- pVals_[oppositeI],
pVals[oppositeI],
pCoords[oppositeI],
snappedPoint[oppositeI],
@@ -411,7 +436,6 @@ void Foam::isoSurfaceCell::generateTriPoints
}
label fp = f.fcIndex(fp0);
-
for (label i = 2; i < f.size(); i++)
{
label nextFp = f.fcIndex(fp);
@@ -425,22 +449,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
snappedPoints,
- pVals_[tri[1]],
pVals[tri[1]],
pCoords[tri[1]],
snappedPoint[tri[1]],
- pVals_[tri[0]],
pVals[tri[0]],
pCoords[tri[0]],
snappedPoint[tri[0]],
- pVals_[tri[2]],
pVals[tri[2]],
pCoords[tri[2]],
snappedPoint[tri[2]],
- cVals_[cellI],
cVals[cellI],
cCoords[cellI],
snappedCc[cellI],
@@ -454,22 +474,18 @@ void Foam::isoSurfaceCell::generateTriPoints
(
snappedPoints,
- pVals_[tri[0]],
pVals[tri[0]],
pCoords[tri[0]],
snappedPoint[tri[0]],
- pVals_[tri[1]],
pVals[tri[1]],
pCoords[tri[1]],
snappedPoint[tri[1]],
- pVals_[tri[2]],
pVals[tri[2]],
pCoords[tri[2]],
snappedPoint[tri[2]],
- cVals_[cellI],
cVals[cellI],
cCoords[cellI],
snappedCc[cellI],
@@ -495,7 +511,7 @@ void Foam::isoSurfaceCell::generateTriPoints
if (countNotFoundTets > 0)
{
- WarningIn("Foam::isoSurfaceCell::generateTriPoints")
+ WarningIn("Foam::isoSurfaceCell::generateTriPoints(..)")
<< "Could not find " << countNotFoundTets
<< " tet base points, which may lead to inverted triangles."
<< endl;
@@ -510,13 +526,11 @@ template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::isoSurfaceCell::interpolate
(
- const scalarField& cVals,
- const scalarField& pVals,
const Field<Type>& cCoords,
const Field<Type>& pCoords
) const
{
- DynamicList<Type> triPoints(nCutCells_);
+ DynamicList<Type> triPoints(3*nCutCells_);
DynamicList<label> triMeshCells(nCutCells_);
// Dummy snap data
@@ -524,59 +538,6 @@ Foam::isoSurfaceCell::interpolate
labelList snappedCc(mesh_.nCells(), -1);
labelList snappedPoint(mesh_.nPoints(), -1);
-
- generateTriPoints
- (
- cVals,
- pVals,
-
- cCoords,
- pCoords,
-
- snappedPoints,
- snappedCc,
- snappedPoint,
-
- triPoints,
- triMeshCells
- );
-
-
- // One value per point
- tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
- Field<Type>& values = tvalues();
-
- forAll(triPoints, i)
- {
- label mergedPointI = triPointMergeMap_[i];
-
- if (mergedPointI >= 0)
- {
- values[mergedPointI] = triPoints[i];
- }
- }
-
- return tvalues;
-}
-
-
-template<class Type>
-Foam::tmp<Foam::Field<Type> >
-Foam::isoSurfaceCell::interpolate
-(
- const Field<Type>& cCoords,
- const Field<Type>& pCoords
-) const
-{
- DynamicList<Type> triPoints(nCutCells_);
- DynamicList<label> triMeshCells(nCutCells_);
-
- // Dummy snap data
- DynamicList<Type> snappedPoints;
- labelList snappedCc(mesh_.nCells(), -1);
- labelList snappedPoint(mesh_.nPoints(), -1);
-
-
generateTriPoints
(
cVals_,
@@ -593,22 +554,15 @@ Foam::isoSurfaceCell::interpolate
triMeshCells
);
-
- // One value per point
- tmp<Field<Type> > tvalues(new Field<Type>(points().size()));
- Field<Type>& values = tvalues();
-
- forAll(triPoints, i)
- {
- label mergedPointI = triPointMergeMap_[i];
-
- if (mergedPointI >= 0)
- {
- values[mergedPointI] = triPoints[i];
- }
- }
-
- return tvalues;
+ return isoSurface::interpolate
+ (
+ points().size(),
+ triPointMergeMap_,
+ interpolatedPoints_,
+ interpolatedOldPoints_,
+ interpolationWeights_,
+ triPoints
+ );
}
diff --git a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
index 6f64c9fe541198ea8524f097b69921702d525264..0c4f0f7811d53181d6b56251dc26a8213ac41416 100644
--- a/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
+++ b/src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -120,7 +120,7 @@ Foam::isoSurface::adaptPatchFields
{
fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
(
- fld.boundaryField()[patchI]
+ sliceFld.boundaryField()[patchI]
);
const scalarField& w = mesh.weights().boundaryField()[patchI];
@@ -189,6 +189,8 @@ Type Foam::isoSurface::generatePoint
}
+// Note: cannot use simpler isoSurfaceCell::generateTriPoints since
+// the need here to sometimes pass in remote 'snappedPoints'
template<class Type>
void Foam::isoSurface::generateTriPoints
(
@@ -240,8 +242,8 @@ void Foam::isoSurface::generateTriPoints
case 0x0F:
break;
- case 0x0E:
case 0x01:
+ case 0x0E:
points.append
(
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1)
@@ -254,10 +256,16 @@ void Foam::isoSurface::generateTriPoints
(
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
);
+ if (triIndex == 0x0E)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x0D:
case 0x02:
+ case 0x0D:
points.append
(
generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0)
@@ -270,97 +278,133 @@ void Foam::isoSurface::generateTriPoints
(
generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
);
+ if (triIndex == 0x0D)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x0C:
case 0x03:
- {
- Type tp1 =
- generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
- Type tp2 =
- generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
+ case 0x0C:
+ {
+ Type p0p2 =
+ generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
+ Type p1p3 =
+ generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
- points.append
- (
- generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
- );
- points.append(tp1);
- points.append(tp2);
- points.append(tp2);
- points.append
- (
- generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
- );
- points.append(tp1);
- }
+ points.append
+ (
+ generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+ );
+ points.append(p1p3);
+ points.append(p0p2);
+
+ points.append(p1p3);
+ points.append
+ (
+ generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+ );
+ points.append(p0p2);
+ }
+ if (triIndex == 0x0C)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-5], points[sz-4]);
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x0B:
case 0x04:
- {
- points.append
- (
- generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
- );
- points.append
- (
- generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
- );
- points.append
- (
- generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
- );
- }
+ case 0x0B:
+ {
+ points.append
+ (
+ generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
+ );
+ points.append
+ (
+ generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
+ );
+ points.append
+ (
+ generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3)
+ );
+ }
+ if (triIndex == 0x0B)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x0A:
case 0x05:
- {
- Type tp0 =
- generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
- Type tp1 =
- generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+ case 0x0A:
+ {
+ Type p0p1 =
+ generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+ Type p2p3 =
+ generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+
+ points.append(p0p1);
+ points.append(p2p3);
+ points.append
+ (
+ generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
+ );
- points.append(tp0);
- points.append(tp1);
- points.append
- (
- generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
- );
- points.append(tp0);
- points.append
- (
- generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
- );
- points.append(tp1);
- }
+ points.append(p0p1);
+ points.append
+ (
+ generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
+ );
+ points.append(p2p3);
+ }
+ if (triIndex == 0x0A)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-5], points[sz-4]);
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x09:
case 0x06:
- {
- Type tp0 =
- generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
- Type tp1 =
- generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
+ case 0x09:
+ {
+ Type p0p1 =
+ generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
+ Type p2p3 =
+ generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
- points.append(tp0);
- points.append
- (
- generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
- );
- points.append(tp1);
- points.append(tp0);
- points.append
- (
- generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
- );
- points.append(tp1);
- }
+ points.append(p0p1);
+ points.append
+ (
+ generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
+ );
+ points.append(p2p3);
+
+ points.append(p0p1);
+ points.append(p2p3);
+ points.append
+ (
+ generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
+ );
+ }
+ if (triIndex == 0x09)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-5], points[sz-4]);
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
- case 0x07:
case 0x08:
+ case 0x07:
points.append
(
generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0)
@@ -373,6 +417,12 @@ void Foam::isoSurface::generateTriPoints
(
generatePoint(s3,p3,hasSnap3,snapP3,s1,p1,hasSnap1,snapP1)
);
+ if (triIndex == 0x07)
+ {
+ // Flip normals
+ label sz = points.size();
+ Swap(points[sz-2], points[sz-1]);
+ }
break;
}
}
@@ -689,6 +739,69 @@ void Foam::isoSurface::generateTriPoints
//}
+template<class Type>
+Foam::tmp<Foam::Field<Type> >
+Foam::isoSurface::interpolate
+(
+ const label nPoints,
+ const labelList& triPointMergeMap,
+ const labelList& interpolatedPoints,
+ const List<FixedList<label, 3> >& interpolatedOldPoints,
+ const List<FixedList<scalar, 3> >& interpolationWeights,
+ const DynamicList<Type>& unmergedValues
+)
+{
+ // One value per point
+ tmp<Field<Type> > tvalues(new Field<Type>(nPoints, pTraits<Type>::zero));
+ Field<Type>& values = tvalues();
+
+
+ // Pass1: unweighted average of merged point values
+ {
+ labelList nValues(values.size(), 0);
+
+ forAll(unmergedValues, i)
+ {
+ label mergedPointI = triPointMergeMap[i];
+
+ if (mergedPointI >= 0)
+ {
+ values[mergedPointI] += unmergedValues[i];
+ nValues[mergedPointI]++;
+ }
+ }
+
+ forAll(values, i)
+ {
+ if (nValues[i] > 0)
+ {
+ values[i] /= scalar(nValues[i]);
+ }
+ }
+ }
+
+
+ // Pass2: weighted average for remaining values (from clipped triangles)
+
+ forAll(interpolatedPoints, i)
+ {
+ label pointI = interpolatedPoints[i];
+ const FixedList<label, 3>& oldPoints = interpolatedOldPoints[i];
+ const FixedList<scalar, 3>& w = interpolationWeights[i];
+
+ // Note: zeroing should not be necessary if interpolation only done
+ // for newly introduced points (i.e. not in triPointMergeMap)
+ values[pointI] = pTraits<Type>::zero;
+ forAll(oldPoints, j)
+ {
+ values[pointI] = w[j]*unmergedValues[oldPoints[j]];
+ }
+ }
+
+ return tvalues;
+}
+
+
template<class Type>
Foam::tmp<Foam::Field<Type> >
Foam::isoSurface::interpolate
@@ -707,7 +820,7 @@ Foam::isoSurface::interpolate
> > c2(adaptPatchFields(cCoords));
- DynamicList<Type> triPoints(nCutCells_);
+ DynamicList<Type> triPoints(3*nCutCells_);
DynamicList<label> triMeshCells(nCutCells_);
// Dummy snap data
@@ -731,52 +844,15 @@ Foam::isoSurface::interpolate
triMeshCells
);
-
- // One value per point
- tmp<Field<Type> > tvalues
+ return interpolate
(
- new Field<Type>(points().size(), pTraits<Type>::zero)
+ points().size(),
+ triPointMergeMap_,
+ interpolatedPoints_,
+ interpolatedOldPoints_,
+ interpolationWeights_,
+ triPoints
);
- Field<Type>& values = tvalues();
- labelList nValues(values.size(), 0);
-
- forAll(triPoints, i)
- {
- label mergedPointI = triPointMergeMap_[i];
-
- if (mergedPointI >= 0)
- {
- values[mergedPointI] += triPoints[i];
- nValues[mergedPointI]++;
- }
- }
-
- if (debug)
- {
- Pout<< "nValues:" << values.size() << endl;
- label nMult = 0;
- forAll(nValues, i)
- {
- if (nValues[i] == 0)
- {
- FatalErrorIn("isoSurface::interpolate(..)")
- << "point:" << i << " nValues:" << nValues[i]
- << abort(FatalError);
- }
- else if (nValues[i] > 1)
- {
- nMult++;
- }
- }
- Pout<< "Of which mult:" << nMult << endl;
- }
-
- forAll(values, i)
- {
- values[i] /= scalar(nValues[i]);
- }
-
- return tvalues;
}
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
index bb171c2c0e98a814e752dad4fe5865f2118b8f47..0399eeae06e824ce8b4dd8fd7ca0db77c6668fdd 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
@@ -122,18 +122,52 @@ void Foam::sampledIsoSurface::getIsoFields() const
// Get pointField
// ~~~~~~~~~~~~~~
+ // In case of multiple iso values we don't want to calculate multiple e.g.
+ // "volPointInterpolate(p)" so register it and re-use it. This is the
+ // same as the 'cache' functionality from volPointInterpolate but
+ // unfortunately that one does not guarantee that the field pointer
+ // remain: e.g. some other
+ // functionObject might delete the cached version.
+ // (volPointInterpolation::interpolate with cache=false deletes any
+ // registered one or if mesh.changing())
+
if (!subMeshPtr_.valid())
{
- word pointFldName = "volPointInterpolate(" + isoField_ + ')';
+ const word pointFldName =
+ "volPointInterpolate_"
+ + type()
+ + "("
+ + isoField_
+ + ')';
if (fvm.foundObject<pointScalarField>(pointFldName))
{
if (debug)
{
- Info<< "sampledIsoSurface::getIsoFields() : lookup pointField "
- << pointFldName << endl;
+ Info<< "sampledIsoSurface::getIsoFields() :"
+ << " lookup pointField " << pointFldName << endl;
}
- pointFieldPtr_ = &fvm.lookupObject<pointScalarField>(pointFldName);
+ const pointScalarField& pfld = fvm.lookupObject<pointScalarField>
+ (
+ pointFldName
+ );
+
+ if (!pfld.upToDate(*volFieldPtr_))
+ {
+ if (debug)
+ {
+ Info<< "sampledIsoSurface::getIsoFields() :"
+ << " updating pointField " << pointFldName << endl;
+ }
+ // Update the interpolated value
+ volPointInterpolation::New(fvm).interpolate
+ (
+ *volFieldPtr_,
+ const_cast<pointScalarField&>(pfld)
+ );
+ }
+
+ pointFieldPtr_ = &pfld;
}
else
{
@@ -141,33 +175,23 @@ void Foam::sampledIsoSurface::getIsoFields() const
if (debug)
{
- Info<< "sampledIsoSurface::getIsoFields() : "
- << "checking pointField " << pointFldName
- << " for same time " << fvm.time().timeName()
- << endl;
+ Info<< "sampledIsoSurface::getIsoFields() :"
+ << " creating and storing pointField "
+ << pointFldName << " for time "
+ << fvm.time().timeName() << endl;
}
- if
+ tmp<pointScalarField> tpfld
(
- storedPointFieldPtr_.empty()
- || (fvm.time().timeName() != storedPointFieldPtr_().instance())
- )
- {
- if (debug)
- {
- Info<< "sampledIsoSurface::getIsoFields() :"
- << " interpolating volField " << volFieldPtr_->name()
- << " to get pointField " << pointFldName << endl;
- }
-
- storedPointFieldPtr_.reset
+ volPointInterpolation::New(fvm).interpolate
(
- volPointInterpolation::New(fvm)
- .interpolate(*volFieldPtr_).ptr()
- );
- storedPointFieldPtr_->checkOut();
- pointFieldPtr_ = storedPointFieldPtr_.operator->();
- }
+ *volFieldPtr_,
+ pointFldName,
+ false
+ )
+ );
+ pointFieldPtr_ = tpfld.ptr();
+ const_cast<pointScalarField*>(pointFieldPtr_)->store();
}
@@ -233,8 +257,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
// Pointfield on submesh
- word pointFldName =
- "volPointInterpolate("
+ const word pointFldName =
+ "volPointInterpolate_"
+ + type()
+ + "("
+ volSubFieldPtr_->name()
+ ')';
@@ -245,11 +271,28 @@ void Foam::sampledIsoSurface::getIsoFields() const
Info<< "sampledIsoSurface::getIsoFields() :"
<< " submesh lookup pointField " << pointFldName << endl;
}
- storedPointSubFieldPtr_.clear();
- pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
+ const pointScalarField& pfld = subFvm.lookupObject<pointScalarField>
(
pointFldName
);
+
+ if (!pfld.upToDate(*volSubFieldPtr_))
+ {
+ if (debug)
+ {
+ Info<< "sampledIsoSurface::getIsoFields() :"
+ << " updating submesh pointField "
+ << pointFldName << endl;
+ }
+ // Update the interpolated value
+ volPointInterpolation::New(subFvm).interpolate
+ (
+ *volSubFieldPtr_,
+ const_cast<pointScalarField&>(pfld)
+ );
+ }
+
+ pointFieldPtr_ = &pfld;
}
else
{
@@ -260,15 +303,15 @@ void Foam::sampledIsoSurface::getIsoFields() const
<< volSubFieldPtr_->name()
<< " to get submesh pointField " << pointFldName << endl;
}
- storedPointSubFieldPtr_.reset
+ tmp<pointScalarField> tpfld
(
volPointInterpolation::New
(
subFvm
- ).interpolate(*volSubFieldPtr_).ptr()
+ ).interpolate(*volSubFieldPtr_)
);
- storedPointSubFieldPtr_->checkOut();
- pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->();
+ pointSubFieldPtr_ = tpfld.ptr();
+ const_cast<pointScalarField*>(pointSubFieldPtr_)->store();
}
@@ -349,28 +392,34 @@ bool Foam::sampledIsoSurface::updateGeometry() const
if (subMeshPtr_.valid())
{
+ const volScalarField& vfld = *volSubFieldPtr_;
+
surfPtr_.reset
(
new isoSurface
(
- *volSubFieldPtr_,
+ vfld,
*pointSubFieldPtr_,
isoVal_,
regularise_,
+ bounds_,
mergeTol_
)
);
}
else
{
+ const volScalarField& vfld = *volFieldPtr_;
+
surfPtr_.reset
(
new isoSurface
(
- *volFieldPtr_,
+ vfld,
*pointFieldPtr_,
isoVal_,
regularise_,
+ bounds_,
mergeTol_
)
);
@@ -412,6 +461,7 @@ Foam::sampledIsoSurface::sampledIsoSurface
sampledSurface(name, mesh, dict),
isoField_(dict.lookup("isoField")),
isoVal_(readScalar(dict.lookup("isoValue"))),
+ bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", false)),
@@ -422,7 +472,6 @@ Foam::sampledIsoSurface::sampledIsoSurface
prevTimeIndex_(-1),
storedVolFieldPtr_(NULL),
volFieldPtr_(NULL),
- storedPointFieldPtr_(NULL),
pointFieldPtr_(NULL)
{
if (!sampledSurface::interpolate())
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
index 8cfd14c3c15cd3b0d26d2e21510d42acc89c96ec..488ba56452544309557a54de64611341025913d8 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
@@ -63,6 +63,9 @@ class sampledIsoSurface
//- Iso value
const scalar isoVal_;
+ //- Optional bounding box to trim triangles against
+ const boundBox bounds_;
+
//- Merge tolerance
const scalar mergeTol_;
@@ -94,7 +97,6 @@ class sampledIsoSurface
mutable const volScalarField* volFieldPtr_;
//- Cached pointfield
- mutable autoPtr<pointScalarField> storedPointFieldPtr_;
mutable const pointScalarField* pointFieldPtr_;
// And on subsetted mesh
@@ -107,7 +109,6 @@ class sampledIsoSurface
mutable const volScalarField* volSubFieldPtr_;
//- Cached pointfield
- mutable autoPtr<pointScalarField> storedPointSubFieldPtr_;
mutable const pointScalarField* pointSubFieldPtr_;
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
index 239217a352c365316da312e91a14992ea8a2fd0c..4c21257ce7b07de49a65b131d687531b03e2d74d 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -148,7 +148,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
cellAvg,
pointFld().internalField(),
isoVal_,
- regularise_
+ regularise_,
+ bounds_
);
const_cast<sampledIsoSurfaceCell&>
@@ -166,7 +167,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
cellFld.internalField(),
pointFld().internalField(),
isoVal_,
- regularise_
+ regularise_,
+ bounds_
);
const_cast<sampledIsoSurfaceCell&>
@@ -185,6 +187,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
<< " average : " << average_ << nl
<< " isoField : " << isoField_ << nl
<< " isoValue : " << isoVal_ << nl
+ << " bounds : " << bounds_ << nl
<< " points : " << points().size() << nl
<< " tris : " << triSurface::size() << nl
<< " cut cells : " << meshCells_.size() << endl;
@@ -206,6 +209,7 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
sampledSurface(name, mesh, dict),
isoField_(dict.lookup("isoField")),
isoVal_(readScalar(dict.lookup("isoValue"))),
+ bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", true)),
zoneKey_(keyType::null),
diff --git a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
index 7021bf0b97e35b941ab743cb83ea1b985a96fb86..0cfa6a7461a9ec02b6b022e989d24c961222f528 100644
--- a/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
+++ b/src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
@@ -62,6 +62,9 @@ class sampledIsoSurfaceCell
//- Iso value
const scalar isoVal_;
+ //- Optional bounding box to trim triangles against
+ const boundBox bounds_;
+
//- Whether to coarse
const Switch regularise_;
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
index 7da10d7a57db7ce97e6e7454106c6f444cbe696f..7030b20d5e62a97e4fe824b8ac2ebcc41b2230a2 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -230,6 +230,7 @@ void Foam::sampledCuttingPlane::createGeometry()
pointDistance_,
0.0,
regularise_,
+ bounds_,
mergeTol_
)
//new isoSurfaceCell
@@ -262,6 +263,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
:
sampledSurface(name, mesh, dict),
plane_(dict),
+ bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", false)),
diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
index a8aa477407fe466a140d6f1b695a6ccb4f0e6a33..308bdcd46a01aeb8bb5fa5aad007141ade2fd886 100644
--- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
+++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
@@ -60,6 +60,9 @@ class sampledCuttingPlane
//- Plane
const plane plane_;
+ //- Optional bounding box to trim triangles against
+ const boundBox bounds_;
+
//- Merge tolerance
const scalar mergeTol_;
diff --git a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
index d07a6391f120f7a2a2154d547de6bb01be7f63a8..faeccd69a3f9669228860950d337026f0a25ba4f 100644
--- a/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
+++ b/src/sampling/sampledSurface/sampledPatch/sampledPatch.C
@@ -2,7 +2,7 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation |
-------------------------------------------------------------------------------
License
@@ -84,11 +84,7 @@ const Foam::labelList& Foam::sampledPatch::patchIDs() const
{
if (patchIDs_.empty())
{
- patchIDs_ = mesh().boundaryMesh().patchSet
- (
- patchNames_,
- false
- ).sortedToc();
+ patchIDs_ = mesh().boundaryMesh().patchSet(patchNames_).sortedToc();
}
return patchIDs_;
}
diff --git a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
index 47a9f9ebbc4d0e0e67bfb7b63e98c683c6eb56de..4e12874c2545f21c1f31f5a2aa4f8d7a4519dfbd 100644
--- a/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
+++ b/src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
@@ -2,8 +2,8 @@
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
- \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
- \\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd
+ \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
+ \\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of OpenFOAM.
@@ -162,19 +162,7 @@ void Foam::sampledSurfaces::write()
const label nFields = classifyFields();
- if (Pstream::master())
- {
- if (debug)
- {
- Pout<< "Creating directory "
- << outputPath_/obr_.time().timeName() << nl << endl;
-
- }
-
- mkDir(outputPath_/obr_.time().timeName());
- }
-
- // Write geometry first if required,
+ // write geometry first if required,
// or when no fields would otherwise be written
if (nFields == 0 || formatter_->separateGeometry())
{
diff --git a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
index 094ae45cfcb488323910c412de721701d661b8af..654b3d438f23d3a100eee1d8499b6bfd3cc19e45 100644
--- a/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
+++ b/src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
@@ -644,8 +644,25 @@ bool Foam::sampledTriSurfaceMesh::update()
return false;
}
+ // Calculate surface and mesh overlap bounding box
+ treeBoundBox bb
+ (
+ surface_.triSurface::points(),
+ surface_.triSurface::meshPoints()
+ );
+ bb.min() = max(bb.min(), mesh().bounds().min());
+ bb.max() = min(bb.max(), mesh().bounds().max());
+
+ // Extend a bit
+ const vector span(bb.span());
+
+ bb.min() -= 0.5*span;
+ bb.max() += 0.5*span;
+
+ bb.inflate(1e-6);
+
// Mesh search engine, no triangulation of faces.
- meshSearch meshSearcher(mesh(), polyMesh::FACE_PLANES);
+ meshSearch meshSearcher(mesh(), bb, polyMesh::FACE_PLANES);
return update(meshSearcher);
}