diff --git a/Allwclean b/Allwclean
index fa17a8986..1f55d5690 100755
--- a/Allwclean
+++ b/Allwclean
@@ -16,7 +16,9 @@ wclean ./src/dfCombustionModels
wclean ./src/TurbulenceModels/compressible
wclean ./src/TurbulenceModels/turbulenceModels
wclean ./src/dfCanteraMixture
+wclean ./src/dfMeshTools
wclean ./applications/utilities/flameSpeed
+wclean ./applications/utilities/dfAutoPatch
wclean ./applications/solvers/dfSprayFoam
wclean ./applications/solvers/df0DFoam
wclean ./applications/solvers/dfLowMachFoam
diff --git a/Allwmake b/Allwmake
index dd7ff8071..74136b067 100755
--- a/Allwmake
+++ b/Allwmake
@@ -19,6 +19,7 @@ wmake src/dynamicMesh
wmake src/dynamicFvMesh
wmake src/fluxSchemes
wmake src/radiationModels
+wmake src/dfMeshTools
wmake applications/solvers/df0DFoam
wmake applications/solvers/dfLowMachFoam
@@ -28,3 +29,4 @@ wmake applications/solvers/dfBuoyancyFoam
wmake applications/solvers/dfSteadyFoam
wmake applications/utilities/flameSpeed
+wmake applications/utilities/dfAutoPatch
diff --git a/applications/utilities/dfAutoPatch/Make/files b/applications/utilities/dfAutoPatch/Make/files
new file mode 100644
index 000000000..9483fe7b2
--- /dev/null
+++ b/applications/utilities/dfAutoPatch/Make/files
@@ -0,0 +1,3 @@
+dfAutoPatch.C
+
+EXE = $(DF_APPBIN)/dfAutoPatch
\ No newline at end of file
diff --git a/applications/utilities/dfAutoPatch/Make/options b/applications/utilities/dfAutoPatch/Make/options
new file mode 100644
index 000000000..7349856ca
--- /dev/null
+++ b/applications/utilities/dfAutoPatch/Make/options
@@ -0,0 +1,7 @@
+EXE_INC = \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/dynamicMesh/lnInclude
+
+EXE_LIBS = \
+ -ldynamicMesh \
+ -lmeshTools
diff --git a/applications/utilities/dfAutoPatch/dfAutoPatch.C b/applications/utilities/dfAutoPatch/dfAutoPatch.C
new file mode 100644
index 000000000..50903401f
--- /dev/null
+++ b/applications/utilities/dfAutoPatch/dfAutoPatch.C
@@ -0,0 +1,357 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2011-2018 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 .
+
+Application
+ dfAutoPatch
+
+Description
+ An extension of autoPatch to work with polyhedral meshes.
+ This utility divides external faces into patches based on (user supplied) feature
+ angle and the source patches can be limited to user supplied patches if desired.
+ The resulting patches are renamed to _autoN if only a subset of
+ patches is processed, otherwise to autoN.
+
+ Example usage:
+ dfAutoPatch -onlyPatches "(box)" 45
+
+Authorship
+ Author: Teng Zhang @ AISI,
+ Date: 2025-10-22
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "polyMesh.H"
+#include "Time.H"
+#include "boundaryMesh.H"
+#include "repatchPolyTopoChanger.H"
+#include "unitConversion.H"
+#include "OFstream.H"
+#include "ListOps.H"
+#include "wordList.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// Get all feature edges.
+void collectFeatureEdges(const boundaryMesh& bMesh, labelList& markedEdges)
+{
+ markedEdges.setSize(bMesh.mesh().nEdges());
+
+ label markedI = 0;
+
+ forAll(bMesh.featureSegments(), i)
+ {
+ const labelList& segment = bMesh.featureSegments()[i];
+
+ forAll(segment, j)
+ {
+ label featEdgeI = segment[j];
+
+ label meshEdgeI = bMesh.featureToEdge()[featEdgeI];
+
+ markedEdges[markedI++] = meshEdgeI;
+ }
+ }
+ markedEdges.setSize(markedI);
+}
+
+
+
+int main(int argc, char *argv[])
+{
+ #include "addOverwriteOption.H"
+ argList::noParallel();
+ argList::validArgs.append("feature angle[0-180]");
+ argList::addOption
+ (
+ "onlyPatches",
+ "wordList of boundary patch names to be repatched"
+ );
+
+ #include "setRootCase.H"
+ #include "createTime.H"
+ runTime.functionObjects().off();
+ #include "createPolyMesh.H"
+ const word oldInstance = mesh.pointsInstance();
+
+ wordList onlyPatches;
+
+ if (args.optionFound("onlyPatches"))
+ {
+ onlyPatches = wordList(args.optionLookup("onlyPatches")());
+
+ Info<< "Limiting autoPatch to patches: " << onlyPatches
+ << nl << endl;
+ }
+
+ Info<< "Mesh read in = "
+ << runTime.cpuTimeIncrement()
+ << " s\n" << endl << endl;
+
+
+ const scalar featureAngle = args.argRead(1);
+ const bool overwrite = args.optionFound("overwrite");
+
+ const scalar minCos = Foam::cos(degToRad(featureAngle));
+
+ Info<< "Feature:" << featureAngle << endl
+ << "minCos :" << minCos << endl
+ << endl;
+
+ //
+ // Use boundaryMesh to reuse all the featureEdge stuff in there.
+ //
+
+ boundaryMesh bMesh;
+ bMesh.read(mesh);
+
+ // Set feature angle (calculate feature edges)
+ bMesh.setFeatureEdges(minCos);
+
+ // Collect all feature edges as edge labels
+ labelList markedEdges;
+
+ collectFeatureEdges(bMesh, markedEdges);
+
+
+
+ const labelList& boundaryFaces = bMesh.meshFace();
+ const polyBoundaryMesh& originalBoundary = mesh.boundaryMesh();
+ const label nBoundaryFaces = boundaryFaces.size();
+
+ // (new) patch ID for every face in mesh.
+ labelList patchIDs(nBoundaryFaces, -1);
+ boolList processFace(nBoundaryFaces, false);
+ labelList originalPatchIDs(nBoundaryFaces, -1);
+
+ // Track per-patch suffix when only a subset is processed.
+ labelList patchSuffixCount(originalBoundary.size(), 0);
+
+ const bool limitToPatches = onlyPatches.size() != 0;
+
+ // Determine which faces are allowed to be repatched.
+ forAll(boundaryFaces, facei)
+ {
+ const label meshFacei = boundaryFaces[facei];
+ const label patchi = originalBoundary.whichPatch(meshFacei);
+ const word& patchName = originalBoundary[patchi].name();
+
+ originalPatchIDs[facei] = patchi;
+ patchIDs[facei] = patchi;
+
+ if (!limitToPatches || findIndex(onlyPatches, patchName) != -1)
+ {
+ patchIDs[facei] = -1;
+ processFace[facei] = true;
+ }
+ }
+
+ if (limitToPatches)
+ {
+ label selectedFaceCount = 0;
+
+ forAll(processFace, facei)
+ {
+ if (processFace[facei])
+ {
+ ++selectedFaceCount;
+ }
+ }
+
+ if (!selectedFaceCount)
+ {
+ Info<< "No faces found for the supplied onlyPatches option." << nl
+ << endl;
+ }
+ }
+
+ //
+ // Fill patchIDs with values for every face by floodfilling without
+ // crossing feature edge.
+ //
+
+ // Current patch number.
+ label newPatchi = bMesh.patches().size();
+
+ label suffix = 0;
+
+ while (true)
+ {
+ // Find first unset face respecting the optional patch filter.
+ label unsetFacei = -1;
+
+ forAll(patchIDs, facei)
+ {
+ if (patchIDs[facei] == -1 && processFace[facei])
+ {
+ unsetFacei = facei;
+ break;
+ }
+ }
+
+ if (unsetFacei == -1)
+ {
+ // All faces have patchID set. Exit.
+ break;
+ }
+
+ // Found unset face. Create patch for it.
+ word patchName;
+ if (limitToPatches)
+ {
+ const label basePatchi = originalPatchIDs[unsetFacei];
+ const word& baseName = originalBoundary[basePatchi].name();
+ label& nextSuffix = patchSuffixCount[basePatchi];
+
+ do
+ {
+ patchName = baseName + "_auto" + name(nextSuffix++);
+ }
+ while (bMesh.findPatchID(patchName) != -1);
+ }
+ else
+ {
+ do
+ {
+ patchName = "auto" + name(suffix++);
+ }
+ while (bMesh.findPatchID(patchName) != -1);
+ }
+
+ bMesh.addPatch(patchName);
+
+ bMesh.changePatchType(patchName, "patch");
+
+
+ // Fill visited with all faces reachable from unsetFacei.
+ boolList visited(nBoundaryFaces, false);
+
+ bMesh.markFaces(markedEdges, unsetFacei, visited);
+
+
+ // Assign all visited faces to current patch
+ label nVisited = 0;
+
+ forAll(visited, facei)
+ {
+ if (visited[facei] && processFace[facei])
+ {
+ nVisited++;
+
+ patchIDs[facei] = newPatchi;
+ }
+ }
+
+ Info<< "Assigned " << nVisited << " faces to patch " << patchName
+ << endl << endl;
+
+ newPatchi++;
+ }
+
+
+
+ const PtrList& patches = bMesh.patches();
+
+ // Create new list of patches with old ones first
+ List newPatchPtrList(patches.size());
+
+ newPatchi = 0;
+
+ // Copy old patches
+ forAll(mesh.boundaryMesh(), patchi)
+ {
+ const polyPatch& patch = mesh.boundaryMesh()[patchi];
+
+ newPatchPtrList[newPatchi] =
+ patch.clone
+ (
+ mesh.boundaryMesh(),
+ newPatchi,
+ patch.size(),
+ patch.start()
+ ).ptr();
+
+ newPatchi++;
+ }
+
+ // Add new ones with empty size.
+ for (label patchi = newPatchi; patchi < patches.size(); patchi++)
+ {
+ const boundaryPatch& bp = patches[patchi];
+
+ newPatchPtrList[newPatchi] = polyPatch::New
+ (
+ polyPatch::typeName,
+ bp.name(),
+ 0,
+ mesh.nFaces(),
+ newPatchi,
+ mesh.boundaryMesh()
+ ).ptr();
+
+ newPatchi++;
+ }
+
+ if (!overwrite)
+ {
+ runTime++;
+ }
+
+
+ // Change patches
+ repatchPolyTopoChanger polyMeshRepatcher(mesh);
+ polyMeshRepatcher.changePatches(newPatchPtrList);
+
+
+ // Change face ordering
+
+ // Since bMesh read from mesh there is one to one mapping so we don't
+ // have to do the geometric stuff.
+ forAll(patchIDs, facei)
+ {
+ label meshFacei = boundaryFaces[facei];
+
+ polyMeshRepatcher.changePatchID(meshFacei, patchIDs[facei]);
+ }
+
+ polyMeshRepatcher.repatch();
+
+ // Write resulting mesh
+ if (overwrite)
+ {
+ mesh.setInstance(oldInstance);
+ }
+
+ // Set the precision of the points data to 10
+ IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+
+ mesh.write();
+
+ Info<< "End\n" << endl;
+
+ return 0;
+}
+
+
+// ************************************************************************* //
diff --git a/src/dfMeshTools/Make/files b/src/dfMeshTools/Make/files
new file mode 100644
index 000000000..ee06a14cb
--- /dev/null
+++ b/src/dfMeshTools/Make/files
@@ -0,0 +1,5 @@
+cellSources/cylinderCutToCell/cylinderCutToCell.C
+
+searchableSurfaces/searchableCone/searchableCone.C
+
+LIB = $(DF_LIBBIN)/libdfMeshTools
diff --git a/src/dfMeshTools/Make/options b/src/dfMeshTools/Make/options
new file mode 100644
index 000000000..28f288385
--- /dev/null
+++ b/src/dfMeshTools/Make/options
@@ -0,0 +1,11 @@
+EXE_INC = \
+ -I$(LIB_SRC)/meshTools/lnInclude \
+ -I$(LIB_SRC)/triSurface/lnInclude \
+ -I$(LIB_SRC)/surfMesh/lnInclude \
+ -I$(LIB_SRC)/fileFormats/lnInclude
+
+LIB_LIBS = \
+ -lmeshTools \
+ -ltriSurface \
+ -lsurfMesh \
+ -lfileFormats
diff --git a/src/dfMeshTools/README.md b/src/dfMeshTools/README.md
new file mode 100644
index 000000000..03d4a9385
--- /dev/null
+++ b/src/dfMeshTools/README.md
@@ -0,0 +1,50 @@
+dfMeshTools
+===========
+A collection of tools for mesh generation and manipulation, developed as part of the DeepFlame project.
+
+## List of tools
+### cellSources
+1. cylinderCutToCell: A tool to cut a cylindrical region from a cell in the mesh
+### searchableSurfaces
+2. searchableCone: A tool to define a searchable cone surface for mesh generation. Typically used for defining ROI (region of interest) in the mesh where finer resolution is needed.
+
+## Usage
+**NOTE**: In `controlDict`, the library should be added to `libs`:
+```cpp
+libs ( "libdfMeshTools.so" );
+```
+
+1. For `cellSources`, the tool can be used in `topoSetDict` as follows:
+ ```cpp
+ actions
+ (
+
+ {
+ name cylinderCut;
+ type cellSet;
+ action new;
+ source cylinderCutToCell;
+ sourceInfo
+ {
+ p1 (0 0 0.015);
+ p2 (0 0 0.025);
+ radius 0.0296;
+ }
+ }
+ );
+ ```
+
+2. For `searchableSurfaces`, the tool can be used in `snappyHexMeshDict` as follows:
+ ```cpp
+ geometry
+ {
+ refineRegion
+ {
+ type searchableCone;
+ point1 (0.0 0 0);
+ point2 (0.015 0 0);
+ radius1 0.005;
+ radius2 0.007;
+ }
+ }
+ ```
\ No newline at end of file
diff --git a/src/dfMeshTools/cellSources/cylinderCutToCell/cylinderCutToCell.C b/src/dfMeshTools/cellSources/cylinderCutToCell/cylinderCutToCell.C
new file mode 100644
index 000000000..8bed45e12
--- /dev/null
+++ b/src/dfMeshTools/cellSources/cylinderCutToCell/cylinderCutToCell.C
@@ -0,0 +1,288 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2011-2018 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 "cylinderCutToCell.H"
+#include "polyMesh.H"
+#include "addToRunTimeSelectionTable.H"
+#include "cellModel.H"
+#include "boundBox.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(cylinderCutToCell, 0);
+ addToRunTimeSelectionTable(topoSetSource, cylinderCutToCell, word);
+ addToRunTimeSelectionTable(topoSetSource, cylinderCutToCell, istream);
+}
+
+
+Foam::topoSetSource::addToUsageTable Foam::cylinderCutToCell::usage_
+(
+ cylinderCutToCell::typeName,
+ "\n Usage: cylinderCutToCell (p1X p1Y p1Z) (p2X p2Y p2Z) radius\n\n"
+ " Select all cells that intersect with the bounding cylinder\n\n"
+);
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+// Helper function: check if a point is inside the cylinder
+bool Foam::cylinderCutToCell::pointInCylinder
+(
+ const point& pt,
+ const vector& axis,
+ const scalar magAxis2,
+ const scalar rad2
+) const
+{
+ vector d = pt - p1_;
+ scalar magD = d & axis;
+
+ if ((magD >= 0) && (magD <= magAxis2))
+ {
+ scalar d2 = (d & d) - sqr(magD)/magAxis2;
+ if (d2 <= rad2)
+ {
+ return true;
+ }
+ }
+ return false;
+}
+
+
+// Helper function: check if cylinder intersects with a bounding box
+// Returns true only if the cylinder (not a capsule) intersects the box
+bool Foam::cylinderCutToCell::cylinderIntersectsBoundBox
+(
+ const boundBox& bb,
+ const vector& axis,
+ const vector& axisDir,
+ const scalar magAxis,
+ const scalar rad2
+) const
+{
+ // Check if the cylinder axis passes through (or close to) the bounding box
+ // We need to ensure the intersection is within the cylinder length (not beyond ends)
+
+ const point& boxMin = bb.min();
+ const point& boxMax = bb.max();
+ const point boxCenter = bb.midpoint();
+ const scalar boxRadius = mag(boxMax - boxMin) * 0.5;
+
+ // Find projection of box center onto cylinder axis
+ vector d = boxCenter - p1_;
+ scalar t = (d & axisDir); // projection along normalized axis
+
+ // Check if the box is completely outside the cylinder's axial range
+ // Consider the box could extend into the cylinder range by boxRadius
+ if (t < -boxRadius || t > magAxis + boxRadius)
+ {
+ return false;
+ }
+
+ // Clamp t to [0, magAxis] to find closest point on axis within cylinder
+ scalar tClamped = max(scalar(0), min(t, magAxis));
+
+ // Closest point on axis (within cylinder length)
+ point closestOnAxis = p1_ + tClamped * axisDir;
+
+ // Distance from closest point to box center
+ scalar dist2 = magSqr(closestOnAxis - boxCenter);
+
+ // If distance is less than radius + box diagonal/2, they might intersect
+ if (dist2 <= sqr(radius_ + boxRadius))
+ {
+ return true;
+ }
+
+ return false;
+}
+
+
+void Foam::cylinderCutToCell::combine(topoSet& set, const bool add) const
+{
+ const vector axis = p2_ - p1_;
+ const scalar rad2 = sqr(radius_);
+ const scalar magAxis2 = magSqr(axis);
+ const scalar magAxis = Foam::sqrt(magAxis2);
+ const vector axisDir = axis / (magAxis + VSMALL);
+
+ const pointField& pts = mesh_.points();
+ const cellList& cells = mesh_.cells();
+ const faceList& faces = mesh_.faces();
+
+ forAll(cells, celli)
+ {
+ bool selected = false;
+
+ // First, quick check: does any vertex of the cell lie inside the cylinder?
+ const cell& c = cells[celli];
+ labelHashSet cellPoints(c.size() * 4);
+
+ // Collect all unique points of this cell
+ forAll(c, facei)
+ {
+ const face& f = faces[c[facei]];
+ forAll(f, pointi)
+ {
+ cellPoints.insert(f[pointi]);
+ }
+ }
+
+ // Check if any vertex is inside the cylinder
+ forAllConstIter(labelHashSet, cellPoints, iter)
+ {
+ if (pointInCylinder(pts[iter.key()], axis, magAxis2, rad2))
+ {
+ selected = true;
+ break;
+ }
+ }
+
+ // If no vertex inside, check if cylinder axis intersects cell bounding box
+ if (!selected)
+ {
+ // Build bounding box for this cell - collect points first
+ pointField cellPts(cellPoints.size());
+ label pi = 0;
+ forAllConstIter(labelHashSet, cellPoints, iter)
+ {
+ cellPts[pi++] = pts[iter.key()];
+ }
+ boundBox cellBB(cellPts, false);
+
+ if (cylinderIntersectsBoundBox(cellBB, axis, axisDir, magAxis, rad2))
+ {
+ // More detailed check: find the closest point on cylinder axis
+ // to the cell center and verify it's within cylinder bounds
+ const point& cellCtr = mesh_.cellCentres()[celli];
+ vector d = cellCtr - p1_;
+ scalar t = (d & axisDir);
+
+ // Only proceed if cell center projection is within cylinder length
+ // (with some tolerance based on cell size)
+ scalar cellRadius = mag(cellBB.max() - cellBB.min()) * 0.5;
+ if (t >= -cellRadius && t <= magAxis + cellRadius)
+ {
+ // Clamp to cylinder range
+ scalar tClamped = max(scalar(0), min(t, magAxis));
+ point closestOnAxis = p1_ + tClamped * axisDir;
+
+ // Check radial distance from axis to cell center
+ scalar radialDist2 = magSqr(cellCtr - closestOnAxis);
+
+ // Select if within radius + cell size tolerance
+ if (radialDist2 <= sqr(radius_ + cellRadius))
+ {
+ selected = true;
+ }
+ }
+ }
+ }
+
+ if (selected)
+ {
+ addOrDelete(set, celli, add);
+ }
+ }
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::cylinderCutToCell::cylinderCutToCell
+(
+ const polyMesh& mesh,
+ const vector& p1,
+ const vector& p2,
+ const scalar radius
+)
+:
+ topoSetSource(mesh),
+ p1_(p1),
+ p2_(p2),
+ radius_(radius)
+{}
+
+
+Foam::cylinderCutToCell::cylinderCutToCell
+(
+ const polyMesh& mesh,
+ const dictionary& dict
+)
+:
+ topoSetSource(mesh),
+ p1_(dict.lookup("p1")),
+ p2_(dict.lookup("p2")),
+ radius_(readScalar(dict.lookup("radius")))
+{}
+
+
+Foam::cylinderCutToCell::cylinderCutToCell
+(
+ const polyMesh& mesh,
+ Istream& is
+)
+:
+ topoSetSource(mesh),
+ p1_(checkIs(is)),
+ p2_(checkIs(is)),
+ radius_(readScalar(checkIs(is)))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
+
+Foam::cylinderCutToCell::~cylinderCutToCell()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+void Foam::cylinderCutToCell::applyToSet
+(
+ const topoSetSource::setAction action,
+ topoSet& set
+) const
+{
+ if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD))
+ {
+ Info<< " Adding cells intersecting with cylinder, with p1 = "
+ << p1_ << ", p2 = " << p2_ << " and radius = " << radius_ << endl;
+
+ combine(set, true);
+ }
+ else if (action == topoSetSource::DELETE)
+ {
+ Info<< " Removing cells intersecting with cylinder, with p1 = "
+ << p1_ << ", p2 = " << p2_ << " and radius = " << radius_ << endl;
+
+ combine(set, false);
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/dfMeshTools/cellSources/cylinderCutToCell/cylinderCutToCell.H b/src/dfMeshTools/cellSources/cylinderCutToCell/cylinderCutToCell.H
new file mode 100644
index 000000000..e5d34f735
--- /dev/null
+++ b/src/dfMeshTools/cellSources/cylinderCutToCell/cylinderCutToCell.H
@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration | Website: https://openfoam.org
+ \\ / A nd | Copyright (C) 2011-2019 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::cylinderCutToCell
+
+Description
+ A topoSetSource to select cells that intersect with a cylinder.
+ A cell is selected if any of its vertices is inside the cylinder,
+ or if the cylinder axis passes through the cell's bounding box.
+
+ Example usage in topoSetDict:
+ \verbatim
+ {
+ name cylinderCut;
+ type cellSet;
+ action new;
+ source cylinderCutToCell;
+ sourceInfo
+ {
+ p1 (0 0 0.015);
+ p2 (0 0 0.025);
+ radius 0.0275;
+ }
+ }
+ \endverbatim
+
+SourceFiles
+ cylinderToCell.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef cylinderCutToCell_H
+#define cylinderCutToCell_H
+
+#include "topoSetSource.H"
+#include "boundBox.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class cylinderCutToCell Declaration
+\*---------------------------------------------------------------------------*/
+
+class cylinderCutToCell
+:
+ public topoSetSource
+{
+
+ // Private Data
+
+ //- Add usage string
+ static addToUsageTable usage_;
+
+ //- First point on cylinder axis
+ vector p1_;
+
+ //- Second point on cylinder axis
+ vector p2_;
+
+ //- Radius
+ scalar radius_;
+
+
+ // Private Member Functions
+
+ //- Check if a point is inside the cylinder
+ bool pointInCylinder
+ (
+ const point& pt,
+ const vector& axis,
+ const scalar magAxis2,
+ const scalar rad2
+ ) const;
+
+ //- Check if cylinder intersects with a bounding box
+ bool cylinderIntersectsBoundBox
+ (
+ const boundBox& bb,
+ const vector& axis,
+ const vector& axisDir,
+ const scalar magAxis,
+ const scalar rad2
+ ) const;
+
+ void combine(topoSet& set, const bool add) const;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("cylinderCutToCell");
+
+
+ // Constructors
+
+ //- Construct from components
+ cylinderCutToCell
+ (
+ const polyMesh& mesh,
+ const vector& p1,
+ const vector& p2,
+ const scalar radius
+ );
+
+ //- Construct from dictionary
+ cylinderCutToCell
+ (
+ const polyMesh& mesh,
+ const dictionary& dict
+ );
+
+ //- Construct from Istream
+ cylinderCutToCell
+ (
+ const polyMesh& mesh,
+ Istream&
+ );
+
+
+ //- Destructor
+ virtual ~cylinderCutToCell();
+
+
+ // Member Functions
+
+ virtual sourceType setType() const
+ {
+ return CELLSETSOURCE;
+ }
+
+ virtual void applyToSet
+ (
+ const topoSetSource::setAction action,
+ topoSet&
+ ) const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
diff --git a/src/dfMeshTools/searchableSurfaces/searchableCone/searchableCone.C b/src/dfMeshTools/searchableSurfaces/searchableCone/searchableCone.C
new file mode 100644
index 000000000..c55815b56
--- /dev/null
+++ b/src/dfMeshTools/searchableSurfaces/searchableCone/searchableCone.C
@@ -0,0 +1,1109 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2015-2020 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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 "searchableCone.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+ defineTypeNameAndDebug(searchableCone, 0);
+ addToRunTimeSelectionTable
+ (
+ searchableSurface,
+ searchableCone,
+ dict
+ );
+ addNamedToRunTimeSelectionTable
+ (
+ searchableSurface,
+ searchableCone,
+ dict,
+ cone
+ );
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
+
+Foam::tmp Foam::searchableCone::coordinates() const
+{
+ tmp tCtrs(new pointField(1, 0.5*(point1_ + point2_)));
+
+ return tCtrs;
+}
+
+
+void Foam::searchableCone::boundingSpheres
+(
+ pointField& centres,
+ scalarField& radiusSqr
+) const
+{
+ centres.setSize(1);
+ centres[0] = 0.5*(point1_ + point2_);
+
+ radiusSqr.setSize(1);
+ if (radius1_ > radius2_)
+ {
+ radiusSqr[0] = Foam::magSqr(point1_-centres[0]) + Foam::sqr(radius1_);
+ }
+ else
+ {
+ radiusSqr[0] = Foam::magSqr(point2_-centres[0]) + Foam::sqr(radius2_);
+ }
+
+ // Add a bit to make sure all points are tested inside
+ radiusSqr += Foam::sqr(SMALL);
+}
+
+
+Foam::tmp Foam::searchableCone::points() const
+{
+ tmp tPts(new pointField(2));
+ pointField& pts = tPts.ref();
+
+ pts[0] = point1_;
+ pts[1] = point2_;
+
+ return tPts;
+}
+
+
+void Foam::searchableCone::findNearestAndNormal
+(
+ const point& sample,
+ const scalar nearestDistSqr,
+ pointIndexHit& info,
+ vector& nearNormal
+) const
+{
+ vector v(sample - point1_);
+
+ // Decompose sample-point1 into normal and parallel component
+ const scalar parallel = (v & unitDir_);
+
+ // Remove the parallel component and normalise
+ v -= parallel*unitDir_;
+
+ const scalar magV = mag(v);
+ v = v/mag(v); //this is based on the Vector primitive in v2006
+
+ // Nearest and normal on disk at point1
+ point disk1Point(point1_ + min(max(magV, innerRadius1_), radius1_)*v);
+ vector disk1Normal(-unitDir_);
+
+ // Nearest and normal on disk at point2
+ point disk2Point(point2_ + min(max(magV, innerRadius2_), radius2_)*v);
+ vector disk2Normal(unitDir_);
+
+ // Nearest and normal on cone. Initialise to far-away point so if not
+ // set it picks one of the disk points
+ point nearCone(point::uniform(GREAT));
+ vector normalCone(1, 0, 0);
+
+ // Nearest and normal on inner cone. Initialise to far-away point
+ point iCnearCone(point::uniform(GREAT));
+ vector iCnormalCone(1, 0, 0);
+
+ point projPt1 = point1_+ radius1_*v;
+ point projPt2 = point2_+ radius2_*v;
+
+ vector p1 = (projPt1 - point1_);
+ if (mag(p1) > ROOTVSMALL)
+ {
+ p1 /= mag(p1);
+
+ // Find vector along the two end of cone
+ const vector b = normalised(projPt2 - projPt1);
+
+ // Find the vector along sample pt and pt at one end of cone
+ vector a = (sample - projPt1);
+
+ if (mag(a) <= ROOTVSMALL)
+ {
+ // Exception: sample on disk1. Redo with projPt2.
+ a = (sample - projPt2);
+
+ // Find normal unitvector
+ nearCone = (a & b)*b + projPt2;
+
+ vector b1 = (p1 & b)*b;
+ normalCone = normalised(p1 - b1);
+ }
+ else
+ {
+ // Find nearest point on cone surface
+ nearCone = (a & b)*b + projPt1;
+
+ // Find projection along surface of cone
+ vector b1 = (p1 & b)*b;
+ normalCone = normalised(p1 - b1);
+ }
+
+ if (innerRadius1_ > 0 || innerRadius2_ > 0)
+ {
+ // Same for inner radius
+ point iCprojPt1 = point1_+ innerRadius1_*v;
+ point iCprojPt2 = point2_+ innerRadius2_*v;
+
+ const vector iCp1 = normalised(iCprojPt1 - point1_);
+
+ // Find vector along the two end of cone
+ const vector iCb = normalised(iCprojPt2 - iCprojPt1);
+
+
+ // Find the vector along sample pt and pt at one end of conde
+ vector iCa(sample - iCprojPt1);
+
+ if (mag(iCa) <= ROOTVSMALL)
+ {
+ iCa = (sample - iCprojPt2);
+
+ // Find normal unitvector
+ iCnearCone = (iCa & iCb)*iCb+iCprojPt2;
+
+ vector b1 = (iCp1 & iCb)*iCb;
+ iCnormalCone = normalised(iCp1 - b1);
+ }
+ else
+ {
+ // Find nearest point on cone surface
+ iCnearCone = (iCa & iCb)*iCb+iCprojPt1;
+
+ // Find projection along surface of cone
+ vector b1 = (iCp1 & iCb)*iCb;
+ iCnormalCone = normalised(iCp1 - b1);
+ }
+ }
+ }
+
+
+ // Select nearest out of the 4 points (outer cone, disk1, disk2, inner
+ // cone)
+
+ FixedList dist;
+ dist[0] = magSqr(nearCone-sample);
+ dist[1] = magSqr(disk1Point-sample);
+ dist[2] = magSqr(disk2Point-sample);
+ dist[3] = magSqr(iCnearCone-sample);
+
+ const label minI = findMin(dist);
+
+
+ // Snap the point to the corresponding surface
+
+ if (minI == 0) // Outer cone
+ {
+ // Closest to (infinite) outer cone. See if needs clipping to end disks
+
+ {
+ vector v1(nearCone-point1_);
+ scalar para = (v1 & unitDir_);
+ // Remove the parallel component and normalise
+ v1 -= para*unitDir_;
+ const scalar magV1 = mag(v1);
+ v1 = v1/magV1;
+
+ if (para < 0.0 && magV1 >= radius1_)
+ {
+ // Near point 1. Set point to intersection of disk and cone.
+ // Keep normal from cone.
+ nearCone = disk1Point;
+ }
+ else if (para < 0.0 && magV1 < radius1_)
+ {
+ // On disk1
+ nearCone = disk1Point;
+ normalCone = disk1Normal;
+ }
+ else if (para > magDir_ && magV1 >= radius2_)
+ {
+ // Near point 2. Set point to intersection of disk and cone.
+ // Keep normal from cone.
+ nearCone = disk2Point;
+ }
+ else if (para > magDir_ && magV1 < radius2_)
+ {
+ // On disk2
+ nearCone = disk2Point;
+ normalCone = disk2Normal;
+ }
+ }
+ info.setPoint(nearCone);
+ nearNormal = normalCone;
+ }
+ else if (minI == 1) // Near to disk1
+ {
+ info.setPoint(disk1Point);
+ nearNormal = disk1Normal;
+ }
+ else if (minI == 2) // Near to disk2
+ {
+ info.setPoint(disk2Point);
+ nearNormal = disk2Normal;
+ }
+ else // Near to inner cone
+ {
+ {
+ vector v1(iCnearCone-point1_);
+ scalar para = (v1 & unitDir_);
+ // Remove the parallel component and normalise
+ v1 -= para*unitDir_;
+
+ const scalar magV1 = mag(v1);
+ v1 = v1/magV1;
+
+ if (para < 0.0 && magV1 >= innerRadius1_)
+ {
+ iCnearCone = disk1Point;
+ }
+ else if (para < 0.0 && magV1 < innerRadius1_)
+ {
+ iCnearCone = disk1Point;
+ iCnormalCone = disk1Normal;
+ }
+ else if (para > magDir_ && magV1 >= innerRadius2_)
+ {
+ iCnearCone = disk2Point;
+ }
+ else if (para > magDir_ && magV1 < innerRadius2_)
+ {
+ iCnearCone = disk2Point;
+ iCnormalCone = disk2Normal;
+ }
+ }
+ info.setPoint(iCnearCone);
+ nearNormal = iCnormalCone;
+ }
+
+
+ if (magSqr(sample - info.rawPoint()) < nearestDistSqr)
+ {
+ info.setHit();
+ info.setIndex(0);
+ }
+}
+
+
+Foam::scalar Foam::searchableCone::radius2
+(
+ const searchableCone& cone,
+ const point& pt
+)
+{
+ const vector x = (pt-cone.point1_) ^ cone.unitDir_;
+ return x&x;
+}
+
+
+// From mrl.nyu.edu/~dzorin/rend05/lecture2.pdf,
+// Ray Tracing II, Infinite cone ray intersection.
+void Foam::searchableCone::findLineAll
+(
+ const searchableCone& cone,
+ const scalar innerRadius1,
+ const scalar innerRadius2,
+ const point& start,
+ const point& end,
+ pointIndexHit& near,
+ pointIndexHit& far
+) const
+{
+ near.setMiss();
+ far.setMiss();
+
+ vector point1Start(start-cone.point1_);
+ vector point2Start(start-cone.point2_);
+ vector point1End(end-cone.point1_);
+
+ // Quick rejection of complete vector outside endcaps
+ scalar s1 = point1Start & (cone.unitDir_);
+ scalar s2 = point1End & (cone.unitDir_);
+
+ if ((s1 < 0.0 && s2 < 0.0) || (s1 > cone.magDir_ && s2 > cone.magDir_))
+ {
+ return;
+ }
+
+ // Line as P = start+t*V where V is unit vector and t=[0..mag(end-start)]
+ vector V(end-start);
+ scalar magV = mag(V);
+ if (magV < ROOTVSMALL)
+ {
+ return;
+ }
+ V /= magV;
+
+
+ // We now get the nearest intersections to start. This can either be
+ // the intersection with the end plane or with the cylinder side.
+
+ // Get the two points (expressed in t) on the end planes. This is to
+ // clip any cylinder intersection against.
+ scalar tPoint1;
+ scalar tPoint2;
+
+ // Maintain the two intersections with the endcaps
+ scalar tNear = VGREAT;
+ scalar tFar = VGREAT;
+
+ scalar radius_sec = cone.radius1_;
+
+ {
+ // Find dot product: mag(s)>VSMALL suggest that it is greater
+ scalar s = (V & unitDir_);
+ if (mag(s) > VSMALL)
+ {
+ tPoint1 = -s1/s;
+ tPoint2 = -(point2Start&(cone.unitDir_))/s;
+
+ if (tPoint2 < tPoint1)
+ {
+ Swap(tPoint1, tPoint2);
+ }
+ if (tPoint1 > magV || tPoint2 < 0)
+ {
+ return;
+ }
+ // See if the points on the endcaps are actually inside the cylinder
+ if (tPoint1 >= 0 && tPoint1 <= magV)
+ {
+ scalar r2 = radius2(cone, start+tPoint1*V);
+ vector p1 = (start+tPoint1*V-point1_);
+ vector p2 = (start+tPoint1*V-point2_);
+ radius_sec = cone.radius1_;
+ scalar inC_radius_sec = innerRadius1_;
+
+ if (mag(p2&(cone.unitDir_)) < mag(p1&(cone.unitDir_)))
+ {
+ radius_sec = cone.radius2_;
+ inC_radius_sec = innerRadius2_;
+ }
+
+ if (r2 <= sqr(radius_sec) && r2 >= sqr(inC_radius_sec))
+ {
+ tNear = tPoint1;
+ }
+ }
+ if (tPoint2 >= 0 && tPoint2 <= magV)
+ {
+ // Decompose sample-point1 into normal and parallel component
+ vector p1 = (start+tPoint2*V-cone.point1_);
+ vector p2 = (start+tPoint2*V-cone.point2_);
+ radius_sec = cone.radius1_;
+ scalar inC_radius_sec = innerRadius1_;
+ if (mag(p2&(cone.unitDir_)) < mag(p1&(cone.unitDir_)))
+ {
+ radius_sec = cone.radius2_;
+ inC_radius_sec = innerRadius2_;
+ }
+ scalar r2 = radius2(cone, start+tPoint2*V);
+ if (r2 <= sqr(radius_sec) && r2 >= sqr(inC_radius_sec))
+ {
+ // Check if already have a near hit from point1
+ if (tNear <= 1)
+ {
+ tFar = tPoint2;
+ }
+ else
+ {
+ tNear = tPoint2;
+ }
+ }
+ }
+ }
+ else
+ {
+ // Vector perpendicular to cylinder. Check for outside already done
+ // above so just set tpoint to allow all.
+ tPoint1 = -VGREAT;
+ tPoint2 = VGREAT;
+ }
+ }
+
+
+ // Second order equation of the form a*t^2 + b*t + c
+ scalar a, b, c;
+
+ scalar deltaRadius = cone.radius2_-cone.radius1_;
+ if (mag(deltaRadius) <= ROOTVSMALL)
+ {
+ vector point1Start(start-cone.point1_);
+ const vector x = point1Start ^ cone.unitDir_;
+ const vector y = V ^ cone.unitDir_;
+ const scalar d = sqr(0.5*(cone.radius1_ + cone.radius2_));
+
+ a = (y&y);
+ b = 2*(x&y);
+ c = (x&x)-d;
+ }
+ else
+ {
+ vector va = cone.unitDir_;
+ vector v1 = normalised(end-start);
+ scalar p = (va&v1);
+ vector a1 = (v1-p*va);
+
+ // Determine the end point of the cone
+ point pa =
+ cone.unitDir_*cone.radius1_*mag(cone.point2_-cone.point1_)
+ /(-deltaRadius)
+ +cone.point1_;
+
+ scalar l2 = sqr(deltaRadius)+sqr(cone.magDir_);
+ scalar sqrCosAlpha = sqr(cone.magDir_)/l2;
+ scalar sqrSinAlpha = sqr(deltaRadius)/l2;
+
+
+ vector delP(start-pa);
+ vector p1 = (delP-(delP&va)*va);
+
+ a = sqrCosAlpha*((v1-p*va)&(v1-p*va))-sqrSinAlpha*p*p;
+ b =
+ 2.0*sqrCosAlpha*(a1&p1)
+ -2.0*sqrSinAlpha*(v1&va)*(delP&va);
+ c =
+ sqrCosAlpha
+ *((delP-(delP&va)*va)&(delP-(delP&va)*va))
+ -sqrSinAlpha*sqr(delP&va);
+ }
+
+ const scalar disc = b*b-4.0*a*c;
+
+ scalar t1 = -VGREAT;
+ scalar t2 = VGREAT;
+
+ if (disc < 0)
+ {
+ // Fully outside
+ return;
+ }
+ else if (disc < ROOTVSMALL)
+ {
+ // Single solution
+ if (mag(a) > ROOTVSMALL)
+ {
+ t1 = -b/(2.0*a);
+
+ if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
+ {
+ // Valid. Insert sorted.
+ if (t1 < tNear)
+ {
+ tFar = tNear;
+ tNear = t1;
+ }
+ else if (t1 < tFar)
+ {
+ tFar = t1;
+ }
+ }
+ else
+ {
+ return;
+ }
+ }
+ else
+ {
+ // Aligned with axis. Check if outside radius
+ if (c > 0.0)
+ {
+ return;
+ }
+ }
+ }
+ else
+ {
+ if (mag(a) > ROOTVSMALL)
+ {
+ scalar sqrtDisc = sqrt(disc);
+
+ t1 = (-b - sqrtDisc)/(2.0*a);
+ t2 = (-b + sqrtDisc)/(2.0*a);
+ if (t2 < t1)
+ {
+ Swap(t1, t2);
+ }
+
+ if (t1 >= 0.0 && t1 <= magV && t1 >= tPoint1 && t1 <= tPoint2)
+ {
+ // Valid. Insert sorted.
+ if (t1 < tNear)
+ {
+ tFar = tNear;
+ tNear = t1;
+ }
+ else if (t1 < tFar)
+ {
+ tFar = t1;
+ }
+ }
+ if (t2>=0 && t2 <= magV && t2 >= tPoint1 && t2 <= tPoint2)
+ {
+ // Valid. Insert sorted.
+ if (t2 < tNear)
+ {
+ tFar = tNear;
+ tNear = t2;
+ }
+ else if (t2 < tFar)
+ {
+ tFar = t2;
+ }
+ }
+ }
+ else
+ {
+ // Aligned with axis. Check if outside radius
+ if (c > 0.0)
+ {
+ return;
+ }
+ }
+ }
+
+ // Check tNear, tFar
+ if (tNear>=0 && tNear <= magV)
+ {
+ near.setPoint(start+tNear*V);
+ near.setHit();
+ near.setIndex(0);
+ if (tFar <= magV)
+ {
+ far.setPoint(start+tFar*V);
+ far.setHit();
+ far.setIndex(0);
+ }
+ }
+ else if (tFar>=0 && tFar <= magV)
+ {
+ near.setPoint(start+tFar*V);
+ near.setHit();
+ near.setIndex(0);
+ }
+}
+
+
+void Foam::searchableCone::insertHit
+(
+ const point& start,
+ const point& end,
+ List& info,
+ const pointIndexHit& hit
+) const
+{
+ scalar smallDistSqr = SMALL*magSqr(end-start);
+
+ scalar hitMagSqr = magSqr(hit.hitPoint()-start);
+
+ forAll(info, i)
+ {
+ scalar d2 = magSqr(info[i].hitPoint()-start);
+
+ if (d2 > hitMagSqr+smallDistSqr)
+ {
+ // Insert at i.
+ label sz = info.size();
+ info.setSize(sz+1);
+ for (label j = sz; j > i; --j)
+ {
+ info[j] = info[j-1];
+ }
+ info[i] = hit;
+ return;
+ }
+ else if (d2 > hitMagSqr-smallDistSqr)
+ {
+ // hit is same point as info[i].
+ return;
+ }
+ }
+ // Append
+ label sz = info.size();
+ info.setSize(sz+1);
+ info[sz] = hit;
+}
+
+
+Foam::boundBox Foam::searchableCone::calcBounds() const
+{
+ // Adapted from
+ // http://www.gamedev.net/community/forums
+ // /topic.asp?topic_id=338522&forum_id=20&gforum_id=0
+
+ // Let cylinder have end points A,B and radius r,
+
+ // Bounds in direction X (same for Y and Z) can be found as:
+ // Let A.X= radius1_)
+ {
+ kr *= radius2_;
+ }
+ else
+ {
+ kr *= radius1_;
+ }
+
+ point min = point1_ - kr;
+ point max = point1_ + kr;
+
+ min = ::Foam::min(min, point2_ - kr);
+ max = ::Foam::max(max, point2_ + kr);
+
+ return boundBox(min, max);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
+
+Foam::searchableCone::searchableCone
+(
+ const IOobject& io,
+ const point& point1,
+ const scalar radius1,
+ const scalar innerRadius1,
+ const point& point2,
+ const scalar radius2,
+ const scalar innerRadius2
+)
+:
+ searchableSurface(io),
+ point1_(point1),
+ radius1_(radius1),
+ innerRadius1_(innerRadius1),
+ point2_(point2),
+ radius2_(radius2),
+ innerRadius2_(innerRadius2),
+ magDir_(mag(point2_-point1_)),
+ unitDir_((point2_-point1_)/magDir_)//wrap up the lookup function with read scalar
+{
+ bounds() = calcBounds();
+}
+
+
+Foam::searchableCone::searchableCone
+(
+ const IOobject& io,
+ const dictionary& dict
+)
+:
+ searchableSurface(io),
+ point1_(dict.lookup("point1")),
+ radius1_(readScalar(dict.lookup("radius1"))),
+ innerRadius1_(dict.lookupOrDefault( "innerRadius1", 0.0)),
+ point2_(dict.lookup("point2")),
+ radius2_(readScalar(dict.lookup("radius2"))),
+ innerRadius2_(dict.lookupOrDefault( "innerRadius2", 0.0)),
+ magDir_(mag(point2_-point1_)),
+ unitDir_((point2_-point1_)/magDir_) //wrap up the lookup function with read scalar
+
+{
+ bounds() = calcBounds();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
+
+const Foam::wordList& Foam::searchableCone::regions() const
+{
+ if (regions_.empty())
+ {
+ regions_.setSize(1);
+ regions_[0] = "region0";
+ }
+ return regions_;
+}
+
+
+void Foam::searchableCone::findNearest
+(
+ const pointField& samples,
+ const scalarField& nearestDistSqr,
+ List& info
+) const
+{
+ info.setSize(samples.size());
+
+ forAll(samples, i)
+ {
+ vector normal;
+ findNearestAndNormal(samples[i], nearestDistSqr[i], info[i], normal);
+ }
+}
+
+
+void Foam::searchableCone::findLine
+(
+ const pointField& start,
+ const pointField& end,
+ List& info
+) const
+{
+ info.setSize(start.size());
+
+ forAll(start, i)
+ {
+ // Pick nearest intersection. If none intersected take second one.
+ pointIndexHit b;
+ findLineAll
+ (
+ *this,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ info[i],
+ b
+ );
+ if (!info[i].hit() && b.hit())
+ {
+ info[i] = b;
+ }
+ }
+
+
+ if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
+ {
+ IOobject io(*this);
+ io.rename(name()+"Inner");
+ searchableCone innerCone
+ (
+ io,
+ point1_,
+ innerRadius1_,
+ 0.0,
+ point2_,
+ innerRadius2_,
+ 0.0
+ );
+
+ forAll(info, i)
+ {
+ point newEnd;
+ if (info[i].hit())
+ {
+ newEnd = info[i].hitPoint();
+ }
+ else
+ {
+ newEnd = end[i];
+ }
+ pointIndexHit near;
+ pointIndexHit far;
+ findLineAll
+ (
+ innerCone,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ newEnd,
+ near,
+ far
+ );
+
+ if (near.hit())
+ {
+ info[i] = near;
+ }
+ else if (far.hit())
+ {
+ info[i] = far;
+ }
+ }
+ }
+}
+
+
+void Foam::searchableCone::findLineAny
+(
+ const pointField& start,
+ const pointField& end,
+ List& info
+) const
+{
+ info.setSize(start.size());
+
+ forAll(start, i)
+ {
+ // Discard far intersection
+ pointIndexHit b;
+ findLineAll
+ (
+ *this,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ info[i],
+ b
+ );
+ if (!info[i].hit() && b.hit())
+ {
+ info[i] = b;
+ }
+ }
+ if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
+ {
+ IOobject io(*this);
+ io.rename(name()+"Inner");
+ searchableCone cone
+ (
+ io,
+ point1_,
+ innerRadius1_,
+ 0.0,
+ point2_,
+ innerRadius2_,
+ 0.0
+ );
+
+ forAll(info, i)
+ {
+ if (!info[i].hit())
+ {
+ pointIndexHit b;
+ findLineAll
+ (
+ cone,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ info[i],
+ b
+ );
+ if (!info[i].hit() && b.hit())
+ {
+ info[i] = b;
+ }
+ }
+ }
+ }
+}
+
+
+void Foam::searchableCone::findLineAll
+(
+ const pointField& start,
+ const pointField& end,
+ List>& info
+) const
+{
+ info.setSize(start.size());
+
+ forAll(start, i)
+ {
+ pointIndexHit near, far;
+ findLineAll
+ (
+ *this,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ near,
+ far
+ );
+
+ if (near.hit())
+ {
+ if (far.hit())
+ {
+ info[i].setSize(2);
+ info[i][0] = near;
+ info[i][1] = far;
+ }
+ else
+ {
+ info[i].setSize(1);
+ info[i][0] = near;
+ }
+ }
+ else
+ {
+ if (far.hit())
+ {
+ info[i].setSize(1);
+ info[i][0] = far;
+ }
+ else
+ {
+ info[i].clear();
+ }
+ }
+ }
+
+ if (innerRadius1_ > 0.0 || innerRadius2_ > 0.0)
+ {
+ IOobject io(*this);
+ io.rename(name()+"Inner");
+ searchableCone cone
+ (
+ io,
+ point1_,
+ innerRadius1_,
+ 0.0,
+ point2_,
+ innerRadius2_,
+ 0.0
+ );
+
+ forAll(info, i)
+ {
+ pointIndexHit near;
+ pointIndexHit far;
+ findLineAll
+ (
+ cone,
+ innerRadius1_,
+ innerRadius2_,
+ start[i],
+ end[i],
+ near,
+ far
+ );
+
+ if (near.hit())
+ {
+ insertHit(start[i], end[i], info[i], near);
+ }
+ if (far.hit())
+ {
+ insertHit(start[i], end[i], info[i], far);
+ }
+ }
+ }
+}
+
+
+void Foam::searchableCone::getRegion
+(
+ const List& info,
+ labelList& region
+) const
+{
+ region.setSize(info.size());
+ region = 0;
+}
+
+
+void Foam::searchableCone::getNormal
+(
+ const List& info,
+ vectorField& normal
+) const
+{
+ normal.setSize(info.size());
+ normal = Zero;
+
+ forAll(info, i)
+ {
+ if (info[i].hit())
+ {
+ pointIndexHit nearInfo;
+ findNearestAndNormal
+ (
+ info[i].hitPoint(),
+ Foam::sqr(GREAT),
+ nearInfo,
+ normal[i]
+ );
+ }
+ }
+}
+
+
+void Foam::searchableCone::getVolumeType
+(
+ const pointField& points,
+ List& volType
+) const
+{
+ volType.setSize(points.size());
+
+ forAll(points, pointi)
+ {
+ const point& pt = points[pointi];
+
+ volType[pointi] = volumeType::outside;
+
+ vector v(pt - point1_);
+
+ // Decompose sample-point1 into normal and parallel component
+ const scalar parallel = (v & unitDir_);
+
+ // Quick rejection. Left of point1 endcap, or right of point2 endcap
+ if (parallel < 0 || parallel > magDir_)
+ {
+ continue;
+ }
+
+ const scalar radius_sec =
+ radius1_ + parallel * (radius2_-radius1_)/magDir_;
+
+ const scalar radius_sec_inner =
+ innerRadius1_ + parallel * (innerRadius2_-innerRadius1_)/magDir_;
+
+ // Remove the parallel component
+ v -= parallel*unitDir_;
+
+ if (mag(v) >= radius_sec_inner && mag(v) <= radius_sec)
+ {
+ volType[pointi] = volumeType::inside;
+ }
+ }
+}
+
+
+// ************************************************************************* //
diff --git a/src/dfMeshTools/searchableSurfaces/searchableCone/searchableCone.H b/src/dfMeshTools/searchableSurfaces/searchableCone/searchableCone.H
new file mode 100644
index 000000000..c64a15d63
--- /dev/null
+++ b/src/dfMeshTools/searchableSurfaces/searchableCone/searchableCone.H
@@ -0,0 +1,298 @@
+/*---------------------------------------------------------------------------*\
+ ========= |
+ \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
+ \\ / O peration |
+ \\ / A nd | www.openfoam.com
+ \\/ M anipulation |
+-------------------------------------------------------------------------------
+ Copyright (C) 2015-2018 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+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::searchableCone
+
+Description
+ Searching on (optionally hollow) cone.
+
+ \heading Dictionary parameters
+ \table
+ Property | Description | Required | Default
+ type | cone | selector |
+ point1 | coordinate of endpoint | yes |
+ radius1 | radius at point1 | yes |
+ innerRadius1| inner radius at point1 | no | 0
+ point2 | coordinate of endpoint | yes |
+ radius2 | radius at point2 | yes |
+ innerRadius2| inner radius at point2 | no | 0
+ \endtable
+
+SourceFiles
+ searchableCone.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef searchableCone_H
+#define searchableCone_H
+
+#include "searchableSurface.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+ Class searchableCone Declaration
+\*---------------------------------------------------------------------------*/
+
+class searchableCone
+:
+ public searchableSurface
+{
+ // Private Member Data
+
+ //- 'Left' point
+ const point point1_;
+
+ //- Outer radius at point1
+ const scalar radius1_;
+
+ //- Inner radius at point1
+ const scalar innerRadius1_;
+
+
+ //- 'Right' point
+ const point point2_;
+
+ //- Outer radius at point2
+ const scalar radius2_;
+
+ //- Inner radius at point2
+ const scalar innerRadius2_;
+
+
+ //- Length of vector point2-point1
+ const scalar magDir_;
+
+ //- Normalised vector point2-point1
+ const vector unitDir_;
+
+
+ //- Names of regions
+ mutable wordList regions_;
+
+
+ // Private Member Functions
+
+ //- Find nearest point on cylinder.
+ void findNearestAndNormal
+ (
+ const point& sample,
+ const scalar nearestDistSqr,
+ pointIndexHit & nearInfo,
+ vector& normal
+ ) const;
+
+ //- Determine radial coordinate (squared)
+ static scalar radius2(const searchableCone& cone, const point& pt);
+
+ //- Find both intersections with cone. innerRadii supplied externally.
+ void findLineAll
+ (
+ const searchableCone& cone,
+ const scalar innerRadius1,
+ const scalar innerRadius2,
+ const point& start,
+ const point& end,
+ pointIndexHit& near,
+ pointIndexHit& far
+ ) const;
+
+ //- Insert a hit if it differs (by a tolerance) from the existing ones
+ void insertHit
+ (
+ const point& start,
+ const point& end,
+ List& info,
+ const pointIndexHit& hit
+ ) const;
+
+ //- Return the boundBox of the cylinder
+ boundBox calcBounds() const;
+
+ //- No copy construct
+ searchableCone(const searchableCone&) = delete;
+
+ //- No copy assignment
+ void operator=(const searchableCone&) = delete;
+
+
+public:
+
+ //- Runtime type information
+ TypeName("searchableCone");
+
+
+ // Constructors
+
+ //- Construct from components
+ searchableCone
+ (
+ const IOobject& io,
+ const point& point1,
+ const scalar radius1,
+ const scalar innerRadius1,
+ const point& point2,
+ const scalar radius2,
+ const scalar innerRadius2
+ );
+
+ //- Construct from dictionary (used by searchableSurface)
+ searchableCone
+ (
+ const IOobject& io,
+ const dictionary& dict
+ );
+
+
+ //- Destructor
+ virtual ~searchableCone() = default;
+
+
+ // Member Functions
+
+ //- Names of regions
+ virtual const wordList& regions() const;
+
+ //- Whether supports volume type (below)
+ virtual bool hasVolumeType() const
+ {
+ return true;
+ }
+
+ //- What is type of points outside bounds
+ virtual volumeType outsideVolumeType() const
+ {
+ return volumeType::outside;
+ }
+
+ //- Range of local indices that can be returned.
+ virtual label size() const
+ {
+ return 1;
+ }
+
+ //- Get representative set of element coordinates
+ // Usually the element centres (should be of length size()).
+ virtual tmp coordinates() const;
+
+ //- Get bounding spheres (centre and radius squared), one per element.
+ // Any point on element is guaranteed to be inside.
+ virtual void boundingSpheres
+ (
+ pointField& centres,
+ scalarField& radiusSqr
+ ) const;
+
+ //- Get the points that define the surface.
+ virtual tmp points() const;
+
+ //- Does any part of the surface overlap the supplied bound box?
+ virtual bool overlaps(const boundBox& bb) const
+ {
+ NotImplemented;
+ return false;
+ }
+
+
+ // Multiple point queries.
+
+ //- Find nearest point on cylinder
+ virtual void findNearest
+ (
+ const pointField& sample,
+ const scalarField& nearestDistSqr,
+ List&
+ ) const;
+
+ //- Find nearest intersection on line from start to end
+ virtual void findLine
+ (
+ const pointField& start,
+ const pointField& end,
+ List&
+ ) const;
+
+ //- Find all intersections in order from start to end
+ virtual void findLineAll
+ (
+ const pointField& start,
+ const pointField& end,
+ List>&
+ ) const;
+
+ //- Find any intersection on line from start to end
+ virtual void findLineAny
+ (
+ const pointField& start,
+ const pointField& end,
+ List&
+ ) const;
+
+ //- From a set of points and indices get the region
+ virtual void getRegion
+ (
+ const List&,
+ labelList& region
+ ) const;
+
+ //- From a set of points and indices get the normal
+ virtual void getNormal
+ (
+ const List&,
+ vectorField& normal
+ ) const;
+
+ //- Determine type (inside/outside/mixed) for point.
+ // Unknown if cannot be determined (e.g. non-manifold surface)
+ virtual void getVolumeType
+ (
+ const pointField& points,
+ List& volType
+ ) const;
+
+
+ // regIOobject implementation
+
+ virtual bool writeData(Ostream&) const
+ {
+ NotImplemented;
+ return false;
+ }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //