mirror of
https://develop.openfoam.com/Development/openfoam.git
synced 2025-11-28 03:28:01 +00:00
Merge branch 'feature-sampledSurfaces' into 'develop'
Feature sampled surfaces Bounds handling for sampledSurfaces only (so not for streamlines) See merge request !13
This commit is contained in:
Andrew Heather
@ -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)
|
||||
// - 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));
|
||||
}
|
||||
|
||||
);
|
||||
@ -252,6 +310,8 @@ surfaces
|
||||
// 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
|
||||
// to merge points (default=1e-6)
|
||||
@ -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)
|
||||
}
|
||||
@ -284,6 +347,8 @@ surfaces
|
||||
// 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
|
||||
// to merge points (default=1e-6)
|
||||
@ -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)
|
||||
}
|
||||
|
||||
242
src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C
Normal file
242
src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C
Normal 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]);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// ************************************************************************* //
|
||||
@ -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,13 +83,14 @@ inline Ostream& operator<<
|
||||
template<class Point, class PointRef>
|
||||
class triangle
|
||||
{
|
||||
// Private data
|
||||
|
||||
PointRef a_, b_, c_;
|
||||
|
||||
|
||||
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
|
||||
{
|
||||
@ -95,11 +100,81 @@ public:
|
||||
};
|
||||
|
||||
|
||||
// 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_;
|
||||
|
||||
|
||||
// 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:
|
||||
|
||||
// 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
|
||||
|
||||
// ************************************************************************* //
|
||||
|
||||
@ -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>
|
||||
|
||||
@ -33,8 +33,8 @@ SeeAlso
|
||||
|
||||
\*---------------------------------------------------------------------------*/
|
||||
|
||||
#ifndef plane_H
|
||||
#define plane_H
|
||||
#ifndef planeExtrusion_H
|
||||
#define planeExtrusion_H
|
||||
|
||||
#include "linearNormal.H"
|
||||
|
||||
|
||||
@ -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();
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
|
||||
@ -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),
|
||||
|
||||
@ -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
|
||||
);
|
||||
|
||||
|
||||
|
||||
@ -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,6 +1053,220 @@ Foam::triSurface Foam::isoSurface::stitchTriPoints
|
||||
}
|
||||
|
||||
|
||||
void Foam::isoSurface::trimToPlanes
|
||||
(
|
||||
const PtrList<plane>& planes,
|
||||
const triPointRef& tri,
|
||||
DynamicList<point>& newTriPoints
|
||||
)
|
||||
{
|
||||
// Buffer for generated triangles
|
||||
DynamicList<triPoints> insideTrisA;
|
||||
storeOp insideOpA(insideTrisA);
|
||||
|
||||
// Buffer for generated triangles
|
||||
DynamicList<triPoints> insideTrisB;
|
||||
storeOp insideOpB(insideTrisB);
|
||||
|
||||
triPointRef::dummyOp dop;
|
||||
|
||||
// Store starting triangle in insideTrisA
|
||||
insideOpA(triPoints(tri.a(), tri.b(), tri.c()));
|
||||
|
||||
|
||||
bool useA = true;
|
||||
|
||||
forAll(planes, faceI)
|
||||
{
|
||||
const plane& pl = planes[faceI];
|
||||
|
||||
if (useA)
|
||||
{
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
// Transfer
|
||||
if (useA)
|
||||
{
|
||||
forAll(insideOpA.tris_, i)
|
||||
{
|
||||
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]);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::isoSurface::trimToBox
|
||||
(
|
||||
const treeBoundBox& bb,
|
||||
DynamicList<point>& triPoints,
|
||||
DynamicList<label>& triMap // trimmed to original tris
|
||||
)
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "isoSurface : trimming to " << bb << endl;
|
||||
}
|
||||
|
||||
// Generate inwards pointing planes
|
||||
PtrList<plane> planes(6);
|
||||
const pointField pts(bb.treeBoundBox::points());
|
||||
forAll(treeBoundBox::faces, faceI)
|
||||
{
|
||||
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;
|
||||
|
||||
DynamicList<point> newTriPoints(triPoints.size()/16);
|
||||
triMap.setCapacity(nTris/16);
|
||||
|
||||
label vertI = 0;
|
||||
for (label triI = 0; triI < nTris; triI++)
|
||||
{
|
||||
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
|
||||
);
|
||||
|
||||
label nCells = (newTriPoints.size()-oldNPoints)/3;
|
||||
for (label i = 0; i < nCells; i++)
|
||||
{
|
||||
triMap.append(triI);
|
||||
}
|
||||
}
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "isoSurface : trimmed from " << nTris
|
||||
<< " down to " << triMap.size()
|
||||
<< " triangles." << endl;
|
||||
}
|
||||
|
||||
triPoints.transfer(newTriPoints);
|
||||
}
|
||||
|
||||
|
||||
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);
|
||||
|
||||
// Trim triPoints, return map
|
||||
trimToBox(bb, triPoints, triMap);
|
||||
|
||||
|
||||
// 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);
|
||||
|
||||
|
||||
triPointMap.setSize(triPoints.size());
|
||||
forAll(triMap, triI)
|
||||
{
|
||||
label oldTriI = triMap[triI];
|
||||
|
||||
// Find point correspondence. Assumes coordinate is bit-copy.
|
||||
for (label i = 0; i < 3; i++)
|
||||
{
|
||||
label pointI = 3*triI+i;
|
||||
const point& pt = triPoints[pointI];
|
||||
|
||||
// Compare to old-triangle's points
|
||||
label matchPointI = -1;
|
||||
for (label j = 0; j < 3; j++)
|
||||
{
|
||||
label oldPointI = 3*oldTriI+j;
|
||||
if (pt == oldTriPoints[oldPointI])
|
||||
{
|
||||
matchPointI = oldPointI;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
triPointMap[pointI] = matchPointI;
|
||||
|
||||
// If new point: calculate and store interpolation
|
||||
if (matchPointI == -1)
|
||||
{
|
||||
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);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
interpolatedPoints.transfer(dynInterpolatedPoints);
|
||||
interpolatedOldPoints.transfer(dynInterpolatedOldPoints);
|
||||
interpolationWeights.transfer(dynInterpolationWeights);
|
||||
}
|
||||
|
||||
|
||||
// Does face use valid vertices?
|
||||
bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
|
||||
{
|
||||
@ -1204,437 +1338,6 @@ bool Foam::isoSurface::validTri(const triSurface& surf, const label faceI)
|
||||
}
|
||||
|
||||
|
||||
void Foam::isoSurface::calcAddressing
|
||||
(
|
||||
const triSurface& surf,
|
||||
List<FixedList<label, 3> >& faceEdges,
|
||||
labelList& edgeFace0,
|
||||
labelList& edgeFace1,
|
||||
Map<labelList>& edgeFacesRest
|
||||
) const
|
||||
{
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
// Determine faceEdges
|
||||
faceEdges.setSize(surf.size());
|
||||
edgeI = 0;
|
||||
forAll(surf, triI)
|
||||
{
|
||||
faceEdges[triI][0] = oldToMerged[edgeI++];
|
||||
faceEdges[triI][1] = oldToMerged[edgeI++];
|
||||
faceEdges[triI][2] = oldToMerged[edgeI++];
|
||||
}
|
||||
|
||||
|
||||
// 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);
|
||||
|
||||
forAll(oldToMerged, oldEdgeI)
|
||||
{
|
||||
label triI = oldEdgeI / 3;
|
||||
label edgeI = oldToMerged[oldEdgeI];
|
||||
|
||||
if (edgeFace0[edgeI] == -1)
|
||||
{
|
||||
// First triangle for edge
|
||||
edgeFace0[edgeI] = 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;
|
||||
|
||||
edge e(tri[fp], tri[tri.fcIndex(fp)]);
|
||||
|
||||
label prevTriIndex = -1;
|
||||
|
||||
forAll(prevTri, i)
|
||||
{
|
||||
if (edge(prevTri[i], prevTri[prevTri.fcIndex(i)]) == e)
|
||||
{
|
||||
prevTriIndex = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
if (prevTriIndex == -1)
|
||||
{
|
||||
// Different edge. Store for later.
|
||||
EdgeMap<labelList>::iterator iter = extraEdgeFaces.find(e);
|
||||
|
||||
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();
|
||||
|
||||
edgeFace0.setSize(edgeI + extraEdgeFaces.size());
|
||||
edgeFace1.setSize(edgeI + extraEdgeFaces.size(), -1);
|
||||
|
||||
forAllConstIter(EdgeMap<labelList>, extraEdgeFaces, iter)
|
||||
{
|
||||
const labelList& eFaces = iter();
|
||||
|
||||
// The current edge will become edgeI. Replace all occurrences in
|
||||
// faceEdges
|
||||
forAll(eFaces, i)
|
||||
{
|
||||
label triI = eFaces[i];
|
||||
const labelledTri& tri = surf[triI];
|
||||
|
||||
FixedList<label, 3>& fEdges = faceEdges[triI];
|
||||
forAll(tri, fp)
|
||||
{
|
||||
edge e(tri[fp], tri[tri.fcIndex(fp)]);
|
||||
|
||||
if (e == iter.key())
|
||||
{
|
||||
fEdges[fp] = edgeI;
|
||||
break;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// Add face to edgeFaces
|
||||
|
||||
edgeFace0[edgeI] = eFaces[0];
|
||||
|
||||
if (eFaces.size() >= 2)
|
||||
{
|
||||
edgeFace1[edgeI] = eFaces[1];
|
||||
|
||||
if (eFaces.size() > 2)
|
||||
{
|
||||
edgeFacesRest.insert
|
||||
(
|
||||
edgeI,
|
||||
SubList<label>(eFaces, eFaces.size()-2, 2)
|
||||
);
|
||||
}
|
||||
}
|
||||
|
||||
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);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
void Foam::isoSurface::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
|
||||
(
|
||||
"isoSurface::orientSurface(triSurface&, const label)"
|
||||
) << "problem" << abort(FatalError);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// 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++;
|
||||
}
|
||||
}
|
||||
|
||||
if (nOpen == 1 || nOpen == 2 || nOpen == 3)
|
||||
{
|
||||
return true;
|
||||
}
|
||||
else
|
||||
{
|
||||
return false;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
// 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;
|
||||
|
||||
// Remove any dangling triangles
|
||||
forAllConstIter(Map<labelList>, edgeFacesRest, iter)
|
||||
{
|
||||
// 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();
|
||||
|
||||
// Remove all dangling triangles
|
||||
if (danglingTriangle(faceEdges[edgeFace0[edgeI]], edgeFace1))
|
||||
{
|
||||
keepTriangles[edgeFace0[edgeI]] = false;
|
||||
nDangling++;
|
||||
}
|
||||
if (danglingTriangle(faceEdges[edgeFace1[edgeI]], edgeFace1))
|
||||
{
|
||||
keepTriangles[edgeFace1[edgeI]] = false;
|
||||
nDangling++;
|
||||
}
|
||||
forAll(otherEdgeFaces, i)
|
||||
{
|
||||
label triI = otherEdgeFaces[i];
|
||||
if (danglingTriangle(faceEdges[triI], edgeFace1))
|
||||
{
|
||||
keepTriangles[triI] = false;
|
||||
nDangling++;
|
||||
}
|
||||
}
|
||||
}
|
||||
return nDangling;
|
||||
}
|
||||
|
||||
|
||||
Foam::triSurface Foam::isoSurface::subsetMesh
|
||||
(
|
||||
const triSurface& s,
|
||||
@ -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,8 +1662,8 @@ Foam::isoSurface::isoSurface
|
||||
}
|
||||
|
||||
|
||||
|
||||
DynamicList<point> triPoints(nCutCells_);
|
||||
{
|
||||
DynamicList<point> triPoints(3*nCutCells_);
|
||||
DynamicList<label> triMeshCells(nCutCells_);
|
||||
|
||||
generateTriPoints
|
||||
@ -1973,8 +1678,8 @@ Foam::isoSurface::isoSurface
|
||||
snappedCc,
|
||||
snappedPoint,
|
||||
|
||||
triPoints,
|
||||
triMeshCells
|
||||
triPoints, // 3 points of the triangle
|
||||
triMeshCells // per triangle the originating cell
|
||||
);
|
||||
|
||||
if (debug)
|
||||
@ -1984,6 +1689,27 @@ Foam::isoSurface::isoSurface
|
||||
<< " unmerged points." << endl;
|
||||
}
|
||||
|
||||
label nOldPoints = triPoints.size();
|
||||
|
||||
// Trimmed to original triangle
|
||||
DynamicList<label> trimTriMap;
|
||||
// Trimmed to original point
|
||||
labelList trimTriPointMap;
|
||||
if (bounds_ != boundBox::greatBox)
|
||||
{
|
||||
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; // merged to unmerged triangle
|
||||
@ -2004,11 +1730,35 @@ Foam::isoSurface::isoSurface
|
||||
<< " 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]];
|
||||
}
|
||||
}
|
||||
|
||||
if (debug)
|
||||
{
|
||||
@ -2020,75 +1770,7 @@ Foam::isoSurface::isoSurface
|
||||
// Copied from surfaceCheck
|
||||
validTri(*this, triI);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (false)
|
||||
{
|
||||
List<FixedList<label, 3> > faceEdges;
|
||||
labelList edgeFace0, edgeFace1;
|
||||
Map<labelList> edgeFacesRest;
|
||||
|
||||
|
||||
while (true)
|
||||
{
|
||||
// Calculate addressing
|
||||
calcAddressing
|
||||
(
|
||||
*this,
|
||||
faceEdges,
|
||||
edgeFace0,
|
||||
edgeFace1,
|
||||
edgeFacesRest
|
||||
);
|
||||
|
||||
// See if any dangling triangles
|
||||
boolList keepTriangles;
|
||||
label nDangling = markDanglingTriangles
|
||||
(
|
||||
faceEdges,
|
||||
edgeFace0,
|
||||
edgeFace1,
|
||||
edgeFacesRest,
|
||||
keepTriangles
|
||||
);
|
||||
|
||||
if (debug)
|
||||
{
|
||||
Pout<< "isoSurface : detected " << nDangling
|
||||
<< " dangling triangles." << endl;
|
||||
}
|
||||
|
||||
if (nDangling == 0)
|
||||
{
|
||||
break;
|
||||
}
|
||||
|
||||
// Create face map (new to old)
|
||||
labelList subsetTriMap(findIndices(keepTriangles, true));
|
||||
|
||||
labelList subsetPointMap;
|
||||
labelList reversePointMap;
|
||||
triSurface::operator=
|
||||
(
|
||||
subsetMesh
|
||||
(
|
||||
*this,
|
||||
subsetTriMap,
|
||||
reversePointMap,
|
||||
subsetPointMap
|
||||
)
|
||||
);
|
||||
meshCells_ = labelField(meshCells_, subsetTriMap);
|
||||
inplaceRenumber(reversePointMap, triPointMergeMap_);
|
||||
}
|
||||
|
||||
orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
|
||||
}
|
||||
|
||||
|
||||
if (debug)
|
||||
{
|
||||
fileName stlFile = mesh_.time().path() + ".stl";
|
||||
Pout<< "Dumping surface to " << stlFile << endl;
|
||||
triSurface::write(stlFile);
|
||||
|
||||
@ -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,57 +343,38 @@ class isoSurface
|
||||
labelList& triMap // merged to unmerged triangle
|
||||
) 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
|
||||
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
|
||||
(
|
||||
const triSurface& s,
|
||||
@ -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>
|
||||
|
||||
@ -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,7 +1474,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
|
||||
}
|
||||
|
||||
|
||||
|
||||
{
|
||||
DynamicList<point> triPoints(nCutCells_);
|
||||
DynamicList<label> triMeshCells(nCutCells_);
|
||||
|
||||
@ -1633,8 +1500,32 @@ Foam::isoSurfaceCell::isoSurfaceCell
|
||||
<< " 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; // merged to unmerged triangle
|
||||
labelList triMap;
|
||||
triSurface::operator=
|
||||
(
|
||||
stitchTriPoints
|
||||
@ -1642,7 +1533,7 @@ Foam::isoSurfaceCell::isoSurfaceCell
|
||||
regularise, // check for duplicate tris
|
||||
triPoints,
|
||||
triPointMergeMap_, // unmerged to merged point
|
||||
triMap
|
||||
triMap // merged to unmerged triangle
|
||||
)
|
||||
);
|
||||
|
||||
@ -1652,11 +1543,35 @@ Foam::isoSurfaceCell::isoSurfaceCell
|
||||
<< " 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]];
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
if (debug)
|
||||
{
|
||||
@ -1730,8 +1645,6 @@ Foam::isoSurfaceCell::isoSurfaceCell
|
||||
meshCells_ = labelField(meshCells_, subsetTriMap);
|
||||
inplaceRenumber(reversePointMap, triPointMergeMap_);
|
||||
}
|
||||
|
||||
//orientSurface(*this, faceEdges, edgeFace0, edgeFace1, edgeFacesRest);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -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
|
||||
|
||||
@ -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)
|
||||
);
|
||||
if (triIndex == 0x0E)
|
||||
{
|
||||
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,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)
|
||||
);
|
||||
|
||||
if (triIndex == 0x0D)
|
||||
{
|
||||
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,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);
|
||||
Type p0p2 =
|
||||
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));
|
||||
pts.append(s02);
|
||||
pts.append(s13);
|
||||
pts.append(s13);
|
||||
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);
|
||||
// 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);
|
||||
Type p0p1 =
|
||||
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);
|
||||
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);
|
||||
Type p0p1 =
|
||||
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);
|
||||
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
|
||||
);
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -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,33 +278,47 @@ 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:
|
||||
case 0x0C:
|
||||
{
|
||||
Type tp1 =
|
||||
Type p0p2 =
|
||||
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2);
|
||||
Type tp2 =
|
||||
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(p1p3);
|
||||
points.append(p0p2);
|
||||
|
||||
points.append(p1p3);
|
||||
points.append
|
||||
(
|
||||
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;
|
||||
|
||||
case 0x0B:
|
||||
case 0x04:
|
||||
case 0x0B:
|
||||
{
|
||||
points.append
|
||||
(
|
||||
@ -311,56 +333,78 @@ void Foam::isoSurface::generateTriPoints
|
||||
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:
|
||||
case 0x0A:
|
||||
{
|
||||
Type tp0 =
|
||||
Type p0p1 =
|
||||
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
|
||||
Type tp1 =
|
||||
Type p2p3 =
|
||||
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
|
||||
|
||||
points.append(tp0);
|
||||
points.append(tp1);
|
||||
points.append(p0p1);
|
||||
points.append(p2p3);
|
||||
points.append
|
||||
(
|
||||
generatePoint(s0,p0,hasSnap0,snapP0,s3,p3,hasSnap3,snapP3)
|
||||
);
|
||||
points.append(tp0);
|
||||
|
||||
points.append(p0p1);
|
||||
points.append
|
||||
(
|
||||
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;
|
||||
|
||||
case 0x09:
|
||||
case 0x06:
|
||||
case 0x09:
|
||||
{
|
||||
Type tp0 =
|
||||
Type p0p1 =
|
||||
generatePoint(s0,p0,hasSnap0,snapP0,s1,p1,hasSnap1,snapP1);
|
||||
Type tp1 =
|
||||
Type p2p3 =
|
||||
generatePoint(s2,p2,hasSnap2,snapP2,s3,p3,hasSnap3,snapP3);
|
||||
|
||||
points.append(tp0);
|
||||
points.append(p0p1);
|
||||
points.append
|
||||
(
|
||||
generatePoint(s1,p1,hasSnap1,snapP1,s3,p3,hasSnap3,snapP3)
|
||||
);
|
||||
points.append(tp1);
|
||||
points.append(tp0);
|
||||
points.append(p2p3);
|
||||
|
||||
points.append(p0p1);
|
||||
points.append(p2p3);
|
||||
points.append
|
||||
(
|
||||
generatePoint(s0,p0,hasSnap0,snapP0,s2,p2,hasSnap2,snapP2)
|
||||
);
|
||||
points.append(tp1);
|
||||
}
|
||||
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;
|
||||
}
|
||||
|
||||
|
||||
|
||||
@ -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
|
||||
{
|
||||
@ -142,32 +176,22 @@ void Foam::sampledIsoSurface::getIsoFields() const
|
||||
if (debug)
|
||||
{
|
||||
Info<< "sampledIsoSurface::getIsoFields() :"
|
||||
<< "checking pointField " << pointFldName
|
||||
<< " for same time " << fvm.time().timeName()
|
||||
<< endl;
|
||||
<< " creating and storing pointField "
|
||||
<< pointFldName << " for time "
|
||||
<< fvm.time().timeName() << endl;
|
||||
}
|
||||
|
||||
if
|
||||
tmp<pointScalarField> tpfld
|
||||
(
|
||||
storedPointFieldPtr_.empty()
|
||||
|| (fvm.time().timeName() != storedPointFieldPtr_().instance())
|
||||
volPointInterpolation::New(fvm).interpolate
|
||||
(
|
||||
*volFieldPtr_,
|
||||
pointFldName,
|
||||
false
|
||||
)
|
||||
{
|
||||
if (debug)
|
||||
{
|
||||
Info<< "sampledIsoSurface::getIsoFields() :"
|
||||
<< " interpolating volField " << volFieldPtr_->name()
|
||||
<< " to get pointField " << pointFldName << endl;
|
||||
}
|
||||
|
||||
storedPointFieldPtr_.reset
|
||||
(
|
||||
volPointInterpolation::New(fvm)
|
||||
.interpolate(*volFieldPtr_).ptr()
|
||||
);
|
||||
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
|
||||
|
||||
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())
|
||||
|
||||
@ -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_;
|
||||
|
||||
|
||||
|
||||
@ -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),
|
||||
|
||||
@ -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_;
|
||||
|
||||
|
||||
@ -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)),
|
||||
|
||||
@ -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_;
|
||||
|
||||
|
||||
@ -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_;
|
||||
}
|
||||
|
||||
@ -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())
|
||||
{
|
||||
|
||||
@ -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);
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user