zhyang-dev/openfoam
1
0
Fork 0
mirror of https://develop.openfoam.com/Development/openfoam.git synced 2025-11-28 03:28:01 +00:00
Code Issues Actions Packages Projects Releases Wiki Activity

ENH: sampledSurfaces: added 'bounds' option

Browse Source 
- bounds option (see $FOAM_UTILITIES/postProcessing/sampling/sample/sampleDict)
- fixes memory error http://www.openfoam.org/mantisbt/view.php?id=1487
- cleans up iso surface normal orientation
This commit is contained in:
mattijs
2015-11-26 10:41:34 +00:00
parent 907c9390e3
commit cee6887d68
24 changed files with 2249 additions and 1370 deletions
Download Patch File Download Diff File Expand all files Collapse all files

85
applications/utilities/postProcessing/sampling/sample/sampleDict
View File

@ -31,18 +31,66 @@ setFormat raw;
// dx : DX scalar or vector format // dx : DX scalar or vector format
// vtk : VTK ascii format // vtk : VTK ascii format
// raw : x y z value format for use with e.g. gnuplot 'splot'. // 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: // Note:
// other formats such as obj, stl, etc can also be written (by proxy) // - other formats such as obj, stl, etc can also be written (by proxy)
// but without any values! // 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; surfaceFormat vtk;
// optionally define extra controls for the output formats // Optionally define extra controls for the output formats
formatOptions formatOptions
{ {
ensight ensight
{ {
// ascii/binary format
format ascii; 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; type patchSeed;
axis xyz; axis xyz;
patches (".*Wall.*"); 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; 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; isoValue 0.5;
interpolate true; interpolate true;
//zone ABC; // Optional: zone only // zone ABC; // Optional: zone only
//exposedPatchName fixedWalls; // Optional: zone only // exposedPatchName fixedWalls; // Optional: zone only
// bounds (1 1 1)(2 2 2); // Optional: limit extent
// regularise false; // Optional: do not simplify // regularise false; // Optional: do not simplify
// mergeTol 1e-10; // Optional: fraction of mesh bounding box // mergeTol 1e-10; // Optional: fraction of mesh bounding box
@ -265,6 +325,9 @@ surfaces
isoValue 0.5; isoValue 0.5;
interpolate false; interpolate false;
regularise false; // do not simplify regularise false; // do not simplify
// bounds (1 1 1)(2 2 2); // Optional: limit extent
// mergeTol 1e-10; // Optional: fraction of mesh bounding box // mergeTol 1e-10; // Optional: fraction of mesh bounding box
// to merge points (default=1e-6) // to merge points (default=1e-6)
} }
@ -281,8 +344,10 @@ surfaces
} }
interpolate true; interpolate true;
//zone ABC; // Optional: zone only // zone ABC; // Optional: zone only
//exposedPatchName fixedWalls; // Optional: zone only // exposedPatchName fixedWalls; // Optional: zone only
// bounds (1 1 1)(2 2 2); // Optional: limit extent
// regularise false; // Optional: do not simplify // regularise false; // Optional: do not simplify
// mergeTol 1e-10; // Optional: fraction of mesh bounding box // mergeTol 1e-10; // Optional: fraction of mesh bounding box
@ -305,6 +370,8 @@ surfaces
// of isoSurfaceCell // of isoSurfaceCell
interpolate false; interpolate false;
regularise false; // Optional: do not simplify 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 // mergeTol 1e-10; // Optional: fraction of mesh bounding box
// to merge points (default=1e-6) // to merge points (default=1e-6)
} }

242
src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C Normal file
View File

@ -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]);
}
}
// ************************************************************************* //

122
src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -52,6 +52,8 @@ namespace Foam
class Istream; class Istream;
class Ostream; class Ostream;
class triPoints;
class plane;
// Forward declaration of friend functions and operators // Forward declaration of friend functions and operators
@ -71,6 +73,8 @@ inline Ostream& operator<<
const triangle<Point, PointRef>& const triangle<Point, PointRef>&
); );
typedef triangle<point, const point&> triPointRef;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class triangle Declaration Class triangle Declaration
@ -79,27 +83,98 @@ inline Ostream& operator<<
template<class Point, class PointRef> template<class Point, class PointRef>
class triangle 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 // Private data
PointRef a_, b_, c_; PointRef a_, b_, c_;
// Private Member Functions
//- 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: public:
//- Return types for classify
enum proxType
{
NONE,
POINT, // Close to point
EDGE // Close to edge
};
// Constructors // Constructors
//- Construct from three points //- Construct from three points
inline triangle(const Point& a, const Point& b, const Point& c); 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 //- Construct from three points in the list of points
// The indices could be from triFace etc. // The indices could be from triFace etc.
inline triangle inline triangle
@ -237,6 +312,27 @@ public:
pointHit& edgePoint pointHit& edgePoint
) const; ) 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 // IOstream operators
@ -264,6 +360,12 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "triangle.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif #endif
// ************************************************************************* // // ************************************************************************* //

75
src/OpenFOAM/meshes/primitiveShapes/triangle/triangleI.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -25,6 +25,7 @@ License
#include "IOstreams.H" #include "IOstreams.H"
#include "pointHit.H" #include "pointHit.H"
#include "triPoints.H"
#include "mathematicalConstants.H" #include "mathematicalConstants.H"
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * 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> template<class Point, class PointRef>
inline Foam::triangle<Point, PointRef>::triangle 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 * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * Ostream Operator * * * * * * * * * * * * * //
template<class Point, class PointRef> template<class Point, class PointRef>

4
src/mesh/extrudeModel/planeExtrusion/planeExtrusion.H
View File

@ -33,8 +33,8 @@ SeeAlso
\*---------------------------------------------------------------------------*/ \*---------------------------------------------------------------------------*/
#ifndef plane_H #ifndef planeExtrusion_H
#define plane_H #define planeExtrusion_H
#include "linearNormal.H" #include "linearNormal.H"

158
src/sampling/sampledSet/patchSeed/patchSeedSet.C
View File

@ -30,9 +30,8 @@ License
#include "treeDataFace.H" #include "treeDataFace.H"
#include "Time.H" #include "Time.H"
#include "meshTools.H" #include "meshTools.H"
//#include "Random.H"
// For 'facePoint' helper function only
#include "mappedPatchBase.H" #include "mappedPatchBase.H"
#include "indirectPrimitivePatch.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * 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 // Shuffle and truncate if in random mode
label totalSize = returnReduce(patchFaces.size(), sumOp<label>());
if (maxPoints_ < totalSize) if (maxPoints_ < totalSize)
{ {
// Check what fraction of maxPoints_ I need to generate locally. // 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)); labelList subset = identity(patchFaces.size());
Random& rndGen = rndGenPtr_();
labelList subset = identity(sz);
for (label iter = 0; iter < 4; iter++) for (label iter = 0; iter < 4; iter++)
{ {
forAll(subset, i) forAll(subset, i)
@ -115,7 +240,7 @@ void Foam::patchSeedSet::calcSamples
if (debug) if (debug)
{ {
Pout<< "In random mode : selected " << patchFaces.size() 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) forAll(patchFaces, i)
{ {
label faceI = 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 pointIndexHit info = mappedPatchBase::facePoint
( (
mesh(), mesh(),
@ -217,9 +345,15 @@ Foam::patchSeedSet::patchSeedSet
wordReList(dict.lookup("patches")) wordReList(dict.lookup("patches"))
) )
), ),
//searchDist_(readScalar(dict.lookup("maxDistance"))), maxPoints_(readLabel(dict.lookup("maxPoints"))),
//offsetDist_(readScalar(dict.lookup("offsetDist"))), selectedLocations_
maxPoints_(readLabel(dict.lookup("maxPoints"))) (
dict.lookupOrDefault<pointField>
(
"points",
pointField(0)
)
)
{ {
genSamples(); genSamples();

13
src/sampling/sampledSet/patchSeed/patchSeedSet.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2012-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -58,18 +58,15 @@ class patchSeedSet
//- Patches to sample //- Patches to sample
const labelHashSet patchSet_; const labelHashSet patchSet_;
// //- Maximum distance to look for nearest //- Maximum number of patch faces to seed (if in random subset mode)
// const scalar searchDist_;
//
// //- Offset distance
// const scalar offsetDist_;
//- Maximum number of patch faces to seed
const label maxPoints_; const label maxPoints_;
//- Random number generator (if maxPoints < num patch faces) //- Random number generator (if maxPoints < num patch faces)
autoPtr<Random> rndGenPtr_; autoPtr<Random> rndGenPtr_;
//- Patch faces to seed selected based on nearness to supplied points
const pointField selectedLocations_;
// Private Member Functions // Private Member Functions

13
src/sampling/sampledSurface/distanceSurface/distanceSurface.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -278,7 +278,8 @@ void Foam::distanceSurface::createGeometry()
cellDistance, cellDistance,
pointDistance_, pointDistance_,
distance_, distance_,
regularise_ regularise_,
bounds_
) )
); );
} }
@ -291,7 +292,8 @@ void Foam::distanceSurface::createGeometry()
cellDistance, cellDistance,
pointDistance_, pointDistance_,
distance_, distance_,
regularise_ regularise_,
bounds_
) )
); );
} }
@ -336,6 +338,7 @@ Foam::distanceSurface::distanceSurface
cell_(dict.lookupOrDefault("cell", true)), cell_(dict.lookupOrDefault("cell", true)),
regularise_(dict.lookupOrDefault("regularise", true)), regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", false)), average_(dict.lookupOrDefault("average", false)),
bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
zoneKey_(keyType::null), zoneKey_(keyType::null),
needsUpdate_(true), needsUpdate_(true),
isoSurfCellPtr_(NULL), isoSurfCellPtr_(NULL),
@ -364,7 +367,8 @@ Foam::distanceSurface::distanceSurface
const bool signedDistance, const bool signedDistance,
const bool cell, const bool cell,
const Switch regularise, const Switch regularise,
const Switch average const Switch average,
const boundBox& bounds
) )
: :
sampledSurface(name, mesh, interpolate), sampledSurface(name, mesh, interpolate),
@ -390,6 +394,7 @@ Foam::distanceSurface::distanceSurface
cell_(cell), cell_(cell),
regularise_(regularise), regularise_(regularise),
average_(average), average_(average),
bounds_(bounds),
zoneKey_(keyType::null), zoneKey_(keyType::null),
needsUpdate_(true), needsUpdate_(true),
isoSurfCellPtr_(NULL), isoSurfCellPtr_(NULL),

6
src/sampling/sampledSurface/distanceSurface/distanceSurface.H
View File

@ -75,6 +75,9 @@ class distanceSurface
//- Whether to recalculate cell values as average of point values //- Whether to recalculate cell values as average of point values
const Switch average_; 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 //- If restricted to zones, name of this zone or a regular expression
keyType zoneKey_; keyType zoneKey_;
@ -144,7 +147,8 @@ public:
const bool signedDistance, const bool signedDistance,
const bool cell, const bool cell,
const Switch regularise, const Switch regularise,
const Switch average const Switch average,
const boundBox& bounds = boundBox::greatBox
); );

1561
src/sampling/sampledSurface/isoSurface/isoSurface.C
View File

File diff suppressed because it is too large Load Diff

130
src/sampling/sampledSurface/isoSurface/isoSurface.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -76,6 +76,8 @@ namespace Foam
{ {
class fvMesh; class fvMesh;
class plane;
class treeBoundBox;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class isoSurface Declaration Class isoSurface Declaration
@ -116,10 +118,12 @@ class isoSurface
//- Regularise? //- Regularise?
const Switch regularise_; const Switch regularise_;
//- Optional bounds
const boundBox bounds_;
//- When to merge points //- When to merge points
const scalar mergeDistance_; const scalar mergeDistance_;
//- Whether face might be cut //- Whether face might be cut
List<cellCutType> faceCutType_; List<cellCutType> faceCutType_;
@ -135,6 +139,15 @@ class isoSurface
//- For every unmerged triangle point the point in the triSurface //- For every unmerged triangle point the point in the triSurface
labelList triPointMergeMap_; 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 // Private Member Functions
@ -193,20 +206,8 @@ class isoSurface
const scalarField& pVals const scalarField& pVals
); );
static labelPair findCommonPoints
(
const labelledTri&,
const labelledTri&
);
static point calcCentre(const triSurface&); 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 //- Determine per cc whether all near cuts can be snapped to single
// point. // point.
void calcSnappedCc void calcSnappedCc
@ -323,6 +324,17 @@ class isoSurface
DynamicList<label>& triMeshCells DynamicList<label>& triMeshCells
) const; ) 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 triSurface stitchTriPoints
( (
const bool checkDuplicates, const bool checkDuplicates,
@ -331,57 +343,38 @@ class isoSurface
labelList& triMap // merged to unmerged triangle labelList& triMap // merged to unmerged triangle
) const; ) const;
//- Trim triangle to planes
static void trimToPlanes
(
const PtrList<plane>& planes,
const triPointRef& tri,
DynamicList<point>& newTriPoints
);
//- Trim all triangles to box
static void trimToBox
(
const treeBoundBox& bb,
DynamicList<point>& triPoints,
DynamicList<label>& triMeshCells
);
//- Trim all triangles to box. Determine interpolation
// for existing and new points
static void trimToBox
(
const treeBoundBox& bb,
DynamicList<point>& triPoints,
DynamicList<label>& triMap,
labelList& triPointMap,
labelList& interpolatedPoints,
List<FixedList<label, 3> >& interpolatedOldPoints,
List<FixedList<scalar, 3> >& interpolationWeights
);
//- Check single triangle for (topological) validity //- Check single triangle for (topological) validity
static bool validTri(const triSurface&, const label); 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
(
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
(
const FixedList<label, 3>& fEdges,
const labelList& edgeFace1
);
//- 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
);
static triSurface subsetMesh static triSurface subsetMesh
( (
const triSurface& s, const triSurface& s,
@ -392,6 +385,10 @@ class isoSurface
public: public:
//- Declare friendship with isoSurfaceCell to share some functionality
friend class isoSurfaceCell;
//- Runtime type information //- Runtime type information
TypeName("isoSurface"); TypeName("isoSurface");
@ -407,24 +404,19 @@ public:
const scalarField& pointIsoVals, const scalarField& pointIsoVals,
const scalar iso, const scalar iso,
const bool regularise, const bool regularise,
const boundBox& bounds = boundBox::greatBox,
const scalar mergeTol = 1e-6 // fraction of bounding box const scalar mergeTol = 1e-6 // fraction of bounding box
); );
// Member Functions // Member Functions
//- For every face original cell in mesh //- For every triangle the original cell in mesh
const labelList& meshCells() const const labelList& meshCells() const
{ {
return meshCells_; 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 //- Interpolates cCoords,pCoords. Uses the references to the original
// fields used to create the iso surface. // fields used to create the iso surface.
template<class Type> template<class Type>

279
src/sampling/sampledSurface/isoSurface/isoSurfaceCell.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -30,6 +30,9 @@ License
#include "tetMatcher.H" #include "tetMatcher.H"
#include "syncTools.H" #include "syncTools.H"
#include "addToRunTimeSelectionTable.H" #include "addToRunTimeSelectionTable.H"
#include "Time.H"
#include "triPoints.H"
#include "isoSurface.H"
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // // * * * * * * * * * * * * * * 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. // Checks if triangle is connected through edgeI only.
bool Foam::isoSurfaceCell::danglingTriangle bool Foam::isoSurfaceCell::danglingTriangle
( (
@ -1517,6 +1382,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
const scalarField& pVals, const scalarField& pVals,
const scalar iso, const scalar iso,
const bool regularise, const bool regularise,
const boundBox& bounds,
const scalar mergeTol const scalar mergeTol
) )
: :
@ -1524,6 +1390,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
cVals_(cVals), cVals_(cVals),
pVals_(pVals), pVals_(pVals),
iso_(iso), iso_(iso),
bounds_(bounds),
mergeDistance_(mergeTol*mesh.bounds().mag()) mergeDistance_(mergeTol*mesh.bounds().mag())
{ {
if (debug) if (debug)
@ -1607,56 +1474,104 @@ Foam::isoSurfaceCell::isoSurfaceCell
} }
DynamicList<point> triPoints(nCutCells_);
DynamicList<label> triMeshCells(nCutCells_);
generateTriPoints
(
cVals,
pVals,
mesh_.cellCentres(),
mesh_.points(),
snappedPoints,
snappedCc,
snappedPoint,
triPoints,
triMeshCells
);
if (debug)
{ {
Pout<< "isoSurfaceCell : generated " << triMeshCells.size() DynamicList<point> triPoints(nCutCells_);
<< " unmerged triangles." << endl; DynamicList<label> triMeshCells(nCutCells_);
}
// Merge points and compact out non-valid triangles generateTriPoints
labelList triMap; // merged to unmerged triangle
triSurface::operator=
(
stitchTriPoints
( (
regularise, // check for duplicate tris cVals,
pVals,
mesh_.cellCentres(),
mesh_.points(),
snappedPoints,
snappedCc,
snappedPoint,
triPoints, triPoints,
triPointMergeMap_, // unmerged to merged point triMeshCells
triMap );
)
);
if (debug) if (debug)
{ {
Pout<< "isoSurfaceCell : generated " << triMap.size() Pout<< "isoSurfaceCell : generated " << triMeshCells.size()
<< " merged triangles." << endl; << " unmerged triangles." << endl;
}
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=
(
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 (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]];
}
} }
meshCells_.setSize(triMap.size());
forAll(triMap, i)
{
meshCells_[i] = triMeshCells[triMap[i]];
}
if (debug) if (debug)
{ {

62
src/sampling/sampledSurface/isoSurface/isoSurfaceCell.H
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -48,6 +48,7 @@ SourceFiles
#include "labelPair.H" #include "labelPair.H"
#include "pointIndexHit.H" #include "pointIndexHit.H"
#include "PackedBoolList.H" #include "PackedBoolList.H"
#include "boundBox.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@ -55,6 +56,7 @@ namespace Foam
{ {
class polyMesh; class polyMesh;
class plane;
/*---------------------------------------------------------------------------*\ /*---------------------------------------------------------------------------*\
Class isoSurfaceCell Declaration Class isoSurfaceCell Declaration
@ -91,6 +93,9 @@ class isoSurfaceCell
//- isoSurfaceCell value //- isoSurfaceCell value
const scalar iso_; const scalar iso_;
//- Optional bounds
const boundBox bounds_;
//- When to merge points //- When to merge points
const scalar mergeDistance_; const scalar mergeDistance_;
@ -106,6 +111,15 @@ class isoSurfaceCell
//- For every unmerged triangle point the point in the triSurface //- For every unmerged triangle point the point in the triSurface
labelList triPointMergeMap_; 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 // Private Member Functions
@ -214,19 +228,15 @@ class isoSurfaceCell
void generateTriPoints void generateTriPoints
( (
const DynamicList<Type>& snapped, const DynamicList<Type>& snapped,
const scalar isoVal0,
const scalar s0, const scalar s0,
const Type& p0, const Type& p0,
const label p0Index, const label p0Index,
const scalar isoVal1,
const scalar s1, const scalar s1,
const Type& p1, const Type& p1,
const label p1Index, const label p1Index,
const scalar isoVal2,
const scalar s2, const scalar s2,
const Type& p2, const Type& p2,
const label p2Index, const label p2Index,
const scalar isoVal3,
const scalar s3, const scalar s3,
const Type& p3, const Type& p3,
const label p3Index, const label p3Index,
@ -271,27 +281,6 @@ class isoSurfaceCell
Map<labelList>& edgeFacesRest Map<labelList>& edgeFacesRest
) const; ) 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? //- Is triangle (given by 3 edges) not fully connected?
static bool danglingTriangle static bool danglingTriangle
( (
@ -317,8 +306,6 @@ class isoSurfaceCell
labelList& newToOldPoints labelList& newToOldPoints
); );
//- Combine all triangles inside a cell into a minimal triangulation
void combineCellTriangles();
public: public:
@ -336,6 +323,7 @@ public:
const scalarField& pointValues, const scalarField& pointValues,
const scalar iso, const scalar iso,
const bool regularise, const bool regularise,
const boundBox& bounds = boundBox::greatBox,
const scalar mergeTol = 1e-6 // fraction of bounding box const scalar mergeTol = 1e-6 // fraction of bounding box
); );
@ -348,24 +336,6 @@ public:
return meshCells_; 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. //- Interpolates cCoords,pCoords.
template<class Type> template<class Type>
tmp<Field<Type> > interpolate tmp<Field<Type> > interpolate

354
src/sampling/sampledSurface/isoSurface/isoSurfaceCellTemplates.C
View File

@ -26,6 +26,7 @@ License
#include "isoSurfaceCell.H" #include "isoSurfaceCell.H"
#include "polyMesh.H" #include "polyMesh.H"
#include "tetMatcher.H" #include "tetMatcher.H"
#include "isoSurface.H"
// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
@ -76,22 +77,18 @@ void Foam::isoSurfaceCell::generateTriPoints
( (
const DynamicList<Type>& snapped, const DynamicList<Type>& snapped,
const scalar isoVal0,
const scalar s0, const scalar s0,
const Type& p0, const Type& p0,
const label p0Index, const label p0Index,
const scalar isoVal1,
const scalar s1, const scalar s1,
const Type& p1, const Type& p1,
const label p1Index, const label p1Index,
const scalar isoVal2,
const scalar s2, const scalar s2,
const Type& p2, const Type& p2,
const label p2Index, const label p2Index,
const scalar isoVal3,
const scalar s3, const scalar s3,
const Type& p3, const Type& p3,
const label p3Index, const label p3Index,
@ -124,160 +121,196 @@ void Foam::isoSurfaceCell::generateTriPoints
case 0x0F: case 0x0F:
break; break;
case 0x0E:
case 0x01: case 0x01:
case 0x0E:
{ {
// 0 is common point. Orient such that normal points in positive pts.append
// gradient direction (
if (isoVal0 >= isoVal1) 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,s1,p1,p1Index)); // Flip normals
pts.append(generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index)); label sz = pts.size();
pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); Swap(pts[sz-2], pts[sz-1]);
}
else
{
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));
} }
} }
break; break;
case 0x0D:
case 0x02: case 0x02:
case 0x0D:
{ {
// 1 is common point pts.append
if (isoVal1 >= isoVal0) (
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,s0,p0,p0Index)); // Flip normals
pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)); label sz = pts.size();
pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index)); Swap(pts[sz-2], pts[sz-1]);
}
else
{
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));
} }
} }
break; break;
case 0x0C:
case 0x03: case 0x03:
case 0x0C:
{ {
Type s02 = generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index); Type p0p2 =
Type s13 = generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index); generatePoint(snapped,s0,p0,p0Index,s2,p2,p2Index);
Type p1p3 =
generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index);
if (isoVal0 >= isoVal3) 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(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); // Flip normals
pts.append(s02); label sz = pts.size();
pts.append(s13); Swap(pts[sz-5], pts[sz-4]);
pts.append(s13); Swap(pts[sz-2], pts[sz-1]);
pts.append(generatePoint(snapped,s1,p1,p1Index,s2,p2,p2Index));
pts.append(s02);
}
else
{
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);
} }
} }
break; break;
case 0x0B:
case 0x04: case 0x04:
case 0x0B:
{ {
// 2 is common point pts.append
if (isoVal2 >= isoVal0) (
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)); // Flip normals
pts.append(generatePoint(snapped,s2,p2,p2Index,s1,p1,p1Index)); label sz = pts.size();
pts.append(generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index)); Swap(pts[sz-2], pts[sz-1]);
}
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));
} }
} }
break; break;
case 0x0A:
case 0x05: case 0x05:
case 0x0A:
{ {
Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index); Type p0p1 =
Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index); generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
Type p2p3 =
generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
if (isoVal3 >= isoVal0) 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); // Flip normals
pts.append(s23); label sz = pts.size();
pts.append(generatePoint(snapped,s0,p0,p0Index,s3,p3,p3Index)); Swap(pts[sz-5], pts[sz-4]);
pts.append(s01); Swap(pts[sz-2], pts[sz-1]);
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);
} }
} }
break; break;
case 0x09:
case 0x06: case 0x06:
case 0x09:
{ {
Type s01 = generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index); Type p0p1 =
Type s23 = generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index); generatePoint(snapped,s0,p0,p0Index,s1,p1,p1Index);
Type p2p3 =
generatePoint(snapped,s2,p2,p2Index,s3,p3,p3Index);
if (isoVal3 >= isoVal1) 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); // Flip normals
pts.append(generatePoint(snapped,s1,p1,p1Index,s3,p3,p3Index)); label sz = pts.size();
pts.append(s23); Swap(pts[sz-5], pts[sz-4]);
pts.append(s01); Swap(pts[sz-2], pts[sz-1]);
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);
} }
} }
break; break;
case 0x07:
case 0x08: case 0x08:
case 0x07:
{ {
// 3 is common point pts.append
if (isoVal3 >= isoVal0) (
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)); // Flip normals
pts.append(generatePoint(snapped,s3,p3,p3Index,s2,p2,p2Index)); label sz = pts.size();
pts.append(generatePoint(snapped,s3,p3,p3Index,s1,p1,p1Index)); Swap(pts[sz-2], pts[sz-1]);
}
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));
} }
} }
break; break;
@ -341,22 +374,18 @@ void Foam::isoSurfaceCell::generateTriPoints
( (
snappedPoints, snappedPoints,
pVals_[f0[1]],
pVals[f0[1]], pVals[f0[1]],
pCoords[f0[1]], pCoords[f0[1]],
snappedPoint[f0[1]], snappedPoint[f0[1]],
pVals_[f0[0]],
pVals[f0[0]], pVals[f0[0]],
pCoords[f0[0]], pCoords[f0[0]],
snappedPoint[f0[0]], snappedPoint[f0[0]],
pVals_[f0[2]],
pVals[f0[2]], pVals[f0[2]],
pCoords[f0[2]], pCoords[f0[2]],
snappedPoint[f0[2]], snappedPoint[f0[2]],
pVals_[oppositeI],
pVals[oppositeI], pVals[oppositeI],
pCoords[oppositeI], pCoords[oppositeI],
snappedPoint[oppositeI], snappedPoint[oppositeI],
@ -370,22 +399,18 @@ void Foam::isoSurfaceCell::generateTriPoints
( (
snappedPoints, snappedPoints,
pVals_[f0[0]],
pVals[f0[0]], pVals[f0[0]],
pCoords[f0[0]], pCoords[f0[0]],
snappedPoint[f0[0]], snappedPoint[f0[0]],
pVals_[f0[1]],
pVals[f0[1]], pVals[f0[1]],
pCoords[f0[1]], pCoords[f0[1]],
snappedPoint[f0[1]], snappedPoint[f0[1]],
pVals_[f0[2]],
pVals[f0[2]], pVals[f0[2]],
pCoords[f0[2]], pCoords[f0[2]],
snappedPoint[f0[2]], snappedPoint[f0[2]],
pVals_[oppositeI],
pVals[oppositeI], pVals[oppositeI],
pCoords[oppositeI], pCoords[oppositeI],
snappedPoint[oppositeI], snappedPoint[oppositeI],
@ -411,7 +436,6 @@ void Foam::isoSurfaceCell::generateTriPoints
} }
label fp = f.fcIndex(fp0); label fp = f.fcIndex(fp0);
for (label i = 2; i < f.size(); i++) for (label i = 2; i < f.size(); i++)
{ {
label nextFp = f.fcIndex(fp); label nextFp = f.fcIndex(fp);
@ -425,22 +449,18 @@ void Foam::isoSurfaceCell::generateTriPoints
( (
snappedPoints, snappedPoints,
pVals_[tri[1]],
pVals[tri[1]], pVals[tri[1]],
pCoords[tri[1]], pCoords[tri[1]],
snappedPoint[tri[1]], snappedPoint[tri[1]],
pVals_[tri[0]],
pVals[tri[0]], pVals[tri[0]],
pCoords[tri[0]], pCoords[tri[0]],
snappedPoint[tri[0]], snappedPoint[tri[0]],
pVals_[tri[2]],
pVals[tri[2]], pVals[tri[2]],
pCoords[tri[2]], pCoords[tri[2]],
snappedPoint[tri[2]], snappedPoint[tri[2]],
cVals_[cellI],
cVals[cellI], cVals[cellI],
cCoords[cellI], cCoords[cellI],
snappedCc[cellI], snappedCc[cellI],
@ -454,22 +474,18 @@ void Foam::isoSurfaceCell::generateTriPoints
( (
snappedPoints, snappedPoints,
pVals_[tri[0]],
pVals[tri[0]], pVals[tri[0]],
pCoords[tri[0]], pCoords[tri[0]],
snappedPoint[tri[0]], snappedPoint[tri[0]],
pVals_[tri[1]],
pVals[tri[1]], pVals[tri[1]],
pCoords[tri[1]], pCoords[tri[1]],
snappedPoint[tri[1]], snappedPoint[tri[1]],
pVals_[tri[2]],
pVals[tri[2]], pVals[tri[2]],
pCoords[tri[2]], pCoords[tri[2]],
snappedPoint[tri[2]], snappedPoint[tri[2]],
cVals_[cellI],
cVals[cellI], cVals[cellI],
cCoords[cellI], cCoords[cellI],
snappedCc[cellI], snappedCc[cellI],
@ -495,7 +511,7 @@ void Foam::isoSurfaceCell::generateTriPoints
if (countNotFoundTets > 0) if (countNotFoundTets > 0)
{ {
WarningIn("Foam::isoSurfaceCell::generateTriPoints") WarningIn("Foam::isoSurfaceCell::generateTriPoints(..)")
<< "Could not find " << countNotFoundTets << "Could not find " << countNotFoundTets
<< " tet base points, which may lead to inverted triangles." << " tet base points, which may lead to inverted triangles."
<< endl; << endl;
@ -510,13 +526,11 @@ template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::tmp<Foam::Field<Type> >
Foam::isoSurfaceCell::interpolate Foam::isoSurfaceCell::interpolate
( (
const scalarField& cVals,
const scalarField& pVals,
const Field<Type>& cCoords, const Field<Type>& cCoords,
const Field<Type>& pCoords const Field<Type>& pCoords
) const ) const
{ {
DynamicList<Type> triPoints(nCutCells_); DynamicList<Type> triPoints(3*nCutCells_);
DynamicList<label> triMeshCells(nCutCells_); DynamicList<label> triMeshCells(nCutCells_);
// Dummy snap data // Dummy snap data
@ -524,59 +538,6 @@ Foam::isoSurfaceCell::interpolate
labelList snappedCc(mesh_.nCells(), -1); labelList snappedCc(mesh_.nCells(), -1);
labelList snappedPoint(mesh_.nPoints(), -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 generateTriPoints
( (
cVals_, cVals_,
@ -593,22 +554,15 @@ Foam::isoSurfaceCell::interpolate
triMeshCells triMeshCells
); );
return isoSurface::interpolate
// One value per point (
tmp<Field<Type> > tvalues(new Field<Type>(points().size())); points().size(),
Field<Type>& values = tvalues(); triPointMergeMap_,
interpolatedPoints_,
forAll(triPoints, i) interpolatedOldPoints_,
{ interpolationWeights_,
label mergedPointI = triPointMergeMap_[i]; triPoints
);
if (mergedPointI >= 0)
{
values[mergedPointI] = triPoints[i];
}
}
return tvalues;
} }

320
src/sampling/sampledSurface/isoSurface/isoSurfaceTemplates.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -120,7 +120,7 @@ Foam::isoSurface::adaptPatchFields
{ {
fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&> fvPatchField<Type>& pfld = const_cast<fvPatchField<Type>&>
( (
fld.boundaryField()[patchI] sliceFld.boundaryField()[patchI]
); );
const scalarField& w = mesh.weights().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> template<class Type>
void Foam::isoSurface::generateTriPoints void Foam::isoSurface::generateTriPoints
( (
@ -240,8 +242,8 @@ void Foam::isoSurface::generateTriPoints
case 0x0F: case 0x0F:
break; break;
case 0x0E:
case 0x01: case 0x01:
case 0x0E:
points.append points.append
( (
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1) 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) 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; break;
case 0x0D:
case 0x02: case 0x02:
case 0x0D:
points.append points.append
( (
generatePoint(s1,p1,hasSnap1,snapP1,s0,p0,hasSnap0,snapP0) 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) 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; break;
case 0x0C:
case 0x03: case 0x03:
{ case 0x0C:
Type tp1 = {
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2); Type p0p2 =
Type tp2 = generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3); Type p1p3 =
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3);
points.append points.append
( (
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3) generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
); );
points.append(tp1); points.append(p1p3);
points.append(tp2); points.append(p0p2);
points.append(tp2);
points.append points.append(p1p3);
( points.append
generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2) (
); generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
points.append(tp1); );
} 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; break;
case 0x0B:
case 0x04: case 0x04:
{ case 0x0B:
points.append {
( points.append
generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0) (
); generatePoint(s2,p2,hasSnap2,snapP2,s0,p0,hasSnap0,snapP0)
points.append );
( points.append
generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1) (
); generatePoint(s2,p2,hasSnap2,snapP2,s1,p1,hasSnap1,snapP1)
points.append );
( points.append
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3) (
); 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; break;
case 0x0A:
case 0x05: case 0x05:
{ case 0x0A:
Type tp0 = {
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1); Type p0p1 =
Type tp1 = generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3); Type p2p3 =
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
points.append(tp0); points.append(p0p1);
points.append(tp1); points.append(p2p3);
points.append points.append
( (
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3) generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
); );
points.append(tp0);
points.append points.append(p0p1);
( points.append
generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2) (
); generatePoint(s1,p1,hasSnap1,snapP1,s2,p2,hasSnap2,snapP2)
points.append(tp1); );
} 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; break;
case 0x09:
case 0x06: case 0x06:
{ case 0x09:
Type tp0 = {
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1); Type p0p1 =
Type tp1 = generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3); Type p2p3 =
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
points.append(tp0); points.append(p0p1);
points.append points.append
( (
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3) generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
); );
points.append(tp1); points.append(p2p3);
points.append(tp0);
points.append points.append(p0p1);
( points.append(p2p3);
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2) points.append
); (
points.append(tp1); 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; break;
case 0x07:
case 0x08: case 0x08:
case 0x07:
points.append points.append
( (
generatePoint(s3,p3,hasSnap3,snapP3,s0,p0,hasSnap0,snapP0) 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) 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; 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> template<class Type>
Foam::tmp<Foam::Field<Type> > Foam::tmp<Foam::Field<Type> >
Foam::isoSurface::interpolate Foam::isoSurface::interpolate
@ -707,7 +820,7 @@ Foam::isoSurface::interpolate
> > c2(adaptPatchFields(cCoords)); > > c2(adaptPatchFields(cCoords));
DynamicList<Type> triPoints(nCutCells_); DynamicList<Type> triPoints(3*nCutCells_);
DynamicList<label> triMeshCells(nCutCells_); DynamicList<label> triMeshCells(nCutCells_);
// Dummy snap data // Dummy snap data
@ -731,52 +844,15 @@ Foam::isoSurface::interpolate
triMeshCells triMeshCells
); );
return interpolate
// One value per point
tmp<Field<Type> > tvalues
( (
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;
} }

125
src/sampling/sampledSurface/isoSurface/sampledIsoSurface.C
View File

@ -122,18 +122,52 @@ void Foam::sampledIsoSurface::getIsoFields() const
// Get pointField // 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()) if (!subMeshPtr_.valid())
{ {
word pointFldName = "volPointInterpolate(" + isoField_ + ')'; const word pointFldName =
"volPointInterpolate_"
+ type()
+ "("
+ isoField_
+ ')';
if (fvm.foundObject<pointScalarField>(pointFldName)) if (fvm.foundObject<pointScalarField>(pointFldName))
{ {
if (debug) if (debug)
{ {
Info<< "sampledIsoSurface::getIsoFields() : lookup pointField " Info<< "sampledIsoSurface::getIsoFields() :"
<< pointFldName << endl; << " 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 else
{ {
@ -141,33 +175,23 @@ void Foam::sampledIsoSurface::getIsoFields() const
if (debug) if (debug)
{ {
Info<< "sampledIsoSurface::getIsoFields() : " Info<< "sampledIsoSurface::getIsoFields() :"
<< "checking pointField " << pointFldName << " creating and storing pointField "
<< " for same time " << fvm.time().timeName() << pointFldName << " for time "
<< endl; << fvm.time().timeName() << endl;
} }
if tmp<pointScalarField> tpfld
( (
storedPointFieldPtr_.empty() volPointInterpolation::New(fvm).interpolate
|| (fvm.time().timeName() != storedPointFieldPtr_().instance())
)
{
if (debug)
{
Info<< "sampledIsoSurface::getIsoFields() :"
<< " interpolating volField " << volFieldPtr_->name()
<< " to get pointField " << pointFldName << endl;
}
storedPointFieldPtr_.reset
( (
volPointInterpolation::New(fvm) *volFieldPtr_,
.interpolate(*volFieldPtr_).ptr() pointFldName,
); false
storedPointFieldPtr_->checkOut(); )
pointFieldPtr_ = storedPointFieldPtr_.operator->(); );
} pointFieldPtr_ = tpfld.ptr();
const_cast<pointScalarField*>(pointFieldPtr_)->store();
} }
@ -233,8 +257,10 @@ void Foam::sampledIsoSurface::getIsoFields() const
// Pointfield on submesh // Pointfield on submesh
word pointFldName = const word pointFldName =
"volPointInterpolate(" "volPointInterpolate_"
+ type()
+ "("
+ volSubFieldPtr_->name() + volSubFieldPtr_->name()
+ ')'; + ')';
@ -245,11 +271,28 @@ void Foam::sampledIsoSurface::getIsoFields() const
Info<< "sampledIsoSurface::getIsoFields() :" Info<< "sampledIsoSurface::getIsoFields() :"
<< " submesh lookup pointField " << pointFldName << endl; << " submesh lookup pointField " << pointFldName << endl;
} }
storedPointSubFieldPtr_.clear(); const pointScalarField& pfld = subFvm.lookupObject<pointScalarField>
pointSubFieldPtr_ = &subFvm.lookupObject<pointScalarField>
( (
pointFldName 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 else
{ {
@ -260,15 +303,15 @@ void Foam::sampledIsoSurface::getIsoFields() const
<< volSubFieldPtr_->name() << volSubFieldPtr_->name()
<< " to get submesh pointField " << pointFldName << endl; << " to get submesh pointField " << pointFldName << endl;
} }
storedPointSubFieldPtr_.reset tmp<pointScalarField> tpfld
( (
volPointInterpolation::New volPointInterpolation::New
( (
subFvm subFvm
).interpolate(*volSubFieldPtr_).ptr() ).interpolate(*volSubFieldPtr_)
); );
storedPointSubFieldPtr_->checkOut(); pointSubFieldPtr_ = tpfld.ptr();
pointSubFieldPtr_ = storedPointSubFieldPtr_.operator->(); const_cast<pointScalarField*>(pointSubFieldPtr_)->store();
} }
@ -349,28 +392,34 @@ bool Foam::sampledIsoSurface::updateGeometry() const
if (subMeshPtr_.valid()) if (subMeshPtr_.valid())
{ {
const volScalarField& vfld = *volSubFieldPtr_;
surfPtr_.reset surfPtr_.reset
( (
new isoSurface new isoSurface
( (
*volSubFieldPtr_, vfld,
*pointSubFieldPtr_, *pointSubFieldPtr_,
isoVal_, isoVal_,
regularise_, regularise_,
bounds_,
mergeTol_ mergeTol_
) )
); );
} }
else else
{ {
const volScalarField& vfld = *volFieldPtr_;
surfPtr_.reset surfPtr_.reset
( (
new isoSurface new isoSurface
( (
*volFieldPtr_, vfld,
*pointFieldPtr_, *pointFieldPtr_,
isoVal_, isoVal_,
regularise_, regularise_,
bounds_,
mergeTol_ mergeTol_
) )
); );
@ -412,6 +461,7 @@ Foam::sampledIsoSurface::sampledIsoSurface
sampledSurface(name, mesh, dict), sampledSurface(name, mesh, dict),
isoField_(dict.lookup("isoField")), isoField_(dict.lookup("isoField")),
isoVal_(readScalar(dict.lookup("isoValue"))), isoVal_(readScalar(dict.lookup("isoValue"))),
bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)), mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
regularise_(dict.lookupOrDefault("regularise", true)), regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", false)), average_(dict.lookupOrDefault("average", false)),
@ -422,7 +472,6 @@ Foam::sampledIsoSurface::sampledIsoSurface
prevTimeIndex_(-1), prevTimeIndex_(-1),
storedVolFieldPtr_(NULL), storedVolFieldPtr_(NULL),
volFieldPtr_(NULL), volFieldPtr_(NULL),
storedPointFieldPtr_(NULL),
pointFieldPtr_(NULL) pointFieldPtr_(NULL)
{ {
if (!sampledSurface::interpolate()) if (!sampledSurface::interpolate())

5
src/sampling/sampledSurface/isoSurface/sampledIsoSurface.H
View File

@ -63,6 +63,9 @@ class sampledIsoSurface
//- Iso value //- Iso value
const scalar isoVal_; const scalar isoVal_;
//- Optional bounding box to trim triangles against
const boundBox bounds_;
//- Merge tolerance //- Merge tolerance
const scalar mergeTol_; const scalar mergeTol_;
@ -94,7 +97,6 @@ class sampledIsoSurface
mutable const volScalarField* volFieldPtr_; mutable const volScalarField* volFieldPtr_;
//- Cached pointfield //- Cached pointfield
mutable autoPtr<pointScalarField> storedPointFieldPtr_;
mutable const pointScalarField* pointFieldPtr_; mutable const pointScalarField* pointFieldPtr_;
// And on subsetted mesh // And on subsetted mesh
@ -107,7 +109,6 @@ class sampledIsoSurface
mutable const volScalarField* volSubFieldPtr_; mutable const volScalarField* volSubFieldPtr_;
//- Cached pointfield //- Cached pointfield
mutable autoPtr<pointScalarField> storedPointSubFieldPtr_;
mutable const pointScalarField* pointSubFieldPtr_; mutable const pointScalarField* pointSubFieldPtr_;

10
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -148,7 +148,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
cellAvg, cellAvg,
pointFld().internalField(), pointFld().internalField(),
isoVal_, isoVal_,
regularise_ regularise_,
bounds_
); );
const_cast<sampledIsoSurfaceCell&> const_cast<sampledIsoSurfaceCell&>
@ -166,7 +167,8 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
cellFld.internalField(), cellFld.internalField(),
pointFld().internalField(), pointFld().internalField(),
isoVal_, isoVal_,
regularise_ regularise_,
bounds_
); );
const_cast<sampledIsoSurfaceCell&> const_cast<sampledIsoSurfaceCell&>
@ -185,6 +187,7 @@ bool Foam::sampledIsoSurfaceCell::updateGeometry() const
<< " average : " << average_ << nl << " average : " << average_ << nl
<< " isoField : " << isoField_ << nl << " isoField : " << isoField_ << nl
<< " isoValue : " << isoVal_ << nl << " isoValue : " << isoVal_ << nl
<< " bounds : " << bounds_ << nl
<< " points : " << points().size() << nl << " points : " << points().size() << nl
<< " tris : " << triSurface::size() << nl << " tris : " << triSurface::size() << nl
<< " cut cells : " << meshCells_.size() << endl; << " cut cells : " << meshCells_.size() << endl;
@ -206,6 +209,7 @@ Foam::sampledIsoSurfaceCell::sampledIsoSurfaceCell
sampledSurface(name, mesh, dict), sampledSurface(name, mesh, dict),
isoField_(dict.lookup("isoField")), isoField_(dict.lookup("isoField")),
isoVal_(readScalar(dict.lookup("isoValue"))), isoVal_(readScalar(dict.lookup("isoValue"))),
bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
regularise_(dict.lookupOrDefault("regularise", true)), regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", true)), average_(dict.lookupOrDefault("average", true)),
zoneKey_(keyType::null), zoneKey_(keyType::null),

3
src/sampling/sampledSurface/isoSurface/sampledIsoSurfaceCell.H
View File

@ -62,6 +62,9 @@ class sampledIsoSurfaceCell
//- Iso value //- Iso value
const scalar isoVal_; const scalar isoVal_;
//- Optional bounding box to trim triangles against
const boundBox bounds_;
//- Whether to coarse //- Whether to coarse
const Switch regularise_; const Switch regularise_;

4
src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2012 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -230,6 +230,7 @@ void Foam::sampledCuttingPlane::createGeometry()
pointDistance_, pointDistance_,
0.0, 0.0,
regularise_, regularise_,
bounds_,
mergeTol_ mergeTol_
) )
//new isoSurfaceCell //new isoSurfaceCell
@ -262,6 +263,7 @@ Foam::sampledCuttingPlane::sampledCuttingPlane
: :
sampledSurface(name, mesh, dict), sampledSurface(name, mesh, dict),
plane_(dict), plane_(dict),
bounds_(dict.lookupOrDefault("bounds", boundBox::greatBox)),
mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)), mergeTol_(dict.lookupOrDefault("mergeTol", 1e-6)),
regularise_(dict.lookupOrDefault("regularise", true)), regularise_(dict.lookupOrDefault("regularise", true)),
average_(dict.lookupOrDefault("average", false)), average_(dict.lookupOrDefault("average", false)),

3
src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.H
View File

@ -60,6 +60,9 @@ class sampledCuttingPlane
//- Plane //- Plane
const plane plane_; const plane plane_;
//- Optional bounding box to trim triangles against
const boundBox bounds_;
//- Merge tolerance //- Merge tolerance
const scalar mergeTol_; const scalar mergeTol_;

8
src/sampling/sampledSurface/sampledPatch/sampledPatch.C
View File

@ -2,7 +2,7 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
@ -84,11 +84,7 @@ const Foam::labelList& Foam::sampledPatch::patchIDs() const
{ {
if (patchIDs_.empty()) if (patchIDs_.empty())
{ {
patchIDs_ = mesh().boundaryMesh().patchSet patchIDs_ = mesh().boundaryMesh().patchSet(patchNames_).sortedToc();
(
patchNames_,
false
).sortedToc();
} }
return patchIDs_; return patchIDs_;
} }

18
src/sampling/sampledSurface/sampledSurfaces/sampledSurfaces.C
View File

@ -2,8 +2,8 @@
========= | ========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration | \\ / O peration |
\\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation \\ / A nd | Copyright (C) 2011-2014 OpenFOAM Foundation
\\/ M anipulation | Copyright (C) 2015 OpenCFD Ltd \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -162,19 +162,7 @@ void Foam::sampledSurfaces::write()
const label nFields = classifyFields(); const label nFields = classifyFields();
if (Pstream::master()) // write geometry first if required,
{
if (debug)
{
Pout<< "Creating directory "
<< outputPath_/obr_.time().timeName() << nl << endl;
}
mkDir(outputPath_/obr_.time().timeName());
}
// Write geometry first if required,
// or when no fields would otherwise be written // or when no fields would otherwise be written
if (nFields == 0 || formatter_->separateGeometry()) if (nFields == 0 || formatter_->separateGeometry())
{ {

19
src/sampling/sampledSurface/sampledTriSurfaceMesh/sampledTriSurfaceMesh.C
View File

@ -644,8 +644,25 @@ bool Foam::sampledTriSurfaceMesh::update()
return false; 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. // Mesh search engine, no triangulation of faces.
meshSearch meshSearcher(mesh(), polyMesh::FACE_PLANES); meshSearch meshSearcher(mesh(), bb, polyMesh::FACE_PLANES);
return update(meshSearcher); return update(meshSearcher);
} }