diff --git a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C index 2aee0df775..a4a4d0a8f0 100644 --- a/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C +++ b/applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshFeaturePoints.C @@ -99,11 +99,11 @@ bool Foam::conformalVoronoiMesh::meshableRegion { case extendedFeatureEdgeMesh::INSIDE: { - return (side == plane::FLIP) ? true : false; + return (side == plane::BACK); } case extendedFeatureEdgeMesh::OUTSIDE: { - return (side == plane::NORMAL) ? true : false; + return (side == plane::FRONT); } case extendedFeatureEdgeMesh::BOTH: { @@ -132,12 +132,12 @@ bool Foam::conformalVoronoiMesh::regionIsInside { plane::side sideA ( - ((masterPtVec & normalA) <= 0) ? plane::FLIP : plane::NORMAL + ((masterPtVec & normalA) <= 0) ? plane::BACK : plane::FRONT ); plane::side sideB ( - ((masterPtVec & normalB) <= 0) ? plane::FLIP : plane::NORMAL + ((masterPtVec & normalB) <= 0) ? plane::BACK : plane::FRONT ); const bool meshableRegionA = meshableRegion(sideA, volTypeA); diff --git a/applications/utilities/mesh/manipulation/mirrorMesh/mirrorMeshDict b/applications/utilities/mesh/manipulation/mirrorMesh/mirrorMeshDict index a8610eab1e..94c4f5a2cb 100644 --- a/applications/utilities/mesh/manipulation/mirrorMesh/mirrorMeshDict +++ b/applications/utilities/mesh/manipulation/mirrorMesh/mirrorMeshDict @@ -18,8 +18,8 @@ planeType pointAndNormal; pointAndNormalDict { - basePoint (0 0 0); - normalVector (0 1 0); + point (0 0 0); + normal (0 1 0); } planeTolerance 1e-3; diff --git a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C index 6d097114cc..fd2798fc0c 100644 --- a/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C +++ b/applications/utilities/surface/surfaceFeatureExtract/surfaceFeatureExtract.C @@ -493,7 +493,7 @@ int main(int argc, char *argv[]) Info<< "Only include feature edges that intersect the plane" << " with normal " << cutPlane.normal() - << " and base point " << cutPlane.refPoint() << endl; + << " and origin " << cutPlane.origin() << endl; features().subsetPlane(edgeStat, cutPlane); } diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.C b/src/OpenFOAM/meshes/boundBox/boundBox.C index 49124585b0..ad44f505e6 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBox.C +++ b/src/OpenFOAM/meshes/boundBox/boundBox.C @@ -45,12 +45,22 @@ const Foam::boundBox Foam::boundBox::invertedBox const Foam::faceList Foam::boundBox::faces ({ // Point and face order as per hex cellmodel - face{0, 4, 7, 3}, // x-min - face{1, 2, 6, 5}, // x-max - face{0, 1, 5, 4}, // y-min - face{3, 7, 6, 2}, // y-max - face{0, 3, 2, 1}, // z-min - face{4, 5, 6, 7} // z-max + face({0, 4, 7, 3}), // 0: x-min, left + face({1, 2, 6, 5}), // 1: x-max, right + face({0, 1, 5, 4}), // 2: y-min, bottom + face({3, 7, 6, 2}), // 3: y-max, top + face({0, 3, 2, 1}), // 4: z-min, back + face({4, 5, 6, 7}) // 5: z-max, front +}); + +const Foam::FixedList Foam::boundBox::faceNormals +({ + vector(-1, 0, 0), // 0: x-min, left + vector( 1, 0, 0), // 1: x-max, right + vector( 0, -1, 0), // 2: y-min, bottom + vector( 0, 1, 0), // 3: y-max, top + vector( 0, 0, -1), // 4: z-min, back + vector( 0, 0, 1) // 5: z-max, front }); @@ -107,8 +117,8 @@ Foam::boundBox::boundBox Foam::tmp Foam::boundBox::points() const { - tmp tpoints(new pointField(8)); - pointField& pt = tpoints.ref(); + auto tpt = tmp::New(8); + auto& pt = tpt.ref(); pt[0] = min_; // min-x, min-y, min-z pt[1] = point(max_.x(), min_.y(), min_.z()); // max-x, min-y, min-z @@ -119,7 +129,46 @@ Foam::tmp Foam::boundBox::points() const pt[6] = max_; // max-x, max-y, max-z pt[7] = point(min_.x(), max_.y(), max_.z()); // min-x, max-y, max-z - return tpoints; + return tpt; +} + + +Foam::tmp Foam::boundBox::faceCentres() const +{ + auto tpts = tmp::New(6); + auto& pts = tpts.ref(); + + forAll(pts, facei) + { + pts[facei] = faceCentre(facei); + } + + return tpts; +} + + +Foam::point Foam::boundBox::faceCentre(const direction facei) const +{ + point pt = boundBox::midpoint(); + + if (facei > 5) + { + FatalErrorInFunction + << "face should be [0..5]" + << abort(FatalError); + } + + switch (facei) + { + case 0: pt.x() = min().x(); break; // 0: x-min, left + case 1: pt.x() = max().x(); break; // 1: x-max, right + case 2: pt.y() = min().y(); break; // 2: y-min, bottom + case 3: pt.y() = max().y(); break; // 3: y-max, top + case 4: pt.z() = min().z(); break; // 4: z-min, back + case 5: pt.z() = max().z(); break; // 5: z-max, front + } + + return pt; } @@ -160,11 +209,11 @@ bool Foam::boundBox::intersects(const plane& pln) const bool below = false; tmp tpts(points()); - const pointField& pts = tpts(); + const auto& pts = tpts(); for (const point& p : pts) { - if (pln.sideOfPlane(p) == plane::NORMAL) + if (pln.sideOfPlane(p) == plane::FRONT) { above = true; } @@ -208,9 +257,9 @@ bool Foam::boundBox::contains return true; } - forAll(indices, i) + for (const label pointi : indices) { - if (!contains(points[indices[i]])) + if (!contains(points[pointi])) { return false; } @@ -250,9 +299,9 @@ bool Foam::boundBox::containsAny return true; } - forAll(indices, i) + for (const label pointi : indices) { - if (contains(points[indices[i]])) + if (contains(points[pointi])) { return true; } diff --git a/src/OpenFOAM/meshes/boundBox/boundBox.H b/src/OpenFOAM/meshes/boundBox/boundBox.H index abba73f881..5d0692dbdb 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBox.H +++ b/src/OpenFOAM/meshes/boundBox/boundBox.H @@ -78,6 +78,9 @@ public: //- Faces to point addressing, as per a 'hex' cell static const faceList faces; + //- The unit normal per face + static const FixedList faceNormals; + // Constructors @@ -177,6 +180,12 @@ public: //- Corner points in an order corresponding to a 'hex' cell tmp points() const; + //- Face midpoints + tmp faceCentres() const; + + //- Face centre of given face index + point faceCentre(const direction facei) const; + // Manipulate diff --git a/src/OpenFOAM/meshes/boundBox/boundBoxI.H b/src/OpenFOAM/meshes/boundBox/boundBoxI.H index daeb5d405b..0a117c16ca 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBoxI.H +++ b/src/OpenFOAM/meshes/boundBox/boundBoxI.H @@ -203,9 +203,9 @@ inline void Foam::boundBox::add { if (!points.empty()) { - forAll(indices, i) + for (const label pointi : indices) { - add(points[indices[i]]); + add(points[pointi]); } } } diff --git a/src/OpenFOAM/meshes/boundBox/boundBoxTemplates.C b/src/OpenFOAM/meshes/boundBox/boundBoxTemplates.C index 502d89ad5e..561a22fd25 100644 --- a/src/OpenFOAM/meshes/boundBox/boundBoxTemplates.C +++ b/src/OpenFOAM/meshes/boundBox/boundBoxTemplates.C @@ -74,9 +74,9 @@ void Foam::boundBox::add // points may be empty, but a FixedList is never empty if (!points.empty()) { - for (unsigned i=0; i < Size; ++i) + for (const label pointi : indices) { - add(points[indices[i]]); + add(points[pointi]); } } } @@ -95,9 +95,9 @@ bool Foam::boundBox::contains return false; } - forAll(indices, i) + for (const label pointi : indices) { - if (!contains(points[indices[i]])) + if (!contains(points[pointi])) { return false; } @@ -120,9 +120,9 @@ bool Foam::boundBox::containsAny return false; } - forAll(indices, i) + for (const label pointi : indices) { - if (contains(points[indices[i]])) + if (contains(points[pointi])) { return true; } diff --git a/src/OpenFOAM/meshes/meshShapes/face/face.H b/src/OpenFOAM/meshes/meshShapes/face/face.H index a07d5cb6c3..cedec1f73e 100644 --- a/src/OpenFOAM/meshes/meshShapes/face/face.H +++ b/src/OpenFOAM/meshes/meshShapes/face/face.H @@ -57,8 +57,7 @@ SourceFiles namespace Foam { -// Forward declaration of friend functions and operators - +// Forward declarations class face; class triFace; diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H index 8d46b372ac..c20edf59a0 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H +++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cut.H @@ -21,11 +21,11 @@ License You should have received a copy of the GNU General Public License along with OpenFOAM. If not, see . -Description: - Functions which cut triangles and tetrahedra. Generic operations are - applied to each half of a cut. +Description + Functions for cutting triangles and tetrahedra. + Generic operations are applied to each half of a cut. -SourceFiles: +SourceFiles cutI.H cutTemplates.C @@ -478,9 +478,9 @@ public: // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // -//- Cut a triangle along the zero plane defined by the given levels. Templated -// on aboveOp and belowOp; the operations applied to the regions on either side -// of the cut. +//- Cut a triangle along the zero plane defined by the given levels. +// Templated on aboveOp and belowOp; the operations applied to the regions +// on either side of the cut. template typename cut::opAddResult::type triCut ( @@ -495,7 +495,7 @@ template typename cut::opAddResult::type triCut ( const FixedList& tri, - const plane& s, + const plane& pln, const AboveOp& aboveOp, const BelowOp& belowOp ); @@ -515,7 +515,7 @@ template typename cut::opAddResult::type tetCut ( const FixedList& tet, - const plane& s, + const plane& pln, const AboveOp& aboveOp, const BelowOp& belowOp ); diff --git a/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C b/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C index 1ec2438090..187615bad8 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C +++ b/src/OpenFOAM/meshes/primitiveShapes/cut/cutTemplates.C @@ -51,7 +51,7 @@ typename Foam::cut::opAddResult::type Foam::triCut // opposite the first vertex. This may change the sign of the tri. FixedList indices({0, 1, 2}); label i; - for (i = 0; i < 3; ++ i) + for (i = 0; i < 3; ++i) { if (level[(i + 1)%3]*level[(i + 2)%3] >= 0) { @@ -80,7 +80,7 @@ typename Foam::cut::opAddResult::type Foam::triCut // Slice off one corner to form a tri and a quad Pair f; - for (label i = 0; i < 2; ++ i) + for (label i = 0; i < 2; ++i) { f[i] = l[0]/(l[0] - l[i+1]); } @@ -99,7 +99,7 @@ template typename Foam::cut::opAddResult::type Foam::triCut ( const FixedList& tri, - const plane& p, + const plane& pln, const AboveOp& aboveOp, const BelowOp& belowOp ) @@ -108,8 +108,7 @@ typename Foam::cut::opAddResult::type Foam::triCut FixedList level; for (label i = 0; i < 3; ++i) { - // uses unit-normal - level[i] = (tri[i] - p.refPoint()) & p.normal(); + level[i] = pln.signedDistance(tri[i]); } // Run the level set function @@ -129,7 +128,7 @@ typename Foam::cut::opAddResult::type Foam::tetCut // Get the min and max over all four vertices and quick return if there is // no change of sign scalar levelMin = VGREAT, levelMax = - VGREAT; - for (label i = 0; i < 4; ++ i) + for (label i = 0; i < 4; ++i) { levelMin = min(levelMin, level[i]); levelMax = max(levelMax, level[i]); @@ -176,7 +175,7 @@ typename Foam::cut::opAddResult::type Foam::tetCut if (n > 2) { n = 4 - n; - for (label i = 0; i < 2; ++ i) + for (label i = 0; i < 2; ++i) { Swap(indices[i], indices[3-i]); } @@ -199,7 +198,7 @@ typename Foam::cut::opAddResult::type Foam::tetCut { // Slice off one corner to form a tet and a prism FixedList f; - for (label i = 0; i < 3; ++ i) + for (label i = 0; i < 3; ++i) { f[i] = l[0]/(l[0] - l[i+1]); } @@ -216,9 +215,9 @@ typename Foam::cut::opAddResult::type Foam::tetCut { // Slice off two corners to form two prisms FixedList f; - for (label i = 0; i < 2; ++ i) + for (label i = 0; i < 2; ++i) { - for (label j = 0; j < 2; ++ j) + for (label j = 0; j < 2; ++j) { f[2*i+j] = l[i]/(l[i] - l[j+2]); } @@ -245,7 +244,7 @@ template typename Foam::cut::opAddResult::type Foam::tetCut ( const FixedList& tet, - const plane& p, + const plane& pln, const AboveOp& aboveOp, const BelowOp& belowOp ) @@ -254,8 +253,7 @@ typename Foam::cut::opAddResult::type Foam::tetCut FixedList level; for (label i = 0; i < 4; ++i) { - // uses unit-normal - level[i] = (tet[i] - p.refPoint()) & p.normal(); + level[i] = pln.signedDistance(tet[i]); } // Run the level set function diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C index 348802c69f..ec1cd8f19c 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C +++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.C @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2016-2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2016-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -28,19 +28,48 @@ License // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -void Foam::plane::calcPntAndVec(const scalarList& C) +void Foam::plane::makeUnitNormal +( + const char * const caller, + const bool notTest +) { - if (mag(C[0]) > VSMALL) + const scalar magNormal(Foam::mag(normal_)); + + if (magNormal < VSMALL) { - point_ = vector((-C[3]/C[0]), 0, 0); + FatalErrorInFunction + << "Plane normal has zero length.\nCalled from " << caller + << abort(FatalError); } - else if (mag(C[1]) > VSMALL) + + if (notTest) { - point_ = vector(0, (-C[3]/C[1]), 0); + normal_ /= magNormal; } - else if (mag(C[2]) > VSMALL) +} + + +void Foam::plane::calcFromCoeffs +( + const scalar a, + const scalar b, + const scalar c, + const scalar d, + const char * const caller +) +{ + if (mag(a) > VSMALL) { - point_ = vector(0, 0, (-C[3]/C[2])); + origin_ = vector((-d/a), 0, 0); + } + else if (mag(b) > VSMALL) + { + origin_ = vector(0, (-d/b), 0); + } + else if (mag(c) > VSMALL) + { + origin_ = vector(0, 0, (-d/c)); } else { @@ -49,30 +78,22 @@ void Foam::plane::calcPntAndVec(const scalarList& C) << abort(FatalError); } - normal_ = vector(C[0], C[1], C[2]); - scalar magUnitVector(mag(normal_)); - - if (magUnitVector < VSMALL) - { - FatalErrorInFunction - << "Plane normal defined with zero length" - << abort(FatalError); - } - - normal_ /= magUnitVector; + normal_ = vector(a, b, c); + makeUnitNormal(caller); } -void Foam::plane::calcPntAndVec +void Foam::plane::calcFromEmbeddedPoints ( const point& point1, const point& point2, - const point& point3 + const point& point3, + const char * const caller ) { - point_ = (point1 + point2 + point3)/3; - vector line12 = point1 - point2; - vector line23 = point2 - point3; + origin_ = (point1 + point2 + point3)/3; + const vector line12 = point1 - point2; + const vector line23 = point2 - point3; if ( @@ -87,141 +108,115 @@ void Foam::plane::calcPntAndVec } normal_ = line12 ^ line23; - scalar magUnitVector(mag(normal_)); - if (magUnitVector < VSMALL) - { - FatalErrorInFunction - << "Plane normal defined with zero length" << nl - << "Bad points:" << point1 << ' ' << point2 << ' ' << point3 - << abort(FatalError); - } - - normal_ /= magUnitVector; + makeUnitNormal(caller); } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::plane::plane() -{} - - Foam::plane::plane(const vector& normalVector) : normal_(normalVector), - point_(Zero) + origin_(Zero) { - scalar magUnitVector(mag(normal_)); - - if (magUnitVector > VSMALL) - { - normal_ /= magUnitVector; - } - else - { - FatalErrorInFunction - << "plane normal has zero length. base point:" << point_ - << abort(FatalError); - } + makeUnitNormal(FUNCTION_NAME); } Foam::plane::plane ( - const point& basePoint, + const point& originPoint, const vector& normalVector, - const bool normalise + const bool doNormalise ) : normal_(normalVector), - point_(basePoint) + origin_(originPoint) { - scalar magSqrNormalVector(magSqr(normal_)); - - if (magSqrNormalVector < VSMALL) - { - FatalErrorInFunction - << "plane normal has zero length. base point:" << point_ - << abort(FatalError); - } - - if (normalise) - { - normal_ /= Foam::sqrt(magSqrNormalVector); - } + makeUnitNormal(FUNCTION_NAME, doNormalise); } -Foam::plane::plane(const scalarList& C) +Foam::plane::plane(const scalarList& coeffs) { - calcPntAndVec(C); + calcFromCoeffs + ( + coeffs[0], + coeffs[1], + coeffs[2], + coeffs[3], + FUNCTION_NAME + ); } -Foam::plane::plane -( - const point& a, - const point& b, - const point& c -) +Foam::plane::plane(const FixedList& coeffs) { - calcPntAndVec(a, b, c); + calcFromCoeffs + ( + coeffs[0], + coeffs[1], + coeffs[2], + coeffs[3], + FUNCTION_NAME + ); +} + + +Foam::plane::plane(const point& a, const point& b, const point& c) +{ + calcFromEmbeddedPoints(a, b, c, FUNCTION_NAME); } Foam::plane::plane(const dictionary& dict) : normal_(Zero), - point_(Zero) + origin_(Zero) { - const word planeType(dict.lookup("planeType")); + const word planeType(dict.get("planeType")); if (planeType == "planeEquation") { const dictionary& subDict = dict.subDict("planeEquationDict"); - scalarList C(4); - - C[0] = readScalar(subDict.lookup("a")); - C[1] = readScalar(subDict.lookup("b")); - C[2] = readScalar(subDict.lookup("c")); - C[3] = readScalar(subDict.lookup("d")); - - calcPntAndVec(C); + calcFromCoeffs + ( + subDict.get("a"), + subDict.get("b"), + subDict.get("c"), + subDict.get("d"), + "planeEquationDict" // caller name for makeUnitNormal + ); } else if (planeType == "embeddedPoints") { const dictionary& subDict = dict.subDict("embeddedPointsDict"); - point point1(subDict.lookup("point1")); - point point2(subDict.lookup("point2")); - point point3(subDict.lookup("point3")); + calcFromEmbeddedPoints + ( + subDict.get("point1"), + subDict.get("point2"), + subDict.get("point3"), + "embeddedPointsDict" // caller name for makeUnitNormal + ); - calcPntAndVec(point1, point2, point3); } else if (planeType == "pointAndNormal") { const dictionary& subDict = dict.subDict("pointAndNormalDict"); - point_ = - subDict.found("basePoint") - ? subDict.lookup("basePoint") - : subDict.lookup("point"); + origin_ = subDict.getCompat("point", {{"basePoint", 1612}}); + normal_ = subDict.getCompat("normal", {{"normalVector", 1612}}); - normal_ = - subDict.found("normalVector") - ? subDict.lookup("normalVector") - : subDict.lookup("normal"); - - normal_ /= mag(normal_); + makeUnitNormal("pointAndNormalDict"); } else { FatalIOErrorInFunction(dict) << "Invalid plane type: " << planeType << nl - << "Valid options include: planeEquation, embeddedPoints and " - << "pointAndNormal" + << "Valid options: (planeEquation embeddedPoints pointAndNormal)" << abort(FatalIOError); } } @@ -230,93 +225,61 @@ Foam::plane::plane(const dictionary& dict) Foam::plane::plane(Istream& is) : normal_(is), - point_(is) + origin_(is) { - scalar magUnitVector(mag(normal_)); - - if (magUnitVector > VSMALL) - { - normal_ /= magUnitVector; - } - else - { - FatalErrorInFunction - << "plane normal has zero length. base point:" << point_ - << abort(FatalError); - } + makeUnitNormal(FUNCTION_NAME); } // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -const Foam::vector& Foam::plane::normal() const -{ - return normal_; -} - - -const Foam::point& Foam::plane::refPoint() const -{ - return point_; -} - - Foam::FixedList Foam::plane::planeCoeffs() const { - FixedList C(4); + FixedList coeffs(4); - scalar magX = mag(normal_.x()); - scalar magY = mag(normal_.y()); - scalar magZ = mag(normal_.z()); + const scalar magX = mag(normal_.x()); + const scalar magY = mag(normal_.y()); + const scalar magZ = mag(normal_.z()); if (magX > magY) { if (magX > magZ) { - C[0] = 1; - C[1] = normal_.y()/normal_.x(); - C[2] = normal_.z()/normal_.x(); + coeffs[0] = 1; + coeffs[1] = normal_.y()/normal_.x(); + coeffs[2] = normal_.z()/normal_.x(); } else { - C[0] = normal_.x()/normal_.z(); - C[1] = normal_.y()/normal_.z(); - C[2] = 1; + coeffs[0] = normal_.x()/normal_.z(); + coeffs[1] = normal_.y()/normal_.z(); + coeffs[2] = 1; } } else { if (magY > magZ) { - C[0] = normal_.x()/normal_.y(); - C[1] = 1; - C[2] = normal_.z()/normal_.y(); + coeffs[0] = normal_.x()/normal_.y(); + coeffs[1] = 1; + coeffs[2] = normal_.z()/normal_.y(); } else { - C[0] = normal_.x()/normal_.z(); - C[1] = normal_.y()/normal_.z(); - C[2] = 1; + coeffs[0] = normal_.x()/normal_.z(); + coeffs[1] = normal_.y()/normal_.z(); + coeffs[2] = 1; } } - C[3] = - C[0] * point_.x() - - C[1] * point_.y() - - C[2] * point_.z(); + coeffs[3] = + ( + - coeffs[0] * origin_.x() + - coeffs[1] * origin_.y() + - coeffs[2] * origin_.z() + ); - return C; -} - - -Foam::point Foam::plane::nearestPoint(const point& p) const -{ - return p - normal_*((p - point_) & normal_); -} - - -Foam::scalar Foam::plane::distance(const point& p) const -{ - return mag((p - point_) & normal_); + return coeffs; } @@ -326,9 +289,9 @@ Foam::scalar Foam::plane::normalIntersect const vector& dir ) const { - scalar denom = stabilise((dir & normal_), VSMALL); + const scalar denom = stabilise((dir & normal_), VSMALL); - return ((point_ - pnt0) & normal_)/denom; + return ((origin_ - pnt0) & normal_)/denom; } @@ -339,22 +302,22 @@ Foam::plane::ray Foam::plane::planeIntersect(const plane& plane2) const // for that (now 2x2 equation in x and y) // Better: use either z=0 or x=0 or y=0. - const vector& n1 = normal(); + const vector& n1 = this->normal(); const vector& n2 = plane2.normal(); - const point& p1 = refPoint(); - const point& p2 = plane2.refPoint(); + const point& p1 = this->origin(); + const point& p2 = plane2.origin(); - scalar n1p1 = n1&p1; - scalar n2p2 = n2&p2; + const scalar n1p1 = n1 & p1; + const scalar n2p2 = n2 & p2; - vector dir = n1 ^ n2; + const vector dir = n1 ^ n2; // Determine zeroed out direction (can be x,y or z) by looking at which // has the largest component in dir. - scalar magX = mag(dir.x()); - scalar magY = mag(dir.y()); - scalar magZ = mag(dir.z()); + const scalar magX = mag(dir.x()); + const scalar magY = mag(dir.y()); + const scalar magZ = mag(dir.z()); direction iZero, i1, i2; @@ -389,6 +352,7 @@ Foam::plane::ray Foam::plane::planeIntersect(const plane& plane2) const } } + vector pt; pt[iZero] = 0; @@ -428,7 +392,7 @@ Foam::point Foam::plane::somePointInPlane(const scalar dist) const const FixedList coeff(planeCoeffs()); // Perturb the base-point - point p = refPoint() + point::uniform(dist); + point p = origin_ + point::uniform(dist); if (Foam::mag(coeff[2]) < SMALL) { @@ -450,14 +414,6 @@ Foam::point Foam::plane::somePointInPlane(const scalar dist) const } -Foam::plane::side Foam::plane::sideOfPlane(const point& p) const -{ - const scalar angle((p - point_) & normal_); - - return (angle < 0 ? FLIP : NORMAL); -} - - Foam::point Foam::plane::mirror(const point& p) const { const vector mirroredPtDir = p - nearestPoint(p); @@ -479,32 +435,18 @@ void Foam::plane::writeDict(Ostream& os) const os.beginBlock("pointAndNormalDict"); - os.writeEntry("point", point_); + os.writeEntry("point", origin_); os.writeEntry("normal", normal_); os.endBlock() << flush; } -// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // - -bool Foam::operator==(const plane& a, const plane& b) -{ - return (a.refPoint() == b.refPoint() && a.normal() == b.normal()); -} - - -bool Foam::operator!=(const plane& a, const plane& b) -{ - return !(a == b); -} - // * * * * * * * * * * * * * * * Friend Functions * * * * * * * * * * * * * // Foam::Ostream& Foam::operator<<(Ostream& os, const plane& pln) { - os << pln.normal_ << token::SPACE << pln.point_; - + os << pln.normal() << token::SPACE << pln.origin(); return os; } diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H index 88a39b3162..428e9aabad 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H +++ b/src/OpenFOAM/meshes/primitiveShapes/plane/plane.H @@ -3,7 +3,7 @@ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox \\ / O peration | \\ / A nd | Copyright (C) 2011-2016 OpenFOAM Foundation - \\/ M anipulation | Copyright (C) 2017 OpenCFD Ltd. + \\/ M anipulation | Copyright (C) 2017-2018 OpenCFD Ltd. ------------------------------------------------------------------------------- License This file is part of OpenFOAM. @@ -80,8 +80,7 @@ SourceFiles namespace Foam { -// Forward declaration of friend functions and operators - +// Forward Declarations class plane; Ostream& operator<<(Ostream& os, const plane& pln); @@ -97,28 +96,30 @@ public: //- Side of the plane enum side { - NORMAL = 1, //!< The front (or normal) side of the plane - FLIP = -1 //!< The back (or flip) side of the plane + FRONT = 1, //!< The front (positive normal) side of the plane + BACK = -1, //!< The back (negative normal) side of the plane + NORMAL = 1, //!< Same as FRONT + FLIP = -1 //!< Same as BACK }; - //- A direction and a reference point + //- A reference point and direction class ray { - point refPoint_; + point pt_; vector dir_; public: - ray(const point& refPoint, const vector& dir) + ray(const point& pt, const vector& dir) : - refPoint_(refPoint), + pt_(pt), dir_(dir) {} const point& refPoint() const { - return refPoint_; + return pt_; } const vector& dir() const @@ -130,27 +131,43 @@ public: private: - // Private data + // Private Data - //- Normal + //- The unit normal of the plane vector normal_; - //- Reference point - point point_; + //- The origin or reference point for the plane + point origin_; // Private Member Functions - //- Calculates basePoint and normal vector given plane coefficients - void calcPntAndVec(const scalarList& C); + //- Normalise normal_ and emit error if its mag is less than VSMALL + // Optionally pass as test only, without normalisation. + void makeUnitNormal + ( + const char * const caller, + const bool notTest = true + ); - //- Calculates basePoint and normal vector given three points + //- Calculates point and normal given plane coefficients. + void calcFromCoeffs + ( + const scalar a, + const scalar b, + const scalar c, + const scalar d, + const char * const caller + ); + + //- Calculates point and normal vector given three points //- Normal vector determined using right hand rule - void calcPntAndVec + void calcFromEmbeddedPoints ( const point& point1, const point& point2, - const point& point3 + const point& point3, + const char * const caller ); @@ -159,7 +176,7 @@ public: // Constructors //- Construct null. Does not set normal and point. - plane(); + inline plane(); //- Construct from normal vector through the origin. // The vector is normalised to a unit vector on input. @@ -167,11 +184,11 @@ public: //- Construct from normal vector and point in plane. // By default, the vector is normalised to a unit vector on input. - explicit plane + plane ( - const point& basePoint, + const point& originPoint, const vector& normalVector, - const bool normalise = true + const bool doNormalise = true ); //- Construct from three points in plane @@ -182,37 +199,50 @@ public: const point& point3 ); - //- Construct from coefficients for the - // plane equation: ax + by + cz + d = 0 - explicit plane(const scalarList& C); + //- Construct from coefficients for the plane equation: + //- ax + by + cz + d = 0 + explicit plane(const scalarList& coeffs); + + //- Construct from coefficients for the plane equation: + //- ax + by + cz + d = 0 + explicit plane(const FixedList& coeffs); //- Construct from dictionary explicit plane(const dictionary& dict); - //- Construct from Istream. Assumes the base + normal notation. + //- Construct from Istream. Assumes (normal) (origin) input. explicit plane(Istream& is); // Member Functions - //- Return the plane unit normal - const vector& normal() const; + //- The plane unit normal + inline const vector& normal() const; - //- Return plane base point - const point& refPoint() const; + //- The plane base point (same as refPoint) + inline const point& origin() const; - //- Return coefficients for the - // plane equation: ax + by + cz + d = 0 + //- The plane base point (same as origin) + inline const point& refPoint() const; + + //- Flip the plane by reversing the normal + inline void flip(); + + //- Return coefficients for the plane equation: + //- ax + by + cz + d = 0 FixedList planeCoeffs() const; //- Return nearest point in the plane for the given point - point nearestPoint(const point& p) const; + inline point nearestPoint(const point& p) const; + + //- Return distance (magnitude) from the given point to the plane + inline scalar distance(const point& p) const; //- Return distance from the given point to the plane - scalar distance(const point& p) const; + inline scalar signedDistance(const point& p) const; //- Return cut coefficient for plane and line defined by - // origin and direction + //- origin and direction scalar normalIntersect(const point& pnt0, const vector& dir) const; //- Return cut coefficient for plane and ray @@ -246,7 +276,15 @@ public: //- Return the side of the plane that the point is on. // If the point is on the plane, then returns NORMAL. - side sideOfPlane(const point& p) const; + inline side sideOfPlane(const point& p) const; + + //- The sign for the side of the plane that the point is on. + // Uses the supplied tolerance for rounding around zero. + // \return + // - 0: on plane + // - +1: front-side + // - -1: back-side + inline int sign(const point& p, const scalar tol = SMALL) const; //- Mirror the supplied point in the plane. Return the mirrored point. point mirror(const point& p) const; @@ -265,14 +303,16 @@ public: // Global Operators -bool operator==(const plane& a, const plane& b); -bool operator!=(const plane& a, const plane& b); +inline bool operator==(const plane& a, const plane& b); +inline bool operator!=(const plane& a, const plane& b); // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace Foam +#include "planeI.H" + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // #endif diff --git a/src/OpenFOAM/meshes/primitiveShapes/plane/planeI.H b/src/OpenFOAM/meshes/primitiveShapes/plane/planeI.H new file mode 100644 index 0000000000..83d14f529e --- /dev/null +++ b/src/OpenFOAM/meshes/primitiveShapes/plane/planeI.H @@ -0,0 +1,106 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2018 OpenCFD Ltd. + \\/ 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 . + +\*---------------------------------------------------------------------------*/ + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +inline Foam::plane::plane() +{} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +inline const Foam::vector& Foam::plane::normal() const +{ + return normal_; +} + + +inline const Foam::point& Foam::plane::origin() const +{ + return origin_; +} + + +inline const Foam::point& Foam::plane::refPoint() const +{ + return origin_; +} + + +inline void Foam::plane::flip() +{ + normal_ = -normal_; +} + + +inline Foam::point Foam::plane::nearestPoint(const point& p) const +{ + return p - normal_*((p - origin_) & normal_); +} + + +inline Foam::scalar Foam::plane::distance(const point& p) const +{ + return mag(signedDistance(p)); +} + + +inline Foam::scalar Foam::plane::signedDistance(const point& p) const +{ + return ((p - origin_) & normal_); +} + + +inline Foam::plane::side Foam::plane::sideOfPlane(const point& p) const +{ + const scalar dist = signedDistance(p); + + return (dist < 0 ? BACK : FRONT); +} + + +inline int Foam::plane::sign(const point& p, const scalar tol) const +{ + const scalar dist = signedDistance(p); + + return ((dist < -tol) ? -1 : (dist > tol) ? +1 : 0); +} + + +// * * * * * * * * * * * * * * * Global Operators * * * * * * * * * * * * * // + +inline bool Foam::operator==(const plane& a, const plane& b) +{ + return (a.origin() == b.origin() && a.normal() == b.normal()); +} + + +inline bool Foam::operator!=(const plane& a, const plane& b) +{ + return !(a == b); +} + + +// ************************************************************************* // diff --git a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H index 64e18c3dd7..5fae5cbd91 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H +++ b/src/OpenFOAM/meshes/primitiveShapes/tetrahedron/tetrahedronI.H @@ -564,7 +564,7 @@ template template inline void Foam::tetrahedron::tetSliceWithPlane ( - const plane& pl, + const plane& pln, const tetPoints& tet, AboveTetOp& aboveOp, BelowTetOp& belowOp @@ -575,8 +575,7 @@ inline void Foam::tetrahedron::tetSliceWithPlane label nPos = 0; forAll(tet, i) { - // with plane unit-normal - d[i] = ((tet[i] - pl.refPoint()) & pl.normal()); + d[i] = pln.signedDistance(tet[i]); if (d[i] > 0) { ++nPos; @@ -589,8 +588,8 @@ inline void Foam::tetrahedron::tetSliceWithPlane } else if (nPos == 3) { - // Sliced into below tet and above prism. Prism gets split into - // two tets. + // Sliced into below tet and above prism. + // Prism gets split into two tets. // Find the below tet label i0 = -1; diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C index 8085ccd8f0..68f7d04c33 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.C @@ -33,7 +33,7 @@ template template inline void Foam::triangle::triSliceWithPlane ( - const plane& pl, + const plane& pln, const triPoints& tri, AboveOp& aboveOp, BelowOp& belowOp @@ -48,7 +48,7 @@ inline void Foam::triangle::triSliceWithPlane label negI = -1; forAll(tri, i) { - d[i] = ((tri[i] - pl.refPoint()) & pl.normal()); + d[i] = pln.signedDistance(tri[i]); if (d[i] > 0) { diff --git a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H index 5f09ce6026..d8de8e9609 100644 --- a/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H +++ b/src/OpenFOAM/meshes/primitiveShapes/triangle/triangle.H @@ -50,13 +50,13 @@ SourceFiles namespace Foam { +// Forward declarations + class Istream; class Ostream; class triPoints; class plane; -// Forward declaration of friend functions and operators - template class triangle; template @@ -158,7 +158,7 @@ private: template inline static void triSliceWithPlane ( - const plane& pl, + const plane& pln, const triPoints& tri, AboveOp& aboveOp, BelowOp& belowOp @@ -327,7 +327,7 @@ public: template inline void sliceWithPlane ( - const plane& pl, + const plane& pln, AboveOp& aboveOp, BelowOp& belowOp ) const; diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C index da69df3e2a..ffc6856d80 100644 --- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C +++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.C @@ -29,12 +29,12 @@ License const Foam::faceList Foam::treeBoundBox::faces ({ - face{0, 4, 6, 2}, // left - face{1, 3, 7, 5}, // right - face{0, 1, 5, 4}, // bottom - face{2, 6, 7, 3}, // top - face{0, 2, 3, 1}, // back - face{4, 5, 7, 6} // front + face({0, 4, 6, 2}), // 0: x-min, left + face({1, 3, 7, 5}), // 1: x-max, right + face({0, 1, 5, 4}), // 2: y-min, bottom + face({2, 6, 7, 3}), // 3: y-max, top + face({0, 2, 3, 1}), // 4: z-min, back + face({4, 5, 7, 6}) // 5: z-max, front }); const Foam::edgeList Foam::treeBoundBox::edges @@ -53,16 +53,6 @@ const Foam::edgeList Foam::treeBoundBox::edges {2, 6} }); -const Foam::FixedList Foam::treeBoundBox::faceNormals -({ - vector(-1, 0, 0), // left - vector( 1, 0, 0), // right - vector( 0, -1, 0), // bottom - vector( 0, 1, 0), // top - vector( 0, 0, -1), // back - vector( 0, 0, 1) // front -}); - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // @@ -104,15 +94,15 @@ Foam::treeBoundBox::treeBoundBox Foam::tmp Foam::treeBoundBox::points() const { - tmp tpoints(new pointField(8)); - pointField& pts = tpoints.ref(); + auto tpts = tmp::New(8); + auto& pts = tpts.ref(); forAll(pts, octant) { pts[octant] = corner(octant); } - return tpoints; + return tpts; } diff --git a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H index 7177cb020f..c8248c29e0 100644 --- a/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H +++ b/src/OpenFOAM/meshes/treeBoundBox/treeBoundBox.H @@ -71,11 +71,8 @@ SourceFiles namespace Foam { +// Forward declarations class Random; - - -// Forward declaration of friend functions and operators - class treeBoundBox; Istream& operator>>(Istream& is, treeBoundBox& bb); Ostream& operator<<(Ostream& os, const treeBoundBox& bb); @@ -98,20 +95,21 @@ public: // Every octant/corner point is the combination of three faces. enum octantBit { - RIGHTHALF = 0x1 << 0, - TOPHALF = 0x1 << 1, - FRONTHALF = 0x1 << 2 + RIGHTHALF = 0x1, + TOPHALF = 0x2, + FRONTHALF = 0x4 }; - //- Face codes + //- Face codes. Identical order and meaning as per hex cellmodel + //- and boundBox, but have a different point ordering. enum faceId { - LEFT = 0, - RIGHT = 1, - BOTTOM = 2, - TOP = 3, - BACK = 4, - FRONT = 5 + LEFT = 0, //!< 0: x-min, left + RIGHT = 1, //!< 1: x-max, right + BOTTOM = 2, //!< 2: y-min, bottom + TOP = 3, //!< 3: y-max, top + BACK = 4, //!< 4: z-min, back + FRONT = 5 //!< 5: z-max, front }; //- Bits used for face coding @@ -152,9 +150,6 @@ public: //- Edge to point addressing static const edgeList edges; - //- The unit normal per face - static const FixedList faceNormals; - // Constructors diff --git a/src/OpenFOAM/primitives/Scalar/Scalar.H b/src/OpenFOAM/primitives/Scalar/Scalar.H index f091f42765..ea2564075d 100644 --- a/src/OpenFOAM/primitives/Scalar/Scalar.H +++ b/src/OpenFOAM/primitives/Scalar/Scalar.H @@ -378,15 +378,15 @@ inline Scalar sqrtSumSqr(const Scalar a, const Scalar b) //- Stabilisation around zero for division -inline Scalar stabilise(const Scalar s, const Scalar small) +inline Scalar stabilise(const Scalar s, const Scalar tol) { if (s >= 0) { - return s + small; + return s + tol; } else { - return s - small; + return s - tol; } } diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C index 20104f5b6a..97e8b1a00e 100644 --- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C +++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.C @@ -44,7 +44,7 @@ Foam::scalar Foam::faceAreaIntersect::tol = 1e-6; void Foam::faceAreaIntersect::triSliceWithPlane ( const triPoints& tri, - const plane& p, + const plane& pln, FixedList& tris, label& nTris, const scalar len @@ -61,7 +61,7 @@ void Foam::faceAreaIntersect::triSliceWithPlane label copI = -1; forAll(tri, i) { - d[i] = ((tri[i] - p.refPoint()) & p.normal()); + d[i] = pln.signedDistance(tri[i]); if (mag(d[i]) < tol*len) { diff --git a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H index 53394b962f..8119c0abda 100644 --- a/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H +++ b/src/meshTools/AMIInterpolation/faceAreaIntersect/faceAreaIntersect.H @@ -133,7 +133,7 @@ private: void triSliceWithPlane ( const triPoints& tri, - const plane& p, + const plane& pln, FixedList& tris, label& nTris, const scalar len diff --git a/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.C b/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.C index d75dff8ede..52a38d674a 100644 --- a/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.C +++ b/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.C @@ -57,26 +57,21 @@ Foam::pointIndexHit Foam::searchableDisk::findNearest { pointIndexHit info(false, sample, -1); - vector v(sample - origin_); + vector v(sample - origin()); // Decompose sample-origin into normal and parallel component - const scalar parallel = (v & normal_); + const scalar parallel = (v & normal()); // Remove the parallel component and normalise - v -= parallel*normal_; - scalar magV = mag(v); + v -= parallel * normal(); + + const scalar magV = mag(v); + + v.normalise(); - if (magV < ROOTVSMALL) - { - v = Zero; - } - else - { - v /= magV; - } // Clip to radius. - info.setPoint(origin_ + min(magV, radius_)*v); + info.setPoint(origin() + min(magV, radius_)*v); if (magSqr(sample - info.rawPoint()) < nearestDistSqr) { @@ -97,31 +92,26 @@ void Foam::searchableDisk::findLine { info = pointIndexHit(false, Zero, -1); - vector v(start - origin_); + vector v(start - origin()); // Decompose sample-origin into normal and parallel component - const scalar parallel = (v & normal_); + const scalar parallel = (v & normal()); - if (sign(parallel) == sign((end - origin_) & normal_)) + if (Foam::sign(parallel) == Foam::sign(plane::signedDistance(end))) { return; } // Remove the parallel component and normalise - v -= parallel*normal_; - scalar magV = mag(v); + v -= parallel * normal(); + + const scalar magV = mag(v); + + v.normalise(); - if (magV < ROOTVSMALL) - { - v = Zero; - } - else - { - v /= magV; - } // Set (hit or miss) to intersection of ray and plane of disk - info.setPoint(origin_ + magV*v); + info.setPoint(origin() + magV*v); if (magV <= radius_) { @@ -136,30 +126,29 @@ void Foam::searchableDisk::findLine Foam::searchableDisk::searchableDisk ( const IOobject& io, - const point& origin, - const point& normal, + const point& originPoint, + const vector& normalVector, const scalar radius ) : searchableSurface(io), - origin_(origin), - normal_(normal/mag(normal)), + plane(originPoint, normalVector), radius_(radius) { // Rough approximation of bounding box - //vector span(radius_, radius_, radius_); + // vector span(radius_, radius_, radius_); // See searchableCylinder vector span ( - sqrt(sqr(normal_.y()) + sqr(normal_.z())), - sqrt(sqr(normal_.x()) + sqr(normal_.z())), - sqrt(sqr(normal_.x()) + sqr(normal_.y())) + sqrt(sqr(normal().y()) + sqr(normal().z())), + sqrt(sqr(normal().x()) + sqr(normal().z())), + sqrt(sqr(normal().x()) + sqr(normal().y())) ); span *= radius_; - bounds().min() = origin_ - span; - bounds().max() = origin_ + span; + bounds().min() = origin() - span; + bounds().max() = origin() + span; } @@ -169,28 +158,14 @@ Foam::searchableDisk::searchableDisk const dictionary& dict ) : - searchableSurface(io), - origin_(dict.get("origin")), - normal_(dict.get("normal")), - radius_(dict.get("radius")) -{ - normal_ /= mag(normal_); - - // Rough approximation of bounding box - //vector span(radius_, radius_, radius_); - - // See searchableCylinder - vector span + searchableDisk ( - sqrt(sqr(normal_.y()) + sqr(normal_.z())), - sqrt(sqr(normal_.x()) + sqr(normal_.z())), - sqrt(sqr(normal_.x()) + sqr(normal_.y())) - ); - span *= radius_; - - bounds().min() = origin_ - span; - bounds().max() = origin_ + span; -} + io, + dict.get("origin"), + dict.get("normal"), + dict.get("radius") + ) +{} // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // @@ -213,7 +188,7 @@ void Foam::searchableDisk::boundingSpheres ) const { centres.setSize(1); - centres[0] = origin_; + centres[0] = origin(); radiusSqr.setSize(1); radiusSqr[0] = sqr(radius_); @@ -299,7 +274,7 @@ void Foam::searchableDisk::getRegion labelList& region ) const { - region.setSize(info.size()); + region.resize(info.size()); region = 0; } @@ -307,11 +282,11 @@ void Foam::searchableDisk::getRegion void Foam::searchableDisk::getNormal ( const List& info, - vectorField& normal + vectorField& normals ) const { - normal.setSize(info.size()); - normal = normal_; + normals.resize(info.size()); + normals = normal(); } diff --git a/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H b/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H index dbe3f9d247..cc589aeb4e 100644 --- a/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H +++ b/src/meshTools/searchableSurfaces/searchableDisk/searchableDisk.H @@ -45,7 +45,7 @@ SourceFiles #ifndef searchableDisk_H #define searchableDisk_H -#include "treeBoundBox.H" +#include "plane.H" #include "searchableSurface.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -59,18 +59,13 @@ namespace Foam class searchableDisk : - public searchableSurface + public searchableSurface, + public plane { private: // Private Member Data - //- Origin - const point origin_; - - //- Normal - vector normal_; - //- Radius const scalar radius_; @@ -116,8 +111,8 @@ public: searchableDisk ( const IOobject& io, - const point& origin, - const point& normal, + const point& originPoint, + const vector& normalVector, const scalar radius ); @@ -128,6 +123,7 @@ public: const dictionary& dict ); + //- Destructor virtual ~searchableDisk() = default; @@ -153,7 +149,7 @@ public: // Usually the element centres (should be of length size()). virtual tmp coordinates() const { - return tmp::New(1, origin_); + return tmp::New(one(), origin()); } //- Get bounding spheres (centre and radius squared), one per element. @@ -220,7 +216,7 @@ public: virtual void getNormal ( const List&, - vectorField& normal + vectorField& normals ) const; //- Determine type (inside/outside/mixed) for point. diff --git a/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.C b/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.C index 9a9577c7f9..c420c0d04c 100644 --- a/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.C +++ b/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.C @@ -145,7 +145,7 @@ void Foam::searchablePlane::boundingSpheres ) const { centres.setSize(1); - centres[0] = refPoint(); + centres[0] = origin(); radiusSqr.setSize(1); radiusSqr[0] = Foam::sqr(GREAT); diff --git a/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H b/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H index 33518a4668..4eae9e5d01 100644 --- a/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H +++ b/src/meshTools/searchableSurfaces/searchablePlane/searchablePlane.H @@ -134,7 +134,7 @@ public: // Usually the element centres (should be of length size()). virtual tmp coordinates() const { - return tmp::New(1, refPoint()); + return tmp::New(1, origin()); } //- Get bounding spheres (centre and radius squared), one per element. diff --git a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C index 5b65633892..75808f63ea 100644 --- a/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C +++ b/src/meshTools/triSurface/triSurfaceTools/triSurfaceTools.C @@ -855,7 +855,7 @@ Foam::surfaceLocation Foam::triSurfaceTools::cutEdge scalar norm = 0; forAll(d, fp) { - d[fp] = (points[f[fp]]-cutPlane.refPoint()) & cutPlane.normal(); + d[fp] = cutPlane.signedDistance(points[f[fp]]); norm += mag(d[fp]); } diff --git a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C index 2c000f9b1d..6dc7894ac6 100644 --- a/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C +++ b/src/sampling/sampledSurface/sampledCuttingPlane/sampledCuttingPlane.C @@ -126,8 +126,7 @@ void Foam::sampledCuttingPlane::createGeometry() forAll(cc, i) { - // Signed distance - fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal(); + fld[i] = plane_.signedDistance(cc[i]); } } @@ -163,7 +162,7 @@ void Foam::sampledCuttingPlane::createGeometry() fld.setSize(pp.size()); forAll(fld, i) { - fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal(); + fld[i] = plane_.signedDistance(cc[i]); } } else @@ -173,7 +172,7 @@ void Foam::sampledCuttingPlane::createGeometry() forAll(fld, i) { - fld[i] = (cc[i] - plane_.refPoint()) & plane_.normal(); + fld[i] = plane_.signedDistance(cc[i]); } } } @@ -191,7 +190,7 @@ void Foam::sampledCuttingPlane::createGeometry() forAll(pointDistance_, i) { - pointDistance_[i] = (pts[i] - plane_.refPoint()) & plane_.normal(); + pointDistance_[i] = plane_.signedDistance(pts[i]); } } diff --git a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C index d89c3e50ce..8d2b895ed9 100644 --- a/src/sampling/sampledSurface/sampledPlane/sampledPlane.C +++ b/src/sampling/sampledSurface/sampledPlane/sampledPlane.C @@ -90,7 +90,7 @@ Foam::sampledPlane::sampledPlane { coordinateSystem cs(mesh, dict.subDict("coordinateSystem")); - const point base = cs.globalPosition(planeDesc().refPoint()); + const point base = cs.globalPosition(planeDesc().origin()); const vector norm = cs.globalVector(planeDesc().normal()); // Assign the plane description @@ -334,8 +334,8 @@ Foam::tmp Foam::sampledPlane::interpolate void Foam::sampledPlane::print(Ostream& os) const { os << "sampledPlane: " << name() << " :" - << " base:" << refPoint() - << " normal:" << normal() + << " base:" << plane::origin() + << " normal:" << plane::normal() << " triangulate:" << triangulate_ << " faces:" << faces().size() << " points:" << points().size(); diff --git a/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C b/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C index e1b8369dff..9b05c160c7 100644 --- a/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C +++ b/src/sampling/surfMeshSample/plane/surfMeshSamplePlane.C @@ -91,7 +91,7 @@ Foam::surfMeshSamplePlane::surfMeshSamplePlane { coordinateSystem cs(mesh, dict.subDict("coordinateSystem")); - const point base = cs.globalPosition(planeDesc().refPoint()); + const point base = cs.globalPosition(planeDesc().origin()); const vector norm = cs.globalVector(planeDesc().normal()); // Assign the plane description @@ -263,8 +263,8 @@ bool Foam::surfMeshSamplePlane::sample void Foam::surfMeshSamplePlane::print(Ostream& os) const { os << "surfMeshSamplePlane: " << name() << " :" - << " base:" << cuttingPlane::refPoint() - << " normal:" << cuttingPlane::normal() + << " base:" << plane::origin() + << " normal:" << plane::normal() << " triangulate:" << triangulate_ << " faces:" << SurfaceSource::surfFaces().size() << " points:" << SurfaceSource::points().size(); diff --git a/src/sampling/surface/cuttingPlane/cuttingPlane.C b/src/sampling/surface/cuttingPlane/cuttingPlane.C index d257bf097f..580f817aaa 100644 --- a/src/sampling/surface/cuttingPlane/cuttingPlane.C +++ b/src/sampling/surface/cuttingPlane/cuttingPlane.C @@ -86,7 +86,7 @@ void Foam::cuttingPlane::calcCutCells } // Set correct list size - meshCells_.setSize(cutcelli); + meshCells_.resize(cutcelli); } @@ -102,7 +102,7 @@ void Foam::cuttingPlane::intersectEdges const pointField& points = mesh.points(); // Per edge -1 or the label of the intersection point - edgePoint.setSize(edges.size()); + edgePoint.resize(edges.size()); DynamicList dynCuttingPoints(4*meshCells_.size()); diff --git a/src/sampling/surface/isoSurface/isoSurface.C b/src/sampling/surface/isoSurface/isoSurface.C index c5ae60b669..90d185090b 100644 --- a/src/sampling/surface/isoSurface/isoSurface.C +++ b/src/sampling/surface/isoSurface/isoSurface.C @@ -1080,7 +1080,7 @@ void Foam::isoSurface::trimToPlanes forAll(planes, faceI) { - const plane& pl = planes[faceI]; + const plane& pln = planes[faceI]; if (useA) { @@ -1088,7 +1088,7 @@ void Foam::isoSurface::trimToPlanes forAll(insideOpA.tris_, i) { const triPoints& tri = insideOpA.tris_[i]; - triPointRef(tri).sliceWithPlane(pl, insideOpB, dop); + triPointRef(tri).sliceWithPlane(pln, insideOpB, dop); } } else @@ -1097,7 +1097,7 @@ void Foam::isoSurface::trimToPlanes forAll(insideOpB.tris_, i) { const triPoints& tri = insideOpB.tris_[i]; - triPointRef(tri).sliceWithPlane(pl, insideOpA, dop); + triPointRef(tri).sliceWithPlane(pln, insideOpA, dop); } } useA = !useA; @@ -1141,13 +1141,11 @@ void Foam::isoSurface::trimToBox } // Generate inwards pointing planes - PtrList planes(6); - const pointField pts(bb.treeBoundBox::points()); - forAll(treeBoundBox::faces, faceI) + PtrList planes(treeBoundBox::faceNormals.size()); + forAll(treeBoundBox::faceNormals, faceI) { - const face& f = treeBoundBox::faces[faceI]; const vector& n = treeBoundBox::faceNormals[faceI]; - planes.set(faceI, new plane(pts[f[0]], -n)); + planes.set(faceI, new plane(bb.faceCentre(faceI), -n)); } label nTris = triPoints.size()/3; diff --git a/tutorials/incompressible/pimpleFoam/LES/vortexShed/system/controlDict b/tutorials/incompressible/pimpleFoam/LES/vortexShed/system/controlDict index b47233563a..564afa824e 100644 --- a/tutorials/incompressible/pimpleFoam/LES/vortexShed/system/controlDict +++ b/tutorials/incompressible/pimpleFoam/LES/vortexShed/system/controlDict @@ -94,8 +94,8 @@ functions planeType pointAndNormal; pointAndNormalDict { - basePoint (0 0 -0.01); - normalVector (0 0 1); + point (0 0 -0.01); + normal (0 0 1); } interpolate false; }