sampledSet.C 13.6 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
Mark Olesen's avatar
Mark Olesen committed
9
    Copyright (C) 2018-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 32 33

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

#include "sampledSet.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "meshSearch.H"
#include "writer.H"
34
#include "particle.H"
35
#include "globalIndex.H"
36 37 38 39 40 41 42 43 44 45

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

namespace Foam
{
    defineTypeNameAndDebug(sampledSet, 0);
    defineRunTimeSelectionTable(sampledSet, word);
}


Henry Weller's avatar
Henry Weller committed
46
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
47

48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69
void Foam::sampledSet::checkDimensions() const
{
    if
    (
        (cells_.size() != size())
     || (faces_.size() != size())
     || (segments_.size() != size())
     || (curveDist_.size() != size())
    )
    {
        FatalErrorInFunction
            << "sizes not equal : "
            << "  points:" << size()
            << "  cells:" << cells_.size()
            << "  faces:" << faces_.size()
            << "  segments:" << segments_.size()
            << "  curveDist:" << curveDist_.size()
            << abort(FatalError);
    }
}


70
Foam::label Foam::sampledSet::getBoundaryCell(const label facei) const
71
{
72
    return mesh().faceOwner()[facei];
73 74 75
}


76
Foam::label Foam::sampledSet::getNeighbourCell(const label facei) const
77
{
78
    if (facei >= mesh().nInternalFaces())
79
    {
80
        return mesh().faceOwner()[facei];
81 82 83
    }
    else
    {
84
        return mesh().faceNeighbour()[facei];
85 86 87 88
    }
}


89
Foam::label Foam::sampledSet::pointInCell
90
(
91 92
    const point& p,
    const label samplei
93 94
) const
{
95 96
    // Collect the face owner and neighbour cells of the sample into an array
    // for convenience
97
    const label cells[4] =
98 99 100 101 102 103 104 105 106
    {
        mesh().faceOwner()[faces_[samplei]],
        getNeighbourCell(faces_[samplei]),
        mesh().faceOwner()[faces_[samplei+1]],
        getNeighbourCell(faces_[samplei+1])
    };

    // Find the sampled cell by checking the owners and neighbours of the
    // sampled faces
107
    label cellm =
108 109
        (cells[0] == cells[2] || cells[0] == cells[3]) ? cells[0]
      : (cells[1] == cells[2] || cells[1] == cells[3]) ? cells[1]
110
      : -1;
111

112 113
    if (cellm != -1)
    {
114 115
        // If found the sampled cell check the point is in the cell
        // otherwise ignore
116
        if (!mesh().pointInCell(p, cellm, searchEngine_.decompMode()))
117
        {
118
            cellm = -1;
119

120
            if (debug)
121
            {
122 123 124
                WarningInFunction
                    << "Could not find mid-point " << p
                    << " cell " << cellm << endl;
125 126 127
            }
        }
    }
128
    else
129
    {
130 131
        // If the sample does not pass through a single cell check if the point
        // is in any of the owners or neighbours otherwise ignore
132
        for (label i=0; i<4; ++i)
133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154
        {
            if (mesh().pointInCell(p, cells[i], searchEngine_.decompMode()))
            {
                return cells[i];
            }
        }

        if (debug)
        {
            WarningInFunction
                << "Could not find cell for mid-point" << nl
                << "  samplei: " << samplei
                << "  pts[samplei]: " << operator[](samplei)
                << "  face[samplei]: " << faces_[samplei]
                << "  pts[samplei+1]: " << operator[](samplei+1)
                << "  face[samplei+1]: " << faces_[samplei+1]
                << "  cellio: " << cells[0]
                << "  cellin: " << cells[1]
                << "  celljo: " << cells[2]
                << "  celljn: " << cells[3]
                << endl;
        }
155 156 157
    }

    return cellm;
158 159 160 161 162
}


Foam::scalar Foam::sampledSet::calcSign
(
163
    const label facei,
164 165 166
    const point& sample
) const
{
167
    vector vec = sample - mesh().faceCentres()[facei];
168 169 170 171 172

    scalar magVec = mag(vec);

    if (magVec < VSMALL)
    {
173
        // Sample on face centre. Regard as inside
174 175 176 177 178
        return -1;
    }

    vec /= magVec;

179
    const vector n = normalised(mesh().faceAreas()[facei]);
180 181 182 183 184 185 186

    return n & vec;
}


Foam::label Foam::sampledSet::findNearFace
(
187
    const label celli,
188 189 190 191
    const point& sample,
    const scalar smallDist
) const
{
192
    const cell& myFaces = mesh().cells()[celli];
193

194
    forAll(myFaces, myFacei)
195
    {
196
        const face& f = mesh().faces()[myFaces[myFacei]];
197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212

        pointHit inter = f.nearestPoint(sample, mesh().points());

        scalar dist;

        if (inter.hit())
        {
            dist = mag(inter.hitPoint() - sample);
        }
        else
        {
            dist = mag(inter.missPoint() - sample);
        }

        if (dist < smallDist)
        {
213
            return myFaces[myFacei];
214 215 216 217 218 219 220 221 222
        }
    }
    return -1;
}


Foam::point Foam::sampledSet::pushIn
(
    const point& facePt,
223
    const label facei
224 225
) const
{
226 227
    label celli = mesh().faceOwner()[facei];
    const point& cC = mesh().cellCentres()[celli];
228

229
    point newPosition = facePt;
230

231
    // Taken from particle::initCellFacePt()
232
    label tetFacei;
233
    label tetPtI;
234
    mesh().findTetFacePt(celli, facePt, tetFacei, tetPtI);
235

236 237 238 239 240 241 242 243
    // This is the tolerance that was defined as a static constant of the
    // particle class. It is no longer used by particle, following the switch to
    // barycentric tracking. The only place in which the tolerance is now used
    // is here. I'm not sure what the purpose of this code is, but it probably
    // wants removing. It is doing tet-searches for a particle position. This
    // should almost certainly be left to the particle class.
    const scalar trackingCorrectionTol = 1e-5;

244
    if (tetFacei == -1 || tetPtI == -1)
245 246 247
    {
        newPosition = facePt;

248
        label trap(1.0/trackingCorrectionTol + 1);
249 250 251 252 253

        label iterNo = 0;

        do
        {
254
            newPosition += trackingCorrectionTol*(cC - facePt);
255 256 257

            mesh().findTetFacePt
            (
258
                celli,
259
                newPosition,
260
                tetFacei,
261 262 263
                tetPtI
            );

264
            ++iterNo;
265

266
        } while (tetFacei < 0  && iterNo <= trap);
267 268
    }

269
    if (tetFacei == -1)
270
    {
271 272
        FatalErrorInFunction
            << "After pushing " << facePt << " to " << newPosition
273 274 275
            << " it is still outside face " << facei
            << " at " << mesh().faceCentres()[facei]
            << " of cell " << celli
276
            << " at " << cC << endl
277 278 279
            << "Please change your starting point"
            << abort(FatalError);
    }
280 281

    return newPosition;
282 283 284 285 286 287 288
}


bool Foam::sampledSet::getTrackingPoint
(
    const point& samplePt,
    const point& bPoint,
289
    const label bFacei,
290
    const scalar smallDist,
291 292

    point& trackPt,
293 294
    label& trackCelli,
    label& trackFacei
295 296 297 298
) const
{
    bool isGoodSample = false;

299
    if (bFacei == -1)
300 301
    {
        // No boundary intersection. Try and find cell samplePt is in
302
        trackCelli = mesh().findCell(samplePt, searchEngine_.decompMode());
303 304 305

        if
        (
306
            (trackCelli == -1)
307 308 309
        || !mesh().pointInCell
            (
                samplePt,
310
                trackCelli,
311 312
                searchEngine_.decompMode()
            )
313 314 315 316 317
        )
        {
            // Line samplePt - end_ does not intersect domain at all.
            // (or is along edge)

318 319
            trackCelli = -1;
            trackFacei = -1;
320 321 322 323 324

            isGoodSample = false;
        }
        else
        {
325
            // Start is inside. Use it as tracking point
326 327

            trackPt = samplePt;
328
            trackFacei = -1;
329 330 331 332 333 334 335

            isGoodSample = true;
        }
    }
    else if (mag(samplePt - bPoint) < smallDist)
    {
        // samplePt close to bPoint. Snap to it
336 337 338
        trackPt = pushIn(bPoint, bFacei);
        trackFacei = bFacei;
        trackCelli = getBoundaryCell(trackFacei);
339 340 341 342 343

        isGoodSample = true;
    }
    else
    {
344
        scalar sign = calcSign(bFacei, samplePt);
345 346 347 348 349

        if (sign < 0)
        {
            // samplePt inside or marginally outside.
            trackPt = samplePt;
350 351
            trackFacei = -1;
            trackCelli = mesh().findCell(trackPt, searchEngine_.decompMode());
352 353 354 355 356 357

            isGoodSample = true;
        }
        else
        {
            // samplePt outside. use bPoint
358 359 360
            trackPt = pushIn(bPoint, bFacei);
            trackFacei = bFacei;
            trackCelli = getBoundaryCell(trackFacei);
361 362 363 364 365

            isGoodSample = false;
        }
    }

Mark Olesen's avatar
Mark Olesen committed
366 367 368 369 370 371 372 373 374 375
    DebugInFunction
        << " samplePt:" << samplePt
        << " bPoint:" << bPoint
        << " bFacei:" << bFacei << nl
        << "   Calculated first tracking point :"
        << " trackPt:" << trackPt
        << " trackCelli:" << trackCelli
        << " trackFacei:" << trackFacei
        << " isGoodSample:" << isGoodSample
        << endl;
376 377 378 379 380 381 382 383 384 385 386 387 388 389

    return isGoodSample;
}


void Foam::sampledSet::setSamples
(
    const List<point>& samplingPts,
    const labelList& samplingCells,
    const labelList& samplingFaces,
    const labelList& samplingSegments,
    const scalarList& samplingCurveDist
)
{
390
    setPoints(samplingPts);
391 392
    curveDist_ = samplingCurveDist;

393
    segments_ = samplingSegments;
394 395
    cells_ = samplingCells;
    faces_ = samplingFaces;
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417

    checkDimensions();
}


void Foam::sampledSet::setSamples
(
    List<point>&& samplingPts,
    labelList&& samplingCells,
    labelList&& samplingFaces,
    labelList&& samplingSegments,
    scalarList&& samplingCurveDist
)
{
    setPoints(std::move(samplingPts));
    curveDist_ = std::move(samplingCurveDist);

    segments_ = std::move(samplingSegments);
    cells_ = std::move(samplingCells);
    faces_ = std::move(samplingFaces);

    checkDimensions();
418 419 420
}


421 422
Foam::autoPtr<Foam::coordSet> Foam::sampledSet::gather
(
423 424
    labelList& indexSet,
    labelList& allSegments
425 426 427 428 429 430
) const
{
    // Combine sampleSet from processors. Sort by curveDist. Return
    // ordering in indexSet.
    // Note: only master results are valid

431 432
    List<point> allPts;
    globalIndex::gatherOp(*this, allPts);
433

434
    globalIndex::gatherOp(segments(), allSegments);
435

436 437
    scalarList allCurveDist;
    globalIndex::gatherOp(curveDist(), allCurveDist);
438 439


440
    if (Pstream::master() && allCurveDist.empty())
441 442 443 444 445 446 447
    {
        WarningInFunction
            << "Sample set " << name()
            << " has zero points." << endl;
    }

    // Sort curveDist and use to fill masterSamplePts
448 449
    Foam::sortedOrder(allCurveDist, indexSet);      // uses stable sort
    scalarList sortedDist(allCurveDist, indexSet);  // with indices for mapping
450

451 452
    allSegments = UIndirectList<label>(allSegments, indexSet)();

453
    return autoPtr<coordSet>::New
454
    (
455 456 457 458
        name(),
        axis(),
        List<point>(UIndirectList<point>(allPts, indexSet)),
        sortedDist
459 460 461 462
    );
}


463 464 465 466 467 468
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::sampledSet::sampledSet
(
    const word& name,
    const polyMesh& mesh,
469
    const meshSearch& searchEngine,
470
    const coordSet::coordFormat axisType
471 472
)
:
473
    coordSet(name, axisType),
474 475
    mesh_(mesh),
    searchEngine_(searchEngine),
476 477 478
    segments_(),
    cells_(),
    faces_()
479 480 481 482 483 484 485
{}


Foam::sampledSet::sampledSet
(
    const word& name,
    const polyMesh& mesh,
486
    const meshSearch& searchEngine,
487
    const word& axis
488 489
)
:
490
    coordSet(name, axis),
491 492
    mesh_(mesh),
    searchEngine_(searchEngine),
493 494 495
    segments_(),
    cells_(),
    faces_()
496 497 498
{}


499 500 501 502 503 504 505 506 507 508 509 510 511 512
Foam::sampledSet::sampledSet
(
    const word& name,
    const polyMesh& mesh,
    const meshSearch& searchEngine,
    const dictionary& dict
)
:
    coordSet(name, dict.get<word>("axis")),
    mesh_(mesh),
    searchEngine_(searchEngine),
    segments_(),
    cells_(),
    faces_()
513 514 515 516 517
{}


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

518
Foam::autoPtr<Foam::sampledSet> Foam::sampledSet::New
519 520 521
(
    const word& name,
    const polyMesh& mesh,
522
    const meshSearch& searchEngine,
523 524 525
    const dictionary& dict
)
{
526
    const word sampleType(dict.get<word>("type"));
527

528
    auto cstrIter = wordConstructorTablePtr_->cfind(sampleType);
529

530
    if (!cstrIter.found())
531
    {
532
        FatalIOErrorInLookup
533
        (
534
            dict,
535 536 537
            "sample",
            sampleType,
            *wordConstructorTablePtr_
538
        ) << exit(FatalIOError);
539 540 541 542 543 544 545 546 547
    }

    return autoPtr<sampledSet>
    (
        cstrIter()
        (
            name,
            mesh,
            searchEngine,
548
            dict.optionalSubDict(sampleType + "Coeffs")
549 550 551 552 553 554 555 556 557
        )
    );
}


Foam::Ostream& Foam::sampledSet::write(Ostream& os) const
{
    coordSet::write(os);

558
    os  << nl << "\t(celli)\t(facei)" << nl;
559

560
    forAll(*this, samplei)
561
    {
562 563 564
        os  << '\t' << cells_[samplei]
            << '\t' << faces_[samplei]
            << nl;
565 566 567 568 569 570 571
    }

    return os;
}


// ************************************************************************* //