triSurfaceSearch.C 11.3 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
OpenFOAM bot's avatar
OpenFOAM bot committed
5
    \\  /    A nd           | www.openfoam.com
OpenFOAM bot's avatar
OpenFOAM bot committed
6 7
     \\/     M anipulation  |
-------------------------------------------------------------------------------
OpenFOAM bot's avatar
OpenFOAM bot committed
8
    Copyright (C) 2011-2017 OpenFOAM Foundation
9
    Copyright (C) 2015-2020 OpenCFD Ltd.
10 11 12 13
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

14 15 16 17
    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.
18 19 20 21 22 23 24

    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
25
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
26 27 28 29 30

\*---------------------------------------------------------------------------*/

#include "triSurfaceSearch.H"
#include "triSurface.H"
31 32
#include "PatchTools.H"
#include "volumeType.H"
33

34 35 36 37 38
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

bool Foam::triSurfaceSearch::checkUniqueHit
(
    const pointIndexHit& currHit,
39
    const UList<pointIndexHit>& hits,
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56
    const vector& lineVec
) const
{
    // Classify the hit
    label nearType = -1;
    label nearLabel = -1;

    const labelledTri& f = surface()[currHit.index()];

    f.nearestPointClassify
    (
        currHit.hitPoint(),
        surface().points(),
        nearType,
        nearLabel
    );

laurence's avatar
laurence committed
57
    if (nearType == triPointRef::POINT)
58 59
    {
        // near point
60

61
        const label nearPointi = f[nearLabel];
62 63

        const labelList& pointFaces =
64
            surface().pointFaces()[surface().meshPointMap()[nearPointi]];
65 66 67

        forAll(pointFaces, pI)
        {
68
            const label pointFacei = pointFaces[pI];
69

70
            if (pointFacei != currHit.index())
71 72 73 74 75
            {
                forAll(hits, hI)
                {
                    const pointIndexHit& hit = hits[hI];

76
                    if (hit.index() == pointFacei)
77 78 79 80 81 82
                    {
                        return false;
                    }
                }
            }
        }
83
    }
laurence's avatar
laurence committed
84
    else if (nearType == triPointRef::EDGE)
85 86 87 88 89 90 91 92 93 94 95 96
    {
        // near edge
        // check if the other face of the edge is already hit

        const labelList& fEdges = surface().faceEdges()[currHit.index()];

        const label edgeI = fEdges[nearLabel];

        const labelList& edgeFaces = surface().edgeFaces()[edgeI];

        forAll(edgeFaces, fI)
        {
97
            const label edgeFacei = edgeFaces[fI];
98

99
            if (edgeFacei != currHit.index())
100 101 102 103 104
            {
                forAll(hits, hI)
                {
                    const pointIndexHit& hit = hits[hI];

105
                    if (hit.index() == edgeFacei)
106 107 108 109 110 111
                    {
                        // Check normals
                        const vector currHitNormal =
                            surface().faceNormals()[currHit.index()];

                        const vector existingHitNormal =
112
                            surface().faceNormals()[edgeFacei];
113 114

                        const label signCurrHit =
Henry Weller's avatar
Henry Weller committed
115
                            pos0(currHitNormal & lineVec);
116 117

                        const label signExistingHit =
Henry Weller's avatar
Henry Weller committed
118
                            pos0(existingHitNormal & lineVec);
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133

                        if (signCurrHit == signExistingHit)
                        {
                            return false;
                        }
                    }
                }
            }
        }
    }

    return true;
}


134
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
135

136 137 138 139 140
Foam::triSurfaceSearch::triSurfaceSearch(const triSurface& surface)
:
    surface_(surface),
    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
    maxTreeDepth_(10),
141
    treePtr_(nullptr)
142
{}
143 144


145 146 147 148 149 150 151 152 153
Foam::triSurfaceSearch::triSurfaceSearch
(
    const triSurface& surface,
    const dictionary& dict
)
:
    surface_(surface),
    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
    maxTreeDepth_(10),
154
    treePtr_(nullptr)
155 156 157 158 159 160
{
    // Have optional non-standard search tolerance for gappy surfaces.
    if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
    {
        Info<< "    using intersection tolerance " << tolerance_ << endl;
    }
161

162 163 164 165 166 167 168 169 170 171 172 173 174 175
    // Have optional non-standard tree-depth to limit storage.
    if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
    {
        Info<< "    using maximum tree depth " << maxTreeDepth_ << endl;
    }
}


Foam::triSurfaceSearch::triSurfaceSearch
(
    const triSurface& surface,
    const scalar tolerance,
    const label maxTreeDepth
)
176 177
:
    surface_(surface),
178 179
    tolerance_(tolerance),
    maxTreeDepth_(maxTreeDepth),
180
    treePtr_(nullptr)
181 182 183 184 185 186
{
    if (tolerance_ < 0)
    {
        tolerance_ = indexedOctree<treeDataTriSurface>::perturbTol();
    }
}
187 188 189 190 191


// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //

Foam::triSurfaceSearch::~triSurfaceSearch()
192
{
193 194 195 196 197 198 199
    clearOut();
}


void Foam::triSurfaceSearch::clearOut()
{
    treePtr_.clear();
200 201 202 203 204
}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

205 206 207 208 209 210
const Foam::indexedOctree<Foam::treeDataTriSurface>&
Foam::triSurfaceSearch::tree() const
{
    if (treePtr_.empty())
    {
        // Calculate bb without constructing local point numbering.
Henry Weller's avatar
Henry Weller committed
211
        treeBoundBox bb(Zero, Zero);
212

213
        if (surface().size())
214
        {
215 216
            label nPoints;
            PatchTools::calcBounds(surface(), bb, nPoints);
217

218 219
            if (nPoints != surface().points().size())
            {
220
                WarningInFunction
221 222 223 224 225 226
                    << "Surface does not have compact point numbering."
                    << " Of " << surface().points().size()
                    << " only " << nPoints
                    << " are used. This might give problems in some routines."
                    << endl;
            }
227

228 229 230 231 232 233
            // Random number generator. Bit dodgy since not exactly random ;-)
            Random rndGen(65431);

            // Slightly extended bb. Slightly off-centred just so on symmetric
            // geometry there are less face/edge aligned items.
            bb = bb.extend(rndGen, 1e-4);
234 235
            bb.min() -= point::uniform(ROOTVSMALL);
            bb.max() += point::uniform(ROOTVSMALL);
236
        }
237

238
        const scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
239 240 241 242 243 244
        indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;

        treePtr_.reset
        (
            new indexedOctree<treeDataTriSurface>
            (
245
                treeDataTriSurface(false, surface_, tolerance_),
246 247 248 249 250 251 252 253 254 255
                bb,
                maxTreeDepth_,  // maxLevel
                10,             // leafsize
                3.0             // duplicity
            )
        );

        indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
    }

256
    return *treePtr_;
257 258 259
}


260
// Determine inside/outside for samples
261 262 263 264
Foam::boolList Foam::triSurfaceSearch::calcInside
(
    const pointField& samples
) const
265 266 267 268 269 270 271
{
    boolList inside(samples.size());

    forAll(samples, sampleI)
    {
        const point& sample = samples[sampleI];

mattijs's avatar
mattijs committed
272
        if (!tree().bb().contains(sample))
273 274 275
        {
            inside[sampleI] = false;
        }
276
        else if (tree().getVolumeType(sample) == volumeType::INSIDE)
277 278 279 280 281 282 283 284 285 286 287 288
        {
            inside[sampleI] = true;
        }
        else
        {
            inside[sampleI] = false;
        }
    }
    return inside;
}


289
void Foam::triSurfaceSearch::findNearest
290 291
(
    const pointField& samples,
292 293
    const scalarField& nearestDistSqr,
    List<pointIndexHit>& info
294 295
) const
{
296
    const scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
297
    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
298

299
    const indexedOctree<treeDataTriSurface>& octree = tree();
300

301 302
    const treeDataTriSurface::findNearestOp fOp(octree);

303
    info.setSize(samples.size());
304

305
    forAll(samples, i)
306
    {
307
        info[i] = octree.findNearest
308 309 310
        (
            samples[i],
            nearestDistSqr[i],
311
            fOp
312
        );
313 314
    }

315
    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
316 317 318
}


319
Foam::pointIndexHit Foam::triSurfaceSearch::nearest
320
(
321
    const point& pt,
322
    const vector& span
323 324
)
const
325
{
mattijs's avatar
mattijs committed
326
    const scalar nearestDistSqr = 0.25*magSqr(span);
327

328 329
    return tree().findNearest(pt, nearestDistSqr);
}
330 331


332 333 334 335 336 337 338 339
void Foam::triSurfaceSearch::findLine
(
    const pointField& start,
    const pointField& end,
    List<pointIndexHit>& info
) const
{
    const indexedOctree<treeDataTriSurface>& octree = tree();
340

341 342
    info.setSize(start.size());

343
    const scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
344 345 346 347
    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();

    forAll(start, i)
    {
348
        info[i] = octree.findLine(start[i], end[i]);
349 350
    }

351
    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
352 353 354
}


355
void Foam::triSurfaceSearch::findLineAny
356
(
357 358 359 360
    const pointField& start,
    const pointField& end,
    List<pointIndexHit>& info
) const
361
{
362
    const indexedOctree<treeDataTriSurface>& octree = tree();
363

364 365
    info.setSize(start.size());

366
    const scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
367 368 369 370
    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();

    forAll(start, i)
    {
371
        info[i] = octree.findLineAny(start[i], end[i]);
372 373 374
    }

    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
375 376 377
}


378 379
void Foam::triSurfaceSearch::findLineAll
(
380 381
    const pointField& start,
    const pointField& end,
382
    List<List<pointIndexHit>>& info
383
) const
384
{
385
    const indexedOctree<treeDataTriSurface>& octree = tree();
386

387
    info.setSize(start.size());
388

389
    const scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
390
    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
391

392
    // Work array
393
    DynamicList<pointIndexHit> hits;
394

395
    DynamicList<label> shapeMask;
396

397
    treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
398

399
    forAll(start, pointi)
400 401 402
    {
        hits.clear();
        shapeMask.clear();
403

404 405
        while (true)
        {
406
            // See if any intersection between pt and end
407 408
            pointIndexHit inter = octree.findLine
            (
409 410
                start[pointi],
                end[pointi],
411 412
                allIntersectOp
            );
413

414
            if (inter.hit())
415
            {
416
                const vector lineVec = normalised(end[pointi] - start[pointi]);
417 418 419 420 421 422 423 424 425 426 427 428 429

                if
                (
                    checkUniqueHit
                    (
                        inter,
                        hits,
                        lineVec
                    )
                )
                {
                    hits.append(inter);
                }
430

431
                shapeMask.append(inter.index());
432 433 434
            }
            else
            {
435
                break;
436 437
            }
        }
438

439
        info[pointi].transfer(hits);
440
    }
441 442

    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
443 444 445
}


446
// ************************************************************************* //