sampledSet.C 11.7 KB
Newer Older
1 2 3 4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
6 7 8 9 10
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

11 12 13 14
    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.
15 16 17 18 19 20 21

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

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

#include "sampledSet.H"
#include "polyMesh.H"
#include "primitiveMesh.H"
#include "meshSearch.H"
#include "writer.H"
31
#include "particle.H"
32 33 34 35 36 37 38 39 40 41

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

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


Henry Weller's avatar
Henry Weller committed
42
// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
43

44
Foam::label Foam::sampledSet::getBoundaryCell(const label facei) const
45
{
46
    return mesh().faceOwner()[facei];
47 48 49
}


50
Foam::label Foam::sampledSet::getNeighbourCell(const label facei) const
51
{
52
    if (facei >= mesh().nInternalFaces())
53
    {
54
        return mesh().faceOwner()[facei];
55 56 57
    }
    else
    {
58
        return mesh().faceNeighbour()[facei];
59 60 61 62
    }
}


63
Foam::label Foam::sampledSet::pointInCell
64
(
65 66
    const point& p,
    const label samplei
67 68
) const
{
69 70 71 72 73 74 75 76 77 78 79 80
    // Collect the face owner and neighbour cells of the sample into an array
    // for convenience
    label cells[4] =
    {
        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
81
    label cellm =
82 83
        (cells[0] == cells[2] || cells[0] == cells[3]) ? cells[0]
      : (cells[1] == cells[2] || cells[1] == cells[3]) ? cells[1]
84
      : -1;
85

86 87
    if (cellm != -1)
    {
88 89
        // If found the sampled cell check the point is in the cell
        // otherwise ignore
90
        if (!mesh().pointInCell(p, cellm, searchEngine_.decompMode()))
91
        {
92
           cellm = -1;
93

94
            if (debug)
95
            {
96 97 98
                WarningInFunction
                    << "Could not find mid-point " << p
                    << " cell " << cellm << endl;
99 100 101
            }
        }
    }
102
    else
103
    {
104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128
        // If the sample does not pass through a single cell check if the point
        // is in any of the owners or neighbours otherwise ignore
        for (label i=0; i<4; i++)
        {
            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;
        }
129 130 131
    }

    return cellm;
132 133 134 135 136
}


Foam::scalar Foam::sampledSet::calcSign
(
137
    const label facei,
138 139 140
    const point& sample
) const
{
141
    vector vec = sample - mesh().faceCentres()[facei];
142 143 144 145 146

    scalar magVec = mag(vec);

    if (magVec < VSMALL)
    {
147
        // Sample on face centre. Regard as inside
148 149 150 151 152
        return -1;
    }

    vec /= magVec;

153
    vector n = mesh().faceAreas()[facei];
154 155 156 157 158 159 160 161 162

    n /= mag(n) + VSMALL;

    return n & vec;
}


Foam::label Foam::sampledSet::findNearFace
(
163
    const label celli,
164 165 166 167
    const point& sample,
    const scalar smallDist
) const
{
168
    const cell& myFaces = mesh().cells()[celli];
169

170
    forAll(myFaces, myFacei)
171
    {
172
        const face& f = mesh().faces()[myFaces[myFacei]];
173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188

        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)
        {
189
            return myFaces[myFacei];
190 191 192 193 194 195 196 197 198
        }
    }
    return -1;
}


Foam::point Foam::sampledSet::pushIn
(
    const point& facePt,
199
    const label facei
200 201
) const
{
202 203
    label celli = mesh().faceOwner()[facei];
    const point& cC = mesh().cellCentres()[celli];
204

205
    point newPosition = facePt;
206

207
    // Taken from particle::initCellFacePt()
208
    label tetFacei;
209
    label tetPtI;
210
    mesh().findTetFacePt(celli, facePt, tetFacei, tetPtI);
211

212
    if (tetFacei == -1 || tetPtI == -1)
213 214 215 216 217 218 219 220 221 222 223 224 225
    {
        newPosition = facePt;

        label trap(1.0/particle::trackingCorrectionTol + 1);

        label iterNo = 0;

        do
        {
            newPosition += particle::trackingCorrectionTol*(cC - facePt);

            mesh().findTetFacePt
            (
226
                celli,
227
                newPosition,
228
                tetFacei,
229 230 231 232 233
                tetPtI
            );

            iterNo++;

234
        } while (tetFacei < 0  && iterNo <= trap);
235 236
    }

237
    if (tetFacei == -1)
238
    {
239 240
        FatalErrorInFunction
            << "After pushing " << facePt << " to " << newPosition
241 242 243
            << " it is still outside face " << facei
            << " at " << mesh().faceCentres()[facei]
            << " of cell " << celli
244
            << " at " << cC << endl
245 246 247
            << "Please change your starting point"
            << abort(FatalError);
    }
248 249

    return newPosition;
250 251 252 253 254 255 256
}


bool Foam::sampledSet::getTrackingPoint
(
    const point& samplePt,
    const point& bPoint,
257
    const label bFacei,
258
    const scalar smallDist,
259 260

    point& trackPt,
261 262
    label& trackCelli,
    label& trackFacei
263 264 265 266
) const
{
    bool isGoodSample = false;

267
    if (bFacei == -1)
268 269
    {
        // No boundary intersection. Try and find cell samplePt is in
270
        trackCelli = mesh().findCell(samplePt, searchEngine_.decompMode());
271 272 273

        if
        (
274
            (trackCelli == -1)
275 276 277
        || !mesh().pointInCell
            (
                samplePt,
278
                trackCelli,
279 280
                searchEngine_.decompMode()
            )
281 282 283 284 285
        )
        {
            // Line samplePt - end_ does not intersect domain at all.
            // (or is along edge)

286 287
            trackCelli = -1;
            trackFacei = -1;
288 289 290 291 292

            isGoodSample = false;
        }
        else
        {
293
            // Start is inside. Use it as tracking point
294 295

            trackPt = samplePt;
296
            trackFacei = -1;
297 298 299 300 301 302 303

            isGoodSample = true;
        }
    }
    else if (mag(samplePt - bPoint) < smallDist)
    {
        // samplePt close to bPoint. Snap to it
304 305 306
        trackPt = pushIn(bPoint, bFacei);
        trackFacei = bFacei;
        trackCelli = getBoundaryCell(trackFacei);
307 308 309 310 311

        isGoodSample = true;
    }
    else
    {
312
        scalar sign = calcSign(bFacei, samplePt);
313 314 315 316 317

        if (sign < 0)
        {
            // samplePt inside or marginally outside.
            trackPt = samplePt;
318 319
            trackFacei = -1;
            trackCelli = mesh().findCell(trackPt, searchEngine_.decompMode());
320 321 322 323 324 325

            isGoodSample = true;
        }
        else
        {
            // samplePt outside. use bPoint
326 327 328
            trackPt = pushIn(bPoint, bFacei);
            trackFacei = bFacei;
            trackCelli = getBoundaryCell(trackFacei);
329 330 331 332 333 334 335

            isGoodSample = false;
        }
    }

    if (debug)
    {
336
        InfoInFunction
337 338
            << " samplePt:" << samplePt
            << " bPoint:" << bPoint
339
            << " bFacei:" << bFacei
340 341
            << endl << "   Calculated first tracking point :"
            << " trackPt:" << trackPt
342 343
            << " trackCelli:" << trackCelli
            << " trackFacei:" << trackFacei
344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374
            << " isGoodSample:" << isGoodSample
            << endl;
    }

    return isGoodSample;
}


void Foam::sampledSet::setSamples
(
    const List<point>& samplingPts,
    const labelList& samplingCells,
    const labelList& samplingFaces,
    const labelList& samplingSegments,
    const scalarList& samplingCurveDist
)
{
    setSize(samplingPts.size());
    cells_.setSize(samplingCells.size());
    faces_.setSize(samplingFaces.size());
    segments_.setSize(samplingSegments.size());
    curveDist_.setSize(samplingCurveDist.size());

    if
    (
        (cells_.size() != size())
     || (faces_.size() != size())
     || (segments_.size() != size())
     || (curveDist_.size() != size())
    )
    {
375
        FatalErrorInFunction
376 377 378 379 380 381 382 383 384 385 386 387 388
            << "sizes not equal : "
            << "  points:" << size()
            << "  cells:" << cells_.size()
            << "  faces:" << faces_.size()
            << "  segments:" << segments_.size()
            << "  curveDist:" << curveDist_.size()
            << abort(FatalError);
    }

    forAll(samplingPts, sampleI)
    {
        operator[](sampleI) = samplingPts[sampleI];
    }
389 390
    curveDist_ = samplingCurveDist;

391 392 393 394 395 396 397 398 399 400 401 402
    cells_ = samplingCells;
    faces_ = samplingFaces;
    segments_ = samplingSegments;
}


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

Foam::sampledSet::sampledSet
(
    const word& name,
    const polyMesh& mesh,
403
    const meshSearch& searchEngine,
404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419
    const word& axis
)
:
    coordSet(name, axis),
    mesh_(mesh),
    searchEngine_(searchEngine),
    segments_(0),
    cells_(0),
    faces_(0)
{}


Foam::sampledSet::sampledSet
(
    const word& name,
    const polyMesh& mesh,
420
    const meshSearch& searchEngine,
421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440
    const dictionary& dict
)
:
    coordSet(name, dict.lookup("axis")),
    mesh_(mesh),
    searchEngine_(searchEngine),
    segments_(0),
    cells_(0),
    faces_(0)
{}


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

Foam::sampledSet::~sampledSet()
{}


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

441
Foam::autoPtr<Foam::sampledSet> Foam::sampledSet::New
442 443 444
(
    const word& name,
    const polyMesh& mesh,
445
    const meshSearch& searchEngine,
446 447 448
    const dictionary& dict
)
{
449
    const word sampleType(dict.lookup("type"));
450 451 452 453

    wordConstructorTable::iterator cstrIter =
        wordConstructorTablePtr_->find(sampleType);

454
    if (!cstrIter.found())
455
    {
456 457
        FatalErrorInFunction
            << "Unknown sample type "
458
            << sampleType << nl << nl
459
            << "Valid sample types : " << endl
460
            << wordConstructorTablePtr_->sortedToc()
461 462 463 464 465 466 467 468 469 470
            << exit(FatalError);
    }

    return autoPtr<sampledSet>
    (
        cstrIter()
        (
            name,
            mesh,
            searchEngine,
471
            dict.optionalSubDict(sampleType + "Coeffs")
472 473 474 475 476 477 478 479 480
        )
    );
}


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

481
    os  << endl << "\t(celli)\t(facei)" << endl;
482 483 484 485 486 487 488 489 490 491 492 493 494

    forAll(*this, sampleI)
    {
        os  << '\t' << cells_[sampleI]
            << '\t' << faces_[sampleI]
            << endl;
    }

    return os;
}


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