|
| 1 | +/*---------------------------------------------------------------------------*\ |
| 2 | + ========= | |
| 3 | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox |
| 4 | + \\ / O peration | Website: https://openfoam.org |
| 5 | + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation |
| 6 | + \\/ M anipulation | |
| 7 | +------------------------------------------------------------------------------- |
| 8 | +License |
| 9 | + This file is part of OpenFOAM. |
| 10 | +
|
| 11 | + OpenFOAM is free software: you can redistribute it and/or modify it |
| 12 | + under the terms of the GNU General Public License as published by |
| 13 | + the Free Software Foundation, either version 3 of the License, or |
| 14 | + (at your option) any later version. |
| 15 | +
|
| 16 | + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT |
| 17 | + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
| 18 | + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| 19 | + for more details. |
| 20 | +
|
| 21 | + You should have received a copy of the GNU General Public License |
| 22 | + along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>. |
| 23 | +
|
| 24 | +\*---------------------------------------------------------------------------*/ |
| 25 | + |
| 26 | +#include "cylinderCutToCell.H" |
| 27 | +#include "polyMesh.H" |
| 28 | +#include "addToRunTimeSelectionTable.H" |
| 29 | +#include "cellModel.H" |
| 30 | +#include "boundBox.H" |
| 31 | + |
| 32 | +// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * // |
| 33 | + |
| 34 | +namespace Foam |
| 35 | +{ |
| 36 | + defineTypeNameAndDebug(cylinderCutToCell, 0); |
| 37 | + addToRunTimeSelectionTable(topoSetSource, cylinderCutToCell, word); |
| 38 | + addToRunTimeSelectionTable(topoSetSource, cylinderCutToCell, istream); |
| 39 | +} |
| 40 | + |
| 41 | + |
| 42 | +Foam::topoSetSource::addToUsageTable Foam::cylinderCutToCell::usage_ |
| 43 | +( |
| 44 | + cylinderCutToCell::typeName, |
| 45 | + "\n Usage: cylinderCutToCell (p1X p1Y p1Z) (p2X p2Y p2Z) radius\n\n" |
| 46 | + " Select all cells that intersect with the bounding cylinder\n\n" |
| 47 | +); |
| 48 | + |
| 49 | + |
| 50 | +// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * // |
| 51 | + |
| 52 | +// Helper function: check if a point is inside the cylinder |
| 53 | +bool Foam::cylinderCutToCell::pointInCylinder |
| 54 | +( |
| 55 | + const point& pt, |
| 56 | + const vector& axis, |
| 57 | + const scalar magAxis2, |
| 58 | + const scalar rad2 |
| 59 | +) const |
| 60 | +{ |
| 61 | + vector d = pt - p1_; |
| 62 | + scalar magD = d & axis; |
| 63 | + |
| 64 | + if ((magD >= 0) && (magD <= magAxis2)) |
| 65 | + { |
| 66 | + scalar d2 = (d & d) - sqr(magD)/magAxis2; |
| 67 | + if (d2 <= rad2) |
| 68 | + { |
| 69 | + return true; |
| 70 | + } |
| 71 | + } |
| 72 | + return false; |
| 73 | +} |
| 74 | + |
| 75 | + |
| 76 | +// Helper function: check if cylinder intersects with a bounding box |
| 77 | +// Returns true only if the cylinder (not a capsule) intersects the box |
| 78 | +bool Foam::cylinderCutToCell::cylinderIntersectsBoundBox |
| 79 | +( |
| 80 | + const boundBox& bb, |
| 81 | + const vector& axis, |
| 82 | + const vector& axisDir, |
| 83 | + const scalar magAxis, |
| 84 | + const scalar rad2 |
| 85 | +) const |
| 86 | +{ |
| 87 | + // Check if the cylinder axis passes through (or close to) the bounding box |
| 88 | + // We need to ensure the intersection is within the cylinder length (not beyond ends) |
| 89 | + |
| 90 | + const point& boxMin = bb.min(); |
| 91 | + const point& boxMax = bb.max(); |
| 92 | + const point boxCenter = bb.midpoint(); |
| 93 | + const scalar boxRadius = mag(boxMax - boxMin) * 0.5; |
| 94 | + |
| 95 | + // Find projection of box center onto cylinder axis |
| 96 | + vector d = boxCenter - p1_; |
| 97 | + scalar t = (d & axisDir); // projection along normalized axis |
| 98 | + |
| 99 | + // Check if the box is completely outside the cylinder's axial range |
| 100 | + // Consider the box could extend into the cylinder range by boxRadius |
| 101 | + if (t < -boxRadius || t > magAxis + boxRadius) |
| 102 | + { |
| 103 | + return false; |
| 104 | + } |
| 105 | + |
| 106 | + // Clamp t to [0, magAxis] to find closest point on axis within cylinder |
| 107 | + scalar tClamped = max(scalar(0), min(t, magAxis)); |
| 108 | + |
| 109 | + // Closest point on axis (within cylinder length) |
| 110 | + point closestOnAxis = p1_ + tClamped * axisDir; |
| 111 | + |
| 112 | + // Distance from closest point to box center |
| 113 | + scalar dist2 = magSqr(closestOnAxis - boxCenter); |
| 114 | + |
| 115 | + // If distance is less than radius + box diagonal/2, they might intersect |
| 116 | + if (dist2 <= sqr(radius_ + boxRadius)) |
| 117 | + { |
| 118 | + return true; |
| 119 | + } |
| 120 | + |
| 121 | + return false; |
| 122 | +} |
| 123 | + |
| 124 | + |
| 125 | +void Foam::cylinderCutToCell::combine(topoSet& set, const bool add) const |
| 126 | +{ |
| 127 | + const vector axis = p2_ - p1_; |
| 128 | + const scalar rad2 = sqr(radius_); |
| 129 | + const scalar magAxis2 = magSqr(axis); |
| 130 | + const scalar magAxis = Foam::sqrt(magAxis2); |
| 131 | + const vector axisDir = axis / (magAxis + VSMALL); |
| 132 | + |
| 133 | + const pointField& pts = mesh_.points(); |
| 134 | + const cellList& cells = mesh_.cells(); |
| 135 | + const faceList& faces = mesh_.faces(); |
| 136 | + |
| 137 | + forAll(cells, celli) |
| 138 | + { |
| 139 | + bool selected = false; |
| 140 | + |
| 141 | + // First, quick check: does any vertex of the cell lie inside the cylinder? |
| 142 | + const cell& c = cells[celli]; |
| 143 | + labelHashSet cellPoints(c.size() * 4); |
| 144 | + |
| 145 | + // Collect all unique points of this cell |
| 146 | + forAll(c, facei) |
| 147 | + { |
| 148 | + const face& f = faces[c[facei]]; |
| 149 | + forAll(f, pointi) |
| 150 | + { |
| 151 | + cellPoints.insert(f[pointi]); |
| 152 | + } |
| 153 | + } |
| 154 | + |
| 155 | + // Check if any vertex is inside the cylinder |
| 156 | + forAllConstIter(labelHashSet, cellPoints, iter) |
| 157 | + { |
| 158 | + if (pointInCylinder(pts[iter.key()], axis, magAxis2, rad2)) |
| 159 | + { |
| 160 | + selected = true; |
| 161 | + break; |
| 162 | + } |
| 163 | + } |
| 164 | + |
| 165 | + // If no vertex inside, check if cylinder axis intersects cell bounding box |
| 166 | + if (!selected) |
| 167 | + { |
| 168 | + // Build bounding box for this cell - collect points first |
| 169 | + pointField cellPts(cellPoints.size()); |
| 170 | + label pi = 0; |
| 171 | + forAllConstIter(labelHashSet, cellPoints, iter) |
| 172 | + { |
| 173 | + cellPts[pi++] = pts[iter.key()]; |
| 174 | + } |
| 175 | + boundBox cellBB(cellPts, false); |
| 176 | + |
| 177 | + if (cylinderIntersectsBoundBox(cellBB, axis, axisDir, magAxis, rad2)) |
| 178 | + { |
| 179 | + // More detailed check: find the closest point on cylinder axis |
| 180 | + // to the cell center and verify it's within cylinder bounds |
| 181 | + const point& cellCtr = mesh_.cellCentres()[celli]; |
| 182 | + vector d = cellCtr - p1_; |
| 183 | + scalar t = (d & axisDir); |
| 184 | + |
| 185 | + // Only proceed if cell center projection is within cylinder length |
| 186 | + // (with some tolerance based on cell size) |
| 187 | + scalar cellRadius = mag(cellBB.max() - cellBB.min()) * 0.5; |
| 188 | + if (t >= -cellRadius && t <= magAxis + cellRadius) |
| 189 | + { |
| 190 | + // Clamp to cylinder range |
| 191 | + scalar tClamped = max(scalar(0), min(t, magAxis)); |
| 192 | + point closestOnAxis = p1_ + tClamped * axisDir; |
| 193 | + |
| 194 | + // Check radial distance from axis to cell center |
| 195 | + scalar radialDist2 = magSqr(cellCtr - closestOnAxis); |
| 196 | + |
| 197 | + // Select if within radius + cell size tolerance |
| 198 | + if (radialDist2 <= sqr(radius_ + cellRadius)) |
| 199 | + { |
| 200 | + selected = true; |
| 201 | + } |
| 202 | + } |
| 203 | + } |
| 204 | + } |
| 205 | + |
| 206 | + if (selected) |
| 207 | + { |
| 208 | + addOrDelete(set, celli, add); |
| 209 | + } |
| 210 | + } |
| 211 | +} |
| 212 | + |
| 213 | + |
| 214 | +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // |
| 215 | + |
| 216 | +Foam::cylinderCutToCell::cylinderCutToCell |
| 217 | +( |
| 218 | + const polyMesh& mesh, |
| 219 | + const vector& p1, |
| 220 | + const vector& p2, |
| 221 | + const scalar radius |
| 222 | +) |
| 223 | +: |
| 224 | + topoSetSource(mesh), |
| 225 | + p1_(p1), |
| 226 | + p2_(p2), |
| 227 | + radius_(radius) |
| 228 | +{} |
| 229 | + |
| 230 | + |
| 231 | +Foam::cylinderCutToCell::cylinderCutToCell |
| 232 | +( |
| 233 | + const polyMesh& mesh, |
| 234 | + const dictionary& dict |
| 235 | +) |
| 236 | +: |
| 237 | + topoSetSource(mesh), |
| 238 | + p1_(dict.lookup("p1")), |
| 239 | + p2_(dict.lookup("p2")), |
| 240 | + radius_(readScalar(dict.lookup("radius"))) |
| 241 | +{} |
| 242 | + |
| 243 | + |
| 244 | +Foam::cylinderCutToCell::cylinderCutToCell |
| 245 | +( |
| 246 | + const polyMesh& mesh, |
| 247 | + Istream& is |
| 248 | +) |
| 249 | +: |
| 250 | + topoSetSource(mesh), |
| 251 | + p1_(checkIs(is)), |
| 252 | + p2_(checkIs(is)), |
| 253 | + radius_(readScalar(checkIs(is))) |
| 254 | +{} |
| 255 | + |
| 256 | + |
| 257 | +// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * // |
| 258 | + |
| 259 | +Foam::cylinderCutToCell::~cylinderCutToCell() |
| 260 | +{} |
| 261 | + |
| 262 | + |
| 263 | +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // |
| 264 | + |
| 265 | +void Foam::cylinderCutToCell::applyToSet |
| 266 | +( |
| 267 | + const topoSetSource::setAction action, |
| 268 | + topoSet& set |
| 269 | +) const |
| 270 | +{ |
| 271 | + if ((action == topoSetSource::NEW) || (action == topoSetSource::ADD)) |
| 272 | + { |
| 273 | + Info<< " Adding cells intersecting with cylinder, with p1 = " |
| 274 | + << p1_ << ", p2 = " << p2_ << " and radius = " << radius_ << endl; |
| 275 | + |
| 276 | + combine(set, true); |
| 277 | + } |
| 278 | + else if (action == topoSetSource::DELETE) |
| 279 | + { |
| 280 | + Info<< " Removing cells intersecting with cylinder, with p1 = " |
| 281 | + << p1_ << ", p2 = " << p2_ << " and radius = " << radius_ << endl; |
| 282 | + |
| 283 | + combine(set, false); |
| 284 | + } |
| 285 | +} |
| 286 | + |
| 287 | + |
| 288 | +// ************************************************************************* // |
0 commit comments