treeDataPrimitivePatch.C 15.9 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-2016 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 31

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

#include "treeDataPrimitivePatch.H"
#include "indexedOctree.H"
#include "triangleFuncs.H"
laurence's avatar
laurence committed
32 33
#include "triSurfaceTools.H"
#include "triFace.H"
34 35 36

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

37 38
template<class PatchType>
void Foam::treeDataPrimitivePatch<PatchType>::update()
39 40 41 42 43 44 45
{
    if (cacheBb_)
    {
        bbs_.setSize(patch_.size());

        forAll(patch_, i)
        {
46
            bbs_[i] = treeBoundBox(patch_.points(), patch_[i]);
47 48 49 50 51 52 53
        }
    }
}


// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

54 55
template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::treeDataPrimitivePatch
56 57
(
    const bool cacheBb,
laurence's avatar
laurence committed
58 59
    const PatchType& patch,
    const scalar planarTol
60 61 62
)
:
    patch_(patch),
laurence's avatar
laurence committed
63 64
    cacheBb_(cacheBb),
    planarTol_(planarTol)
65 66 67 68 69
{
    update();
}


laurence's avatar
laurence committed
70 71 72
template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::findNearestOp
(
73
    const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree
laurence's avatar
laurence committed
74 75 76 77 78 79 80 81 82
)
:
    tree_(tree)
{}


template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::findIntersectOp::findIntersectOp
(
83
    const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree
laurence's avatar
laurence committed
84 85 86 87 88 89 90 91 92
)
:
    tree_(tree)
{}


template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::findAllIntersectOp::findAllIntersectOp
(
93
    const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree,
laurence's avatar
laurence committed
94 95 96 97 98 99 100 101
    DynamicList<label>& shapeMask
)
:
    tree_(tree),
    shapeMask_(shapeMask)
{}


102 103 104 105
template<class PatchType>
Foam::treeDataPrimitivePatch<PatchType>::
findSelfIntersectOp::findSelfIntersectOp
(
106
    const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree,
107 108 109 110 111 112 113 114
    const label edgeID
)
:
    tree_(tree),
    edgeID_(edgeID)
{}


115 116
// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

117 118
template<class PatchType>
Foam::pointField Foam::treeDataPrimitivePatch<PatchType>::shapePoints() const
119 120 121 122 123 124 125 126 127 128 129 130
{
    pointField cc(patch_.size());

    forAll(patch_, i)
    {
        cc[i] = patch_[i].centre(patch_.points());
    }

    return cc;
}


131
template<class PatchType>
132
Foam::volumeType Foam::treeDataPrimitivePatch<PatchType>::getVolumeType
133
(
134
    const indexedOctree<treeDataPrimitivePatch<PatchType>>& oc,
135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153
    const point& sample
) const
{
    // Need to determine whether sample is 'inside' or 'outside'
    // Done by finding nearest face. This gives back a face which is
    // guaranteed to contain nearest point. This point can be
    // - in interior of face: compare to face normal
    // - on edge of face: compare to edge normal
    // - on point of face: compare to point normal
    // Unfortunately the octree does not give us back the intersection point
    // or where on the face it has hit so we have to recreate all that
    // information.


    // Find nearest face to sample
    pointIndexHit info = oc.findNearest(sample, sqr(GREAT));

    if (info.index() == -1)
    {
154 155
        FatalErrorInFunction
            << "Could not find " << sample << " in octree."
156 157 158 159
            << abort(FatalError);
    }

    // Get actual intersection point on face
160
    label facei = info.index();
161 162 163 164

    if (debug & 2)
    {
        Pout<< "getSampleType : sample:" << sample
165
            << " nearest face:" << facei;
166 167
    }

168 169
    const typename PatchType::face_type& localF = patch_.localFaces()[facei];
    const typename PatchType::face_type& f = patch_[facei];
170 171
    const pointField& points = patch_.points();
    const labelList& mp = patch_.meshPoints();
172 173 174 175 176

    // Retest to classify where on face info is. Note: could be improved. We
    // already have point.

    pointHit curHit = f.nearestPoint(sample, points);
177
    const vector area = f.areaNormal(points);
178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211
    const point& curPt = curHit.rawPoint();

    //
    // 1] Check whether sample is above face
    //

    if (curHit.hit())
    {
        // Nearest point inside face. Compare to face normal.

        if (debug & 2)
        {
            Pout<< " -> face hit:" << curPt
                << " comparing to face normal " << area << endl;
        }
        return indexedOctree<treeDataPrimitivePatch>::getSide
        (
            area,
            sample - curPt
        );
    }

    if (debug & 2)
    {
        Pout<< " -> face miss:" << curPt;
    }

    //
    // 2] Check whether intersection is on one of the face vertices or
    //    face centre
    //

    const scalar typDimSqr = mag(area) + VSMALL;

laurence's avatar
laurence committed
212

213 214
    forAll(f, fp)
    {
laurence's avatar
laurence committed
215
        if ((magSqr(points[f[fp]] - curPt)/typDimSqr) < planarTol_)
216 217 218 219 220 221 222 223
        {
            // Face intersection point equals face vertex fp

            // Calculate point normal (wrong: uses face normals instead of
            // triangle normals)

            return indexedOctree<treeDataPrimitivePatch>::getSide
            (
224
                patch_.pointNormals()[localF[fp]],
225 226 227 228 229 230 231
                sample - curPt
            );
        }
    }

    const point fc(f.centre(points));

laurence's avatar
laurence committed
232
    if ((magSqr(fc - curPt)/typDimSqr) < planarTol_)
233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255
    {
        // Face intersection point equals face centre. Normal at face centre
        // is already average of face normals

        if (debug & 2)
        {
            Pout<< " -> centre hit:" << fc
                << " distance:" << magSqr(fc - curPt)/typDimSqr << endl;
        }

        return indexedOctree<treeDataPrimitivePatch>::getSide
        (
            area,
            sample - curPt
        );
    }



    //
    // 3] Get the 'real' edge the face intersection is on
    //

256
    const labelList& fEdges = patch_.faceEdges()[facei];
257 258 259 260 261

    forAll(fEdges, fEdgeI)
    {
        label edgeI = fEdges[fEdgeI];
        const edge& e = patch_.edges()[edgeI];
262 263
        const linePointRef ln(points[mp[e.start()]], points[mp[e.end()]]);
        pointHit edgeHit = ln.nearestDist(sample);
264

laurence's avatar
laurence committed
265
        if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
266 267 268 269 270 271 272
        {
            // Face intersection point lies on edge e

            // Calculate edge normal (wrong: uses face normals instead of
            // triangle normals)
            const labelList& eFaces = patch_.edgeFaces()[edgeI];

Henry Weller's avatar
Henry Weller committed
273
            vector edgeNormal(Zero);
274 275 276

            forAll(eFaces, i)
            {
graham's avatar
graham committed
277
                edgeNormal += patch_.faceNormals()[eFaces[i]];
278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303
            }

            if (debug & 2)
            {
                Pout<< " -> real edge hit point:" << edgeHit.rawPoint()
                    << " comparing to edge normal:" << edgeNormal
                    << endl;
            }

            // Found face intersection point on this edge. Compare to edge
            // normal
            return indexedOctree<treeDataPrimitivePatch>::getSide
            (
                edgeNormal,
                sample - curPt
            );
        }
    }


    //
    // 4] Get the internal edge the face intersection is on
    //

    forAll(f, fp)
    {
304
        pointHit edgeHit = linePointRef(points[f[fp]], fc).nearestDist(sample);
305

laurence's avatar
laurence committed
306
        if ((magSqr(edgeHit.rawPoint() - curPt)/typDimSqr) < planarTol_)
307 308 309 310 311 312 313 314
        {
            // Face intersection point lies on edge between two face triangles

            // Calculate edge normal as average of the two triangle normals
            vector e = points[f[fp]] - fc;
            vector ePrev = points[f[f.rcIndex(fp)]] - fc;
            vector eNext = points[f[f.fcIndex(fp)]] - fc;

315 316
            vector nLeft = normalised(ePrev ^ e);
            vector nRight = normalised(e ^ eNext);
317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338

            if (debug & 2)
            {
                Pout<< " -> internal edge hit point:" << edgeHit.rawPoint()
                    << " comparing to edge normal "
                    << 0.5*(nLeft + nRight)
                    << endl;
            }

            // Found face intersection point on this edge. Compare to edge
            // normal
            return indexedOctree<treeDataPrimitivePatch>::getSide
            (
                0.5*(nLeft + nRight),
                sample - curPt
            );
        }
    }

    if (debug & 2)
    {
        Pout<< "Did not find sample " << sample
339
            << " anywhere related to nearest face " << facei << endl
340 341 342 343
            << "Face:";

        forAll(f, fp)
        {
344 345
            Pout<< "    vertex:" << f[fp]
                << "  coord:" << points[f[fp]]
346 347 348 349 350 351 352 353 354
                << endl;
        }
    }

    // Can't determine status of sample with respect to nearest face.
    // Either
    // - tolerances are wrong. (if e.g. face has zero area)
    // - or (more likely) surface is not closed.

355
    return volumeType::UNKNOWN;
356 357 358 359
}


// Check if any point on shape is inside cubeBb.
360 361
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
362 363 364 365 366 367
(
    const label index,
    const treeBoundBox& cubeBb
) const
{
    // 1. Quick rejection: bb does not intersect face bb at all
368 369 370 371 372 373
    if
    (
        cacheBb_
      ? !cubeBb.overlaps(bbs_[index])
      : !cubeBb.overlaps(treeBoundBox(patch_.points(), patch_[index]))
    )
374
    {
375
        return false;
376 377 378 379 380 381
    }


    // 2. Check if one or more face points inside

    const pointField& points = patch_.points();
382
    const typename PatchType::face_type& f = patch_[index];
383

384
    if (cubeBb.containsAny(points, f))
385
    {
386
        return true;
387 388 389 390 391 392
    }

    // 3. Difficult case: all points are outside but connecting edges might
    // go through cube. Use triangle-bounding box intersection.
    const point fc = f.centre(points);

laurence's avatar
laurence committed
393
    if (f.size() == 3)
394
    {
laurence's avatar
laurence committed
395
        return triangleFuncs::intersectBb
396
        (
laurence's avatar
laurence committed
397 398 399
            points[f[0]],
            points[f[1]],
            points[f[2]],
400 401
            cubeBb
        );
laurence's avatar
laurence committed
402 403 404 405
    }
    else
    {
        forAll(f, fp)
406
        {
laurence's avatar
laurence committed
407 408 409 410 411 412 413 414 415 416 417 418
            bool triIntersects = triangleFuncs::intersectBb
            (
                points[f[fp]],
                points[f[f.fcIndex(fp)]],
                fc,
                cubeBb
            );

            if (triIntersects)
            {
                return true;
            }
419 420
        }
    }
laurence's avatar
laurence committed
421

422 423 424 425
    return false;
}


426
// Check if any point on shape is inside sphere.
427 428
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::overlaps
429 430 431 432 433 434 435
(
    const label index,
    const point& centre,
    const scalar radiusSqr
) const
{
    // 1. Quick rejection: sphere does not intersect face bb at all
436 437 438 439 440 441
    if
    (
        cacheBb_
      ? !bbs_[index].overlaps(centre, radiusSqr)
      : !treeBoundBox(patch_.points(),patch_[index]).overlaps(centre, radiusSqr)
    )
442
    {
443
        return false;
444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461
    }

    const pointField& points = patch_.points();
    const face& f = patch_[index];

    pointHit nearHit = f.nearestPoint(centre, points);

    // If the distance to the nearest point on the face from the sphere centres
    // is within the radius, then the sphere touches the face.
    if (sqr(nearHit.distance()) < radiusSqr)
    {
        return true;
    }

    return false;
}


462
template<class PatchType>
laurence's avatar
laurence committed
463
void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
464
(
465
    const labelUList& indices,
466 467 468 469 470 471 472
    const point& sample,

    scalar& nearestDistSqr,
    label& minIndex,
    point& nearestPoint
) const
{
laurence's avatar
laurence committed
473 474 475 476
    const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
    const PatchType& patch = shape.patch();

    const pointField& points = patch.points();
477

478
    for (const label index : indices)
479
    {
480
        const typename PatchType::face_type& f = patch[index];
481

482 483
        const pointHit nearHit = f.nearestPoint(sample, points);
        const scalar distSqr = sqr(nearHit.distance());
484 485 486 487 488 489 490 491 492 493 494

        if (distSqr < nearestDistSqr)
        {
            nearestDistSqr = distSqr;
            minIndex = index;
            nearestPoint = nearHit.rawPoint();
        }
    }
}


495
template<class PatchType>
laurence's avatar
laurence committed
496 497 498 499 500 501 502 503 504 505 506
void Foam::treeDataPrimitivePatch<PatchType>::findNearestOp::operator()
(
    const labelUList& indices,
    const linePointRef& ln,

    treeBoundBox& tightest,
    label& minIndex,
    point& linePoint,
    point& nearestPoint
) const
{
507
    NotImplemented;
laurence's avatar
laurence committed
508 509 510 511 512
}


template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::findIntersectOp::operator()
513 514 515 516 517 518 519
(
    const label index,
    const point& start,
    const point& end,
    point& intersectionPoint
) const
{
laurence's avatar
laurence committed
520 521
    return findIntersection(tree_, index, start, end, intersectionPoint);
}
522 523


laurence's avatar
laurence committed
524 525 526 527 528 529 530 531 532
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::findAllIntersectOp::operator()
(
    const label index,
    const point& start,
    const point& end,
    point& intersectionPoint
) const
{
533
    if (shapeMask_.found(index))
534 535 536
    {
        return false;
    }
laurence's avatar
laurence committed
537 538

    return findIntersection(tree_, index, start, end, intersectionPoint);
539 540 541
}


542 543 544 545 546 547 548 549 550 551 552
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::findSelfIntersectOp::operator()
(
    const label index,
    const point& start,
    const point& end,
    point& intersectionPoint
) const
{
    if (edgeID_ == -1)
    {
553 554
        FatalErrorInFunction
            << "EdgeID not set. Please set edgeID to the index of"
555 556 557 558 559 560 561
            << " the edge you are testing"
            << exit(FatalError);
    }

    const treeDataPrimitivePatch<PatchType>& shape = tree_.shapes();
    const PatchType& patch = shape.patch();

562
    const typename PatchType::face_type& f = patch.localFaces()[index];
563 564
    const edge& e = patch.edges()[edgeID_];

565
    if (!f.found(e[0]) && !f.found(e[1]))
566 567 568
    {
        return findIntersection(tree_, index, start, end, intersectionPoint);
    }
569 570

    return false;
571 572 573
}


574 575 576
template<class PatchType>
bool Foam::treeDataPrimitivePatch<PatchType>::findIntersection
(
577
    const indexedOctree<treeDataPrimitivePatch<PatchType>>& tree,
578 579 580 581 582 583 584 585 586 587
    const label index,
    const point& start,
    const point& end,
    point& intersectionPoint
)
{
    const treeDataPrimitivePatch<PatchType>& shape = tree.shapes();
    const PatchType& patch = shape.patch();

    const pointField& points = patch.points();
588
    const typename PatchType::face_type& f = patch[index];
589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636

    // Do quick rejection test
    if (shape.cacheBb_)
    {
        const treeBoundBox& faceBb = shape.bbs_[index];

        if ((faceBb.posBits(start) & faceBb.posBits(end)) != 0)
        {
            // start and end in same block outside of faceBb.
            return false;
        }
    }

    const vector dir(end - start);
    pointHit inter;

    if (f.size() == 3)
    {
        inter = triPointRef
        (
            points[f[0]],
            points[f[1]],
            points[f[2]]
        ).intersection(start, dir, intersection::HALF_RAY, shape.planarTol_);
    }
    else
    {
        const pointField& faceCentres = patch.faceCentres();

        inter = f.intersection
        (
            start,
            dir,
            faceCentres[index],
            points,
            intersection::HALF_RAY,
            shape.planarTol_
        );
    }

    if (inter.hit() && inter.distance() <= 1)
    {
        // Note: no extra test on whether intersection is in front of us
        // since using half_ray
        intersectionPoint = inter.hitPoint();

        return true;
    }
637 638

    return false;
639 640 641
}


642
// ************************************************************************* //