diff --git a/src/mesh/blockMesh/Make/files b/src/mesh/blockMesh/Make/files index d17f08e04..4d1da7206 100644 --- a/src/mesh/blockMesh/Make/files +++ b/src/mesh/blockMesh/Make/files @@ -1,5 +1,6 @@ blockVertices/blockVertex/blockVertex.C blockVertices/pointVertex/pointVertex.C +blockVertices/projectVertex/projectVertex.C blockEdges/blockEdge/blockEdge.C blockEdges/lineDivide/lineDivide.C @@ -11,6 +12,7 @@ blockEdges/BSplineEdge/BSpline.C blockEdges/BSplineEdge/BSplineEdge.C blockEdges/splineEdge/CatmullRomSpline.C blockEdges/splineEdge/splineEdge.C +blockEdges/projectEdge/projectEdge.C blockFaces/blockFace/blockFace.C blockFaces/projectFace/projectFace.C diff --git a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C index 7ccf24d42..fb9a41301 100644 --- a/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C +++ b/src/mesh/blockMesh/blockDescriptor/blockDescriptor.C @@ -336,7 +336,12 @@ void Foam::blockDescriptor::correctFacePoints { if (curvedFaces_[blockFacei] != -1) { - faces_[curvedFaces_[blockFacei]].project(facePoints[blockFacei]); + faces_[curvedFaces_[blockFacei]].project + ( + *this, + blockFacei, + facePoints[blockFacei] + ); } } } diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C index a2cc74182..4ed2647a7 100644 --- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C +++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.C @@ -125,6 +125,20 @@ Foam::pointField Foam::blockEdge::appendEndPoints } +Foam::tmp +Foam::blockEdge::position(const scalarList& lambdas) const +{ + tmp tpoints(new pointField(lambdas.size())); + pointField& points = tpoints.ref(); + + forAll(lambdas, i) + { + points[i] = position(lambdas[i]); + } + return tpoints; +} + + // * * * * * * * * * * * * * * * Member Operators * * * * * * * * * * * * * // Foam::Ostream& Foam::operator<<(Ostream& os, const blockEdge& p) diff --git a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H index 86561f6ab..e6e3196ea 100644 --- a/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H +++ b/src/mesh/blockMesh/blockEdges/blockEdge/blockEdge.H @@ -188,6 +188,10 @@ public: // 0 <= lambda <= 1 virtual point position(const scalar) const = 0; + //- Return the point positions corresponding to the curve parameters + // 0 <= lambda <= 1 + virtual tmp position(const scalarList&) const; + //- Return the length of the curve virtual scalar length() const = 0; diff --git a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C index f0c6ccee5..5c09c0834 100644 --- a/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C +++ b/src/mesh/blockMesh/blockEdges/lineDivide/lineDivide.C @@ -131,10 +131,7 @@ Foam::lineDivide::lineDivide } // Calculate the points - for (label i = 0; i <= nDiv; i++) - { - points_[i] = cedge.position(divisions_[i]); - } + points_ = cedge.position(divisions_); } diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C new file mode 100644 index 000000000..b604e6ffb --- /dev/null +++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.C @@ -0,0 +1,153 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +\*---------------------------------------------------------------------------*/ + +#include "searchableSurfacesQueries.H" +#include "projectEdge.H" +#include "unitConversion.H" +#include "addToRunTimeSelectionTable.H" +#include "pointConstraint.H" +#include "plane.H" + +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // + +namespace Foam +{ + defineTypeNameAndDebug(projectEdge, 0); + addToRunTimeSelectionTable(blockEdge, projectEdge, Istream); +} + + +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // + +void Foam::projectEdge::findNearest +( + const point& pt, + point& near, + pointConstraint& constraint +) const +{ + if (surfaces_.size()) + { + const scalar distSqr = magSqr(points_[end_]-points_[start_]); + + pointField boundaryNear(1); + List boundaryConstraint(1); + + searchableSurfacesQueries::findNearest + ( + geometry_, + surfaces_, + pointField(1, pt), + scalarField(1, distSqr), + boundaryNear, + boundaryConstraint + ); + near = boundaryNear[0]; + constraint = boundaryConstraint[0]; + } + else + { + near = pt; + constraint = pointConstraint(); + } +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +Foam::projectEdge::projectEdge +( + const searchableSurfaces& geometry, + const pointField& points, + Istream& is +) +: + blockEdge(points, is), + geometry_(geometry) +{ + wordList names(is); + surfaces_.setSize(names.size()); + forAll(names, i) + { + surfaces_[i] = geometry_.findSurfaceID(names[i]); + + if (surfaces_[i] == -1) + { + FatalIOErrorInFunction(is) + << "Cannot find surface " << names[i] << " in geometry" + << exit(FatalIOError); + } + } +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::point Foam::projectEdge::position(const scalar lambda) const +{ + // Initial guess + const point start(points_[start_]+lambda*(points_[end_]-points_[start_])); + + point near(start); + + if (lambda >= SMALL && lambda < 1.0-SMALL) + { + pointConstraint constraint; + findNearest(start, near, constraint); + } + + return near; +} + + +Foam::tmp +Foam::projectEdge::position(const scalarList& lambdas) const +{ + tmp tpoints(new pointField(lambdas.size())); + pointField& points = tpoints.ref(); + + const point& startPt = points_[start_]; + const point& endPt = points_[end_]; + const vector d = endPt-startPt; + + forAll(lambdas, i) + { + points[i] = startPt+lambdas[i]*d; + } + + forAll(lambdas, i) + { + if (lambdas[i] >= SMALL && lambdas[i] < 1.0-SMALL) + { + pointConstraint constraint; + findNearest(points[i], points[i], constraint); + } + } + + return tpoints; +} + + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H new file mode 100644 index 000000000..cf8baebe7 --- /dev/null +++ b/src/mesh/blockMesh/blockEdges/projectEdge/projectEdge.H @@ -0,0 +1,125 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2016 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 . + +Class + Foam::projectEdge + +Description + Defines the edge from the projection onto a surface (single surface) + or intersection of two surfaces. + +SourceFiles + projectEdge.C + +\*---------------------------------------------------------------------------*/ + +#ifndef projectEdge_H +#define projectEdge_H + +#include "blockEdge.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +class pointConstraint; + +/*---------------------------------------------------------------------------*\ + Class projectEdge Declaration +\*---------------------------------------------------------------------------*/ + +class projectEdge +: + public blockEdge +{ + // Private data + + const searchableSurfaces& geometry_; + + //- The indices of surfaces onto which the points are projected + labelList surfaces_; + + + // Private Member Functions + + //- Single point find nearest + void findNearest(const point&, point& near, pointConstraint&) const; + + //- Disallow default bitwise copy construct + projectEdge(const projectEdge&); + + //- Disallow default bitwise assignment + void operator=(const projectEdge&); + + +public: + + //- Runtime type information + TypeName("project"); + + + // Constructors + + //- Construct from Istream setting pointsList + projectEdge + ( + const searchableSurfaces& geometry, + const pointField& points, + Istream& + ); + + + //- Destructor + virtual ~projectEdge() + {} + + + // Member Functions + + //- Return the point positions corresponding to the curve parameters + // 0 <= lambda <= 1 + virtual point position(const scalar) const; + + //- Return the point positions corresponding to the curve parameters + // 0 <= lambda <= 1 + virtual tmp position(const scalarList&) const; + + //- Return the length of the curve + virtual scalar length() const + { + NotImplemented; + return 1; + } +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H index b060031d5..fee6771d7 100644 --- a/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H +++ b/src/mesh/blockMesh/blockFaces/blockFace/blockFace.H @@ -44,11 +44,11 @@ namespace Foam // Forward declaration of friend functions and operators +class blockDescriptor; class blockFace; Ostream& operator<<(Ostream&, const blockFace&); - /*---------------------------------------------------------------------------*\ Class blockFace Declaration \*---------------------------------------------------------------------------*/ @@ -138,7 +138,12 @@ public: //- Compare with the given block and block face inline bool compare(const face& vertices) const; - virtual void project(pointField& points) const = 0; + virtual void project + ( + const blockDescriptor&, + const label blockFacei, + pointField& points + ) const = 0; // Ostream operator diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C index fbc40094d..2cf9c9742 100644 --- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C +++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.C @@ -26,6 +26,7 @@ License #include "projectFace.H" #include "unitConversion.H" #include "addToRunTimeSelectionTable.H" +#include "blockDescriptor.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // @@ -80,7 +81,12 @@ Foam::blockFaces::projectFace::projectFace // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // -void Foam::blockFaces::projectFace::project(pointField& points) const +void Foam::blockFaces::projectFace::project +( + const blockDescriptor& desc, + const label blockFacei, + pointField& points +) const { List hits; scalarField nearestDistSqr diff --git a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H index 44a2a25c5..c1b56bbcd 100644 --- a/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H +++ b/src/mesh/blockMesh/blockFaces/projectFace/projectFace.H @@ -98,7 +98,12 @@ public: // Member Functions //- Project the given points onto the surface - virtual void project(pointField& points) const; + virtual void project + ( + const blockDescriptor&, + const label blockFacei, + pointField& points + ) const; }; diff --git a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H index 9df74a713..1179c66bd 100644 --- a/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H +++ b/src/mesh/blockMesh/blockVertices/pointVertex/pointVertex.H @@ -51,7 +51,9 @@ class pointVertex : public blockVertex { - // Private member data +protected: + + // Protected member data //- The vertex location point vertex_; diff --git a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C index a3d6c28d7..04cff7d25 100644 --- a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C +++ b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.C @@ -23,78 +23,77 @@ License \*---------------------------------------------------------------------------*/ -#include "projectFace.H" +#include "projectVertex.H" #include "unitConversion.H" #include "addToRunTimeSelectionTable.H" +#include "searchableSurfacesQueries.H" +#include "pointConstraint.H" // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // namespace Foam { - defineTypeNameAndDebug(projectFace, 0); - addToRunTimeSelectionTable(blockVertex, projectFace, Istream); -} - - -// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // - -const Foam::searchableSurface& Foam::projectFace::lookupSurface -( - const searchableSurfaces& geometry, - Istream& is -) const +namespace blockVertices { - word name(is); - - forAll(geometry, i) - { - if (geometry[i].name() == name) - { - return geometry[i]; - } - } - - FatalIOErrorInFunction(is) - << "Cannot find surface " << name << " in geometry" - << exit(FatalIOError); - - return geometry[0]; + defineTypeNameAndDebug(projectVertex, 0); + addToRunTimeSelectionTable(blockVertex, projectVertex, Istream); +} } // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // -Foam::projectFace::projectFace +Foam::blockVertices::projectVertex::projectVertex ( const searchableSurfaces& geometry, Istream& is ) : - blockVertex(is), - surface_(lookupSurface(geometry, is)) -{} - - -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -void Foam::projectFace::project(pointField& points) const + pointVertex(geometry, is), + geometry_(geometry) { - List hits; - scalarField nearestDistSqr - ( - points.size(), - magSqr(points[0] - points[points.size()-1]) - ); - surface_.findNearest(points, nearestDistSqr, hits); - - forAll(hits, i) + wordList names(is); + surfaces_.setSize(names.size()); + forAll(names, i) { - if (hits[i].hit()) + surfaces_[i] = geometry_.findSurfaceID(names[i]); + + if (surfaces_[i] == -1) { - points[i] = hits[i].hitPoint(); + FatalIOErrorInFunction(is) + << "Cannot find surface " << names[i] << " in geometry" + << exit(FatalIOError); } } } +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +Foam::blockVertices::projectVertex::operator point() const +{ + pointField start(1, pointVertex::operator point()); + + pointField boundaryNear(start); + List boundaryConstraint; + + + // Note: how far do we need to search? Probably not further than + // span of surfaces themselves. + boundBox bb(searchableSurfacesQueries::bounds(geometry_, surfaces_)); + + searchableSurfacesQueries::findNearest + ( + geometry_, + surfaces_, + start, + scalarField(start.size(), magSqr(bb.span())), + boundaryNear, + boundaryConstraint + ); + + return boundaryNear[0]; +} + + // ************************************************************************* // diff --git a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H index 23e4fd10b..929b6a68e 100644 --- a/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H +++ b/src/mesh/blockMesh/blockVertices/projectVertex/projectVertex.H @@ -22,54 +22,52 @@ License along with OpenFOAM. If not, see . Class - Foam::projectFace + Foam::blockVertices::projectVertex Description - Projects the given set of face points onto the selected surface of the + Projects the vertex onto the selected surfacees of the geometry provided as a searchableSurfaces object. SourceFiles - projectFace.C + projectVertex.C \*---------------------------------------------------------------------------*/ -#ifndef projectFace_H -#define projectFace_H +#ifndef blockVertices_projectVertex_H +#define blockVertices_projectVertex_H -#include "blockVertex.H" +#include "pointVertex.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // namespace Foam { +namespace blockVertices +{ /*---------------------------------------------------------------------------*\ - Class projectFace Declaration + Class projectVertex Declaration \*---------------------------------------------------------------------------*/ -class projectFace +class projectVertex : - public blockVertex + public pointVertex { // Private data - //- The surface onto which the points are projected - const searchableSurface& surface_; + const searchableSurfaces& geometry_; + + //- The indices of surfaces onto which the points are projected + labelList surfaces_; // Private Member Functions - const searchableSurface& lookupSurface - ( - const searchableSurfaces& geometry, - Istream& is - ) const; - //- Disallow default bitwise copy construct - projectFace(const projectFace&); + projectVertex(const projectVertex&); //- Disallow default bitwise assignment - void operator=(const projectFace&); + void operator=(const projectVertex&); public: @@ -81,7 +79,7 @@ public: // Constructors //- Construct from Istream setting pointsList - projectFace + projectVertex ( const searchableSurfaces& geometry, Istream& @@ -89,19 +87,20 @@ public: //- Destructor - virtual ~projectFace() + virtual ~projectVertex() {} // Member Functions //- Project the given points onto the surface - virtual void project(pointField& points) const; + virtual operator point() const; }; // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // +} // End namespace blockVertices } // End namespace Foam // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // diff --git a/src/meshTools/searchableSurface/searchableSurfaces.C b/src/meshTools/searchableSurface/searchableSurfaces.C index e25c1c009..6800034d0 100644 --- a/src/meshTools/searchableSurface/searchableSurfaces.C +++ b/src/meshTools/searchableSurface/searchableSurfaces.C @@ -419,24 +419,6 @@ Foam::boundBox Foam::searchableSurfaces::bounds() const } -Foam::pointIndexHit Foam::searchableSurfaces::facesIntersection -( - const scalar initDistSqr, - const scalar convergenceDistSqr, - const point& start -) const -{ - return searchableSurfacesQueries::facesIntersection - ( - *this, - allSurfaces_, - initDistSqr, - convergenceDistSqr, - start - ); -} - - bool Foam::searchableSurfaces::checkClosed(const bool report) const { if (report) diff --git a/src/meshTools/searchableSurface/searchableSurfaces.H b/src/meshTools/searchableSurface/searchableSurfaces.H index b367a2a0d..4efd019eb 100644 --- a/src/meshTools/searchableSurface/searchableSurfaces.H +++ b/src/meshTools/searchableSurface/searchableSurfaces.H @@ -213,18 +213,6 @@ public: //- Calculate bounding box boundBox bounds() const; - - // Single point queries - - //- Calculate point which is on a set of surfaces. - pointIndexHit facesIntersection - ( - const scalar initialDistSqr, - const scalar convergenceDistSqr, - const point& start - ) const; - - // Checking //- Are all surfaces closed and manifold diff --git a/src/meshTools/searchableSurface/searchableSurfacesQueries.C b/src/meshTools/searchableSurface/searchableSurfacesQueries.C index 26abe9bde..818ff6648 100644 --- a/src/meshTools/searchableSurface/searchableSurfacesQueries.C +++ b/src/meshTools/searchableSurface/searchableSurfacesQueries.C @@ -28,6 +28,8 @@ License #include "OFstream.H" #include "meshTools.H" #include "DynamicField.H" +#include "pointConstraint.H" +#include "plane.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -39,288 +41,6 @@ namespace Foam // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // -Foam::pointIndexHit Foam::searchableSurfacesQueries::tempFindNearest -( - const searchableSurface& surf, - const point& pt, - const scalar initDistSqr -) -{ - pointField onePoint(1, pt); - scalarField oneDist(1, initDistSqr); - List oneHit(1); - surf.findNearest(onePoint, oneDist, oneHit); - return oneHit[0]; -} - - -// Calculate sum of distance to surfaces. -Foam::scalar Foam::searchableSurfacesQueries::sumDistSqr -( - const PtrList& allSurfaces, - const labelList& surfacesToTest, - const scalar initDistSqr, - const point& pt -) -{ - scalar sum = 0; - - forAll(surfacesToTest, testI) - { - label surfI = surfacesToTest[testI]; - - pointIndexHit hit - ( - tempFindNearest(allSurfaces[surfI], pt, initDistSqr) - ); - - // Note: make it fall over if not hit. - sum += magSqr(hit.hitPoint()-pt); - } - return sum; -} - - -// Reflects the point furthest away around the triangle centre by a factor fac. -// (triangle centre is the average of all points but the ihi. pSum is running -// sum of all points) -Foam::scalar Foam::searchableSurfacesQueries::tryMorphTet -( - const PtrList& allSurfaces, - const labelList& surfacesToTest, - const scalar initDistSqr, - List& p, - List& y, - vector& pSum, - const label ihi, - const scalar fac -) -{ - scalar fac1 = (1.0-fac)/vector::nComponents; - scalar fac2 = fac1-fac; - - vector ptry = pSum*fac1-p[ihi]*fac2; - - scalar ytry = sumDistSqr(allSurfaces, surfacesToTest, initDistSqr, ptry); - - if (ytry < y[ihi]) - { - y[ihi] = ytry; - pSum += ptry - p[ihi]; - p[ihi] = ptry; - } - return ytry; -} - - -bool Foam::searchableSurfacesQueries::morphTet -( - const PtrList& allSurfaces, - const labelList& surfacesToTest, - const scalar initDistSqr, - const scalar convergenceDistSqr, - const label maxIter, - List& p, - List& y -) -{ - vector pSum = sum(p); - - autoPtr str; - label vertI = 0; - if (debug) - { - wordList names(surfacesToTest.size()); - forAll(surfacesToTest, i) - { - names[i] = allSurfaces[surfacesToTest[i]].name(); - } - Pout<< "searchableSurfacesQueries::morphTet : intersection of " - << names << " starting from points:" << p << endl; - str.reset(new OFstream("track.obj")); - meshTools::writeOBJ(str(), p[0]); - vertI++; - } - - for (label iter = 0; iter < maxIter; iter++) - { - // Get the indices of lowest, highest and second-highest values. - label ilo, ihi, inhi; - { - labelList sortedIndices; - sortedOrder(y, sortedIndices); - ilo = sortedIndices[0]; - ihi = sortedIndices[sortedIndices.size()-1]; - inhi = sortedIndices[sortedIndices.size()-2]; - } - - if (debug) - { - Pout<< "Iteration:" << iter - << " lowest:" << y[ilo] << " highest:" << y[ihi] - << " points:" << p << endl; - - meshTools::writeOBJ(str(), p[ilo]); - vertI++; - str()<< "l " << vertI-1 << ' ' << vertI << nl; - } - - if (y[ihi] < convergenceDistSqr) - { - // Get point on 0th surface. - Swap(p[0], p[ilo]); - Swap(y[0], y[ilo]); - return true; - } - - // Reflection: point furthest away gets reflected. - scalar ytry = tryMorphTet - ( - allSurfaces, - surfacesToTest, - 10*y[ihi], // search box. - p, - y, - pSum, - ihi, - -1.0 - ); - - if (ytry <= y[ilo]) - { - // If in right direction (y lower) expand by two. - ytry = tryMorphTet - ( - allSurfaces, - surfacesToTest, - 10*y[ihi], - p, - y, - pSum, - ihi, - 2.0 - ); - } - else if (ytry >= y[inhi]) - { - // If inside tet try contraction. - - scalar ysave = y[ihi]; - - ytry = tryMorphTet - ( - allSurfaces, - surfacesToTest, - 10*y[ihi], - p, - y, - pSum, - ihi, - 0.5 - ); - - if (ytry >= ysave) - { - // Contract around lowest point. - forAll(p, i) - { - if (i != ilo) - { - p[i] = 0.5*(p[i] + p[ilo]); - y[i] = sumDistSqr - ( - allSurfaces, - surfacesToTest, - y[ihi], - p[i] - ); - } - } - pSum = sum(p); - } - } - } - - if (debug) - { - meshTools::writeOBJ(str(), p[0]); - vertI++; - str()<< "l " << vertI-1 << ' ' << vertI << nl; - } - - // Failure to converge. Return best guess so far. - label ilo = findMin(y); - Swap(p[0], p[ilo]); - Swap(y[0], y[ilo]); - return false; -} - - -//// Get all intersections (in order) for single surface. -//void Foam::searchableSurfacesQueries::findAllIntersections -//( -// const searchableSurface& s, -// const pointField& start, -// const pointField& end, -// const vectorField& smallVec, -// List>& surfaceHitInfo -//) -//{ -// surfaceHitInfo.setSize(start.size()); -// -// // Current start point of vector -// pointField p0(start); -// -// List intersectInfo(start.size()); -// -// // For test whether finished doing vector. -// const vectorField dirVec(end-start); -// const scalarField magSqrDirVec(magSqr(dirVec)); -// -// while (true) -// { -// // Find first intersection. Synced. -// s.findLine(p0, end, intersectInfo); -// -// label nHits = 0; -// -// forAll(intersectInfo, i) -// { -// if (intersectInfo[i].hit()) -// { -// nHits++; -// -// label sz = surfaceHitInfo[i].size(); -// surfaceHitInfo[i].setSize(sz+1); -// surfaceHitInfo[i][sz] = intersectInfo[i]; -// -// p0[i] = intersectInfo[i].hitPoint() + smallVec[i]; -// -// // If beyond endpoint set to endpoint so as not to pick up -// // any intersections. Could instead just filter out hits. -// if (((p0[i]-start[i])&dirVec[i]) > magSqrDirVec[i]) -// { -// p0[i] = end[i]; -// } -// } -// else -// { -// // Set to endpoint to stop intersection test. See above. -// p0[i] = end[i]; -// } -// } -// -// // returnReduce(nHits) ? -// if (nHits == 0) -// { -// break; -// } -// } -//} - - -// Given current set of hits (allSurfaces, allInfo) merge in those coming from -// surface surfI. void Foam::searchableSurfacesQueries::mergeHits ( const point& start, @@ -333,6 +53,9 @@ void Foam::searchableSurfacesQueries::mergeHits scalarList& allDistSqr ) { + // Given current set of hits (allSurfaces, allInfo) merge in those coming + // from surface surfI. + // Precalculate distances scalarList surfDistSqr(surfHits.size()); forAll(surfHits, i) @@ -532,7 +255,6 @@ void Foam::searchableSurfacesQueries::findAllIntersections } -//Find intersections of edge nearest to both endpoints. void Foam::searchableSurfacesQueries::findNearestIntersection ( const PtrList& allSurfaces, @@ -620,7 +342,6 @@ void Foam::searchableSurfacesQueries::findNearestIntersection } -// Find nearest. Return -1 or nearest point void Foam::searchableSurfacesQueries::findNearest ( const PtrList& allSurfaces, @@ -631,6 +352,8 @@ void Foam::searchableSurfacesQueries::findNearest List& nearestInfo ) { + // Find nearest. Return -1 or nearest point + // Initialise nearestSurfaces.setSize(samples.size()); nearestSurfaces = -1; @@ -667,7 +390,6 @@ void Foam::searchableSurfacesQueries::findNearest } -// Find nearest. Return -1 or nearest point void Foam::searchableSurfacesQueries::findNearest ( const PtrList& allSurfaces, @@ -679,6 +401,8 @@ void Foam::searchableSurfacesQueries::findNearest List& nearestInfo ) { + // Find nearest. Return -1 or nearest point + if (regionIndices.empty()) { findNearest @@ -729,6 +453,105 @@ void Foam::searchableSurfacesQueries::findNearest } +void Foam::searchableSurfacesQueries::findNearest +( + const PtrList& allSurfaces, + const labelList& surfacesToTest, + const pointField& start, + const scalarField& distSqr, + pointField& near, + List& constraint, + const label nIter +) +{ + // Multi-surface findNearest + + vectorField normal; + List info; + + allSurfaces[surfacesToTest[0]].findNearest(start, distSqr, info); + allSurfaces[surfacesToTest[0]].getNormal(info, normal); + + // Extract useful info from initial start point + near = start; + forAll(info, i) + { + if (info[i].hit()) + { + near[i] = info[i].hitPoint(); + } + } + constraint.setSize(near.size()); + + + if (surfacesToTest.size() == 1) + { + constraint = pointConstraint(); + forAll(info, i) + { + if (info[i].hit()) + { + constraint[i].applyConstraint(normal[i]); + } + } + } + else if (surfacesToTest.size() >= 2) + { + // Work space + pointField near1; + vectorField normal1; + + label surfi = 1; + for (label iter = 0; iter < nIter; iter++) + { + constraint = pointConstraint(); + forAll(constraint, i) + { + if (info[i].hit()) + { + constraint[i].applyConstraint(normal[i]); + } + } + + // Find intersection with next surface + const searchableSurface& s = allSurfaces[surfacesToTest[surfi]]; + s.findNearest(near, distSqr, info); + s.getNormal(info, normal1); + near1.setSize(info.size()); + forAll(info, i) + { + if (info[i].hit()) + { + near1[i] = info[i].hitPoint(); + } + } + + // Move to intersection + forAll(near, pointi) + { + if (info[pointi].hit()) + { + plane pl0(near[pointi], normal[pointi]); + plane pl1(near1[pointi], normal1[pointi]); + plane::ray r(pl0.planeIntersect(pl1)); + vector n = r.dir() / mag(r.dir()); + + vector d(r.refPoint()-near[pointi]); + d -= (d&n)*n; + + near[pointi] += d; + normal[pointi] = normal1[pointi]; + constraint[pointi].applyConstraint(normal1[pointi]); + } + } + + // Step to next surface + surfi = surfacesToTest.fcIndex(surfi); + } + } +} + + void Foam::searchableSurfacesQueries::signedDistance ( const PtrList& allSurfaces, @@ -849,126 +672,4 @@ Foam::boundBox Foam::searchableSurfacesQueries::bounds } -Foam::pointIndexHit Foam::searchableSurfacesQueries::facesIntersection -( - const PtrList& allSurfaces, - const labelList& surfacesToTest, - const scalar initDistSqr, - const scalar convergenceDistSqr, - const point& start -) -{ - // Get four starting points. Take these as the projection of the - // starting point onto the surfaces and the mid point - List nearest(surfacesToTest.size()+1); - - point sumNearest = Zero; - - forAll(surfacesToTest, i) - { - pointIndexHit hit - ( - tempFindNearest(allSurfaces[surfacesToTest[i]], start, initDistSqr) - ); - - if (hit.hit()) - { - nearest[i] = hit.hitPoint(); - sumNearest += nearest[i]; - } - else - { - FatalErrorInFunction - << "Did not find point within distance " - << initDistSqr << " of starting point " << start - << " on surface " - << allSurfaces[surfacesToTest[i]].IOobject::name() - << abort(FatalError); - } - } - - nearest.last() = sumNearest / surfacesToTest.size(); - - - // Get the sum of distances (initial evaluation) - List nearestDist(nearest.size()); - - forAll(nearestDist, i) - { - nearestDist[i] = sumDistSqr - ( - allSurfaces, - surfacesToTest, - initDistSqr, - nearest[i] - ); - } - - - //- Downhill Simplex method - - bool converged = morphTet - ( - allSurfaces, - surfacesToTest, - initDistSqr, - convergenceDistSqr, - 2000, - nearest, - nearestDist - ); - - - pointIndexHit intersection; - - if (converged) - { - // Project nearest onto 0th surface. - intersection = tempFindNearest - ( - allSurfaces[surfacesToTest[0]], - nearest[0], - nearestDist[0] - ); - } - - //if (!intersection.hit()) - //{ - // // Restart - // scalar smallDist = Foam::sqr(convergenceDistSqr); - // nearest[0] = intersection.hitPoint(); - // nearest[1] = nearest[0]; - // nearest[1].x() += smallDist; - // nearest[2] = nearest[0]; - // nearest[2].y() += smallDist; - // nearest[3] = nearest[0]; - // nearest[3].z() += smallDist; - // - // forAll(nearestDist, i) - // { - // nearestDist[i] = sumDistSqr - // ( - // surfacesToTest, - // initDistSqr, - // nearest[i] - // ); - // } - // - // intersection = morphTet - // ( - // allSurfaces, - // surfacesToTest, - // initDistSqr, - // convergenceDistSqr, - // 1000, - // nearest, - // nearestDist - // ); - //} - - return intersection; -} - - - // ************************************************************************* // diff --git a/src/meshTools/searchableSurface/searchableSurfacesQueries.H b/src/meshTools/searchableSurface/searchableSurfacesQueries.H index 789584dc7..819f4303e 100644 --- a/src/meshTools/searchableSurface/searchableSurfacesQueries.H +++ b/src/meshTools/searchableSurface/searchableSurfacesQueries.H @@ -44,6 +44,7 @@ namespace Foam // Forward declaration of classes class plane; +class pointConstraint; /*---------------------------------------------------------------------------*\ Class searchableSurfacesQueries Declaration @@ -51,68 +52,8 @@ class plane; class searchableSurfacesQueries { - // Private data - // Private Member Functions - //- Temporary wrapper around findNearest. Used in facesIntersection only - static pointIndexHit tempFindNearest - ( - const searchableSurface&, - const point& pt, - const scalar initDistSqr - ); - - //- Calculate sum of distances to nearest point on surfaces. Is used - // in minimisation to find intersection. Returns sum of (square of) - // distances to the surfaces. - static scalar sumDistSqr - ( - const PtrList&, - const labelList& surfacesToTest, - const scalar initialDistSqr, // search box - const point& pt - ); - - //- Takes the tet (points p) and reflects the point with the - // highest value around the centre (pSum). Checks if it gets closer - // and updates p, y if so. - static scalar tryMorphTet - ( - const PtrList&, - const labelList& surfacesToTest, - const scalar initialDistSqr, - List& p, - List& y, - vector& pSum, - const label ihi, - const scalar fac - ); - - //- Downhill simplex method: find the point with min cumulative - // distance to all surfaces. Does so by morphing a tet (points p). - // Returns the point on the 0th surface or hit if not reached within - // maxIters iterations. - static bool morphTet - ( - const PtrList&, - const labelList& surfacesToTest, - const scalar initialDistSqr, - const scalar convergenceDistSqr, - const label maxIter, - List& p, - List& y - ); - - //static void findAllIntersections - //( - // const searchableSurface& s, - // const pointField& start, - // const pointField& end, - // const vectorField& smallVec, - // List>& - //); - static void mergeHits ( const point& start, @@ -158,7 +99,7 @@ public: List>& surfaceHits ); - //Find intersections of edge nearest to both endpoints. + //- Find intersections of edge nearest to both endpoints. static void findNearestIntersection ( const PtrList& allSurfaces, @@ -195,6 +136,21 @@ public: List& nearestInfo ); + //- Find nearest points that are on all supplied surfaces + // (nearest point if single surface; nearest intersection by + // steepst descent if on multiple surfaces). Returns current + // best guess). Wip. + static void findNearest + ( + const PtrList& allSurfaces, + const labelList& surfacesToTest, + const pointField& start, + const scalarField& distSqr, + pointField& near, + List& constraint, + const label nIter = 20 + ); + //- Find signed distance to nearest surface. Outside is positive. // illegalHandling: how to handle non-inside or outside // OUTSIDE : treat as outside @@ -217,18 +173,6 @@ public: const PtrList& allSurfaces, const labelList& surfacesToTest ); - - // Single point queries - - //- Calculate point which is on a set of surfaces. WIP. - static pointIndexHit facesIntersection - ( - const PtrList& allSurfaces, - const labelList& surfacesToTest, - const scalar initDistSqr, - const scalar convergenceDistSqr, - const point& start - ); }; diff --git a/tutorials/mesh/blockMesh/sphere7ProjectedEdges/Allrun b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/Allrun new file mode 100755 index 000000000..3c4ef3e7f --- /dev/null +++ b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/Allrun @@ -0,0 +1,9 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # Run from this directory + +# Source tutorial run functions +. $WM_PROJECT_DIR/bin/tools/RunFunctions + +runApplication blockMesh + +#------------------------------------------------------------------------------ diff --git a/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict new file mode 100644 index 000000000..9706d2149 --- /dev/null +++ b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/blockMeshDict @@ -0,0 +1,143 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + object blockMeshDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +verbose no; + +geometry +{ + sphere + { + type searchableSphere; + centre (0 0 0); + radius 1; + } + + innerSphere + { + type searchableSphere; + centre (0 0 0); + radius 0.5; + } +} + +scale 1; + +n 10; + +// Rough approximation of corner points +v 0.5773502; +mv -0.5773502; +vh 0.2886751; +mvh -0.2886751; + +vertices +( + // Inner block points + project ($mvh $mvh $mvh) (innerSphere) + project ( $vh $mvh $mvh) (innerSphere) + project ( $vh $vh $mvh) (innerSphere) + project ($mvh $vh $mvh) (innerSphere) + project ($mvh $mvh $vh) (innerSphere) + project ( $vh $mvh $vh) (innerSphere) + project ( $vh $vh $vh) (innerSphere) + project ($mvh $vh $vh) (innerSphere) + + // Outer block points + project ($mv $mv $mv) (sphere) + project ( $v $mv $mv) (sphere) + project ( $v $v $mv) (sphere) + project ($mv $v $mv) (sphere) + project ($mv $mv $v) (sphere) + project ( $v $mv $v) (sphere) + project ( $v $v $v) (sphere) + project ($mv $v $v) (sphere) +); + +blocks +( + hex ( 0 1 2 3 4 5 6 7) ($n $n $n) simpleGrading (1 1 1) + hex ( 9 8 12 13 1 0 4 5) ($n $n $n) simpleGrading (1 1 1) + hex (10 9 13 14 2 1 5 6) ($n $n $n) simpleGrading (1 1 1) + hex (11 10 14 15 3 2 6 7) ($n $n $n) simpleGrading (1 1 1) + hex ( 8 11 15 12 0 3 7 4) ($n $n $n) simpleGrading (1 1 1) + hex ( 8 9 10 11 0 1 2 3) ($n $n $n) simpleGrading (1 1 1) + hex (13 12 15 14 5 4 7 6) ($n $n $n) simpleGrading (1 1 1) +); + +edges +( + // Outer edges + project 8 9 (sphere) + project 10 11 (sphere) + project 14 15 (sphere) + project 12 13 (sphere) + + project 8 11 (sphere) + project 9 10 (sphere) + project 13 14 (sphere) + project 12 15 (sphere) + + project 8 12 (sphere) + project 9 13 (sphere) + project 10 14 (sphere) + project 11 15 (sphere) + + + // Inner edges + project 0 1 (innerSphere) + project 2 3 (innerSphere) + project 6 7 (innerSphere) + project 4 5 (innerSphere) + + project 0 3 (innerSphere) + project 1 2 (innerSphere) + project 5 6 (innerSphere) + project 4 7 (innerSphere) + + project 0 4 (innerSphere) + project 1 5 (innerSphere) + project 2 6 (innerSphere) + project 3 7 (innerSphere) +); + +faces +( + project ( 8 12 15 11) sphere + project (10 14 13 9) sphere + project ( 9 13 12 8) sphere + project (11 15 14 10) sphere + project ( 8 11 10 9) sphere + project (12 13 14 15) sphere +); + +boundary +( + walls + { + type wall; + faces + ( + ( 8 12 15 11) + (10 14 13 9) + ( 9 13 12 8) + (11 15 14 10) + ( 8 11 10 9) + (12 13 14 15) + ); + } +); + +// ************************************************************************* // diff --git a/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/controlDict b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/controlDict new file mode 100644 index 000000000..6b6304814 --- /dev/null +++ b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/controlDict @@ -0,0 +1,49 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object controlDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +application blockMesh; + +startFrom startTime; + +startTime 0; + +stopAt endTime; + +endTime 0; + +deltaT 0; + +writeControl timeStep; + +writeInterval 1; + +purgeWrite 0; + +writeFormat ascii; + +writePrecision 6; + +writeCompression off; + +timeFormat general; + +timePrecision 6; + +runTimeModifiable true; + + +// ************************************************************************* // diff --git a/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/fvSchemes b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/fvSchemes new file mode 100644 index 000000000..606b28ce0 --- /dev/null +++ b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/fvSchemes @@ -0,0 +1,37 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSchemes; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +ddtSchemes +{} + +gradSchemes +{} + +divSchemes +{} + +laplacianSchemes +{} + +interpolationSchemes +{} + +snGradSchemes +{} + + +// ************************************************************************* // diff --git a/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/fvSolution b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/fvSolution new file mode 100644 index 000000000..6db4e1a6b --- /dev/null +++ b/tutorials/mesh/blockMesh/sphere7ProjectedEdges/system/fvSolution @@ -0,0 +1,18 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "system"; + object fvSolution; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +// ************************************************************************* //