uniformSet.C 12.4 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
5
    \\  /    A nd           | Copyright (C) 2011-2016 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
31
32
33
34
35
36
37
38

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

#include "uniformSet.H"
#include "meshSearch.H"
#include "DynamicList.H"
#include "polyMesh.H"

#include "addToRunTimeSelectionTable.H"

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

namespace Foam
{
    defineTypeNameAndDebug(uniformSet, 0);
    addToRunTimeSelectionTable(sampledSet, uniformSet, word);
39
40

    const scalar uniformSet::tol = 1e-3;
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
}


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

bool Foam::uniformSet::nextSample
(
    const point& currentPt,
    const vector& offset,
    const scalar smallDist,
    point& samplePt,
    label& sampleI
) const
{
    bool pointFound = false;

    const vector normOffset = offset/mag(offset);

    samplePt += offset;
    sampleI++;

62
    for (; sampleI < nPoints_; sampleI++)
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
    {
        scalar s = (samplePt - currentPt) & normOffset;

        if (s > -smallDist)
        {
            // samplePt is close to or beyond currentPt -> use it
            pointFound = true;

            break;
        }
        samplePt += offset;
    }

    return pointFound;
}


bool Foam::uniformSet::trackToBoundary
(
82
    passiveParticleCloud& particleCloud,
83
    passiveParticle& singleParticle,
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
    point& samplePt,
    label& sampleI,
    DynamicList<point>& samplingPts,
    DynamicList<label>& samplingCells,
    DynamicList<label>& samplingFaces,
    DynamicList<scalar>& samplingCurveDist
) const
{
    // distance vector between sampling points
    const vector offset = (end_ - start_)/(nPoints_ - 1);
    const vector smallVec = tol*offset;
    const scalar smallDist = mag(smallVec);

    // Alias
    const point& trackPt = singleParticle.position();

100
    particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
101

102
    while(true)
103
104
105
106
107
108
109
    {
        // Find next samplePt on/after trackPt. Update samplePt, sampleI
        if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
        {
            // no more samples.
            if (debug)
            {
110
                Pout<< "trackToBoundary : Reached end : samplePt now:"
111
112
113
114
115
116
117
118
119
120
                    << samplePt << "  sampleI now:" << sampleI << endl;
            }
            return false;
        }

        if (mag(samplePt - trackPt) < smallDist)
        {
            // trackPt corresponds with samplePt. Store and use next samplePt
            if (debug)
            {
121
                Pout<< "trackToBoundary : samplePt corresponds to trackPt : "
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
                    << "  trackPt:" << trackPt << "  samplePt:" << samplePt
                    << endl;
            }

            samplingPts.append(trackPt);
            samplingCells.append(singleParticle.cell());
            samplingFaces.append(-1);
            samplingCurveDist.append(mag(trackPt - start_));

            // go to next samplePt
            if (!nextSample(trackPt, offset, smallDist, samplePt, sampleI))
            {
                // no more samples.
                if (debug)
                {
137
                    Pout<< "trackToBoundary : Reached end : "
138
139
140
141
142
143
144
145
146
147
148
149
                        << "  samplePt now:" << samplePt
                        << "  sampleI now:" << sampleI
                        << endl;
                }

                return false;
            }
        }


        if (debug)
        {
150
            Pout<< "Searching along trajectory from "
151
152
153
154
155
156
157
158
159
160
                << "  trackPt:" << trackPt
                << "  trackCellI:" << singleParticle.cell()
                << "  to:" << samplePt << endl;
        }

        point oldPos = trackPt;
        label facei = -1;
        do
        {
            singleParticle.stepFraction() = 0;
161
            singleParticle.track(samplePt, trackData);
162
163
164

            if (debug)
            {
165
                Pout<< "Result of tracking "
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
                    << "  trackPt:" << trackPt
                    << "  trackCellI:" << singleParticle.cell()
                    << "  trackFaceI:" << singleParticle.face()
                    << "  onBoundary:" << singleParticle.onBoundary()
                    << "  samplePt:" << samplePt
                    << "  smallDist:" << smallDist
                    << endl;
            }
        }
        while
        (
            !singleParticle.onBoundary()
         && (mag(trackPt - oldPos) < smallDist)
        );

        if (singleParticle.onBoundary())
        {
183
            //Pout<< "trackToBoundary : reached boundary" << endl;
184
185
            if (mag(trackPt - samplePt) < smallDist)
            {
186
                //Pout<< "trackToBoundary : boundary is also sampling point"
187
188
189
190
191
192
193
194
195
196
197
                //    << endl;
                // Reached samplePt on boundary
                samplingPts.append(trackPt);
                samplingCells.append(singleParticle.cell());
                samplingFaces.append(facei);
                samplingCurveDist.append(mag(trackPt - start_));
            }

            return true;
        }

198
        //Pout<< "trackToBoundary : reached internal sampling point" << endl;
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
        // Reached samplePt in cell or on internal face
        samplingPts.append(trackPt);
        samplingCells.append(singleParticle.cell());
        samplingFaces.append(-1);
        samplingCurveDist.append(mag(trackPt - start_));

        // go to next samplePt
    }
}


void Foam::uniformSet::calcSamples
(
    DynamicList<point>& samplingPts,
    DynamicList<label>& samplingCells,
    DynamicList<label>& samplingFaces,
    DynamicList<label>& samplingSegments,
    DynamicList<scalar>& samplingCurveDist
) const
{
    // distance vector between sampling points
    if ((nPoints_ < 2) || (mag(end_ - start_) < SMALL))
    {
222
        FatalErrorInFunction
223
224
225
226
227
228
229
230
231
232
233
234
235
            << "Incorrect sample specification. Either too few points or"
            << " start equals end point." << endl
            << "nPoints:" << nPoints_
            << "  start:" << start_
            << "  end:" << end_
            << exit(FatalError);
    }

    const vector offset = (end_ - start_)/(nPoints_ - 1);
    const vector normOffset = offset/mag(offset);
    const vector smallVec = tol*offset;
    const scalar smallDist = mag(smallVec);

236
    // Force calculation of cloud addressing on all processors
237
    const bool oldMoving = const_cast<polyMesh&>(mesh()).moving(false);
238
    passiveParticleCloud particleCloud(mesh());
239

240
    // Get all boundary intersections
241
242
243
244
245
    List<pointIndexHit> bHits = searchEngine().intersections
    (
        start_ - smallVec,
        end_ + smallVec
    );
246
247
248
249

    point bPoint(GREAT, GREAT, GREAT);
    label bFaceI = -1;

250
    if (bHits.size())
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
    {
        bPoint = bHits[0].hitPoint();
        bFaceI = bHits[0].index();
    }

    // Get first tracking point. Use bPoint, bFaceI if provided.

    point trackPt;
    label trackCellI = -1;
    label trackFaceI = -1;

    bool isSample =
        getTrackingPoint
        (
            start_,
            bPoint,
            bFaceI,
268
            smallDist,
269
270
271
272
273
274
275
276
277
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
304
305
306
307

            trackPt,
            trackCellI,
            trackFaceI
        );

    if (trackCellI == -1)
    {
        // Line start_ - end_ does not intersect domain at all.
        // (or is along edge)
        // Set points and cell/face labels to empty lists

        return;
    }

    if (isSample)
    {
        samplingPts.append(start_);
        samplingCells.append(trackCellI);
        samplingFaces.append(trackFaceI);
        samplingCurveDist.append(0.0);
    }

    //
    // Track until hit end of all boundary intersections
    //

    // current segment number
    label segmentI = 0;

    // starting index of current segment in samplePts
    label startSegmentI = 0;

    label sampleI = 0;
    point samplePt = start_;

    // index in bHits; current boundary intersection
    label bHitI = 1;

308
    while(true)
309
310
    {
        // Initialize tracking starting from trackPt
311
        passiveParticle singleParticle(mesh(), trackPt, trackCellI);
312
313
314

        bool reachedBoundary = trackToBoundary
        (
315
            particleCloud,
316
            singleParticle,
317
318
319
320
321
322
323
324
325
            samplePt,
            sampleI,
            samplingPts,
            samplingCells,
            samplingFaces,
            samplingCurveDist
        );

        // fill sampleSegments
326
        for (label i = samplingPts.size() - 1; i >= startSegmentI; --i)
327
328
329
330
331
332
333
334
335
        {
            samplingSegments.append(segmentI);
        }


        if (!reachedBoundary)
        {
            if (debug)
            {
336
                Pout<< "calcSamples : Reached end of samples: "
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
                    << "  samplePt now:" << samplePt
                    << "  sampleI now:" << sampleI
                    << endl;
            }
            break;
        }


        bool foundValidB = false;

        while (bHitI < bHits.size())
        {
            scalar dist =
                (bHits[bHitI].hitPoint() - singleParticle.position())
              & normOffset;

            if (debug)
            {
355
                Pout<< "Finding next boundary : "
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
                    << "bPoint:" << bHits[bHitI].hitPoint()
                    << "  tracking:" << singleParticle.position()
                    << "  dist:" << dist
                    << endl;
            }

            if (dist > smallDist)
            {
                // hitpoint is past tracking position
                foundValidB = true;
                break;
            }
            else
            {
                bHitI++;
            }
        }

        if (!foundValidB)
        {
            // No valid boundary intersection found beyond tracking position
            break;
        }

        // Update starting point for tracking
        trackFaceI = bFaceI;
        trackPt = pushIn(bPoint, trackFaceI);
        trackCellI = getBoundaryCell(trackFaceI);

        segmentI++;

        startSegmentI = samplingPts.size();
    }
389
390

    const_cast<polyMesh&>(mesh()).moving(oldMoving);
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
}


void Foam::uniformSet::genSamples()
{
    // Storage for sample points
    DynamicList<point> samplingPts;
    DynamicList<label> samplingCells;
    DynamicList<label> samplingFaces;
    DynamicList<label> samplingSegments;
    DynamicList<scalar> samplingCurveDist;

    calcSamples
    (
        samplingPts,
        samplingCells,
        samplingFaces,
        samplingSegments,
        samplingCurveDist
    );

    samplingPts.shrink();
    samplingCells.shrink();
    samplingFaces.shrink();
    samplingSegments.shrink();
    samplingCurveDist.shrink();

    setSamples
    (
        samplingPts,
        samplingCells,
        samplingFaces,
        samplingSegments,
        samplingCurveDist
    );
}

428

429
430
431
432
433
434
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

Foam::uniformSet::uniformSet
(
    const word& name,
    const polyMesh& mesh,
435
    const meshSearch& searchEngine,
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
    const word& axis,
    const point& start,
    const point& end,
    const label nPoints
)
:
    sampledSet(name, mesh, searchEngine, axis),
    start_(start),
    end_(end),
    nPoints_(nPoints)
{
    genSamples();

    if (debug)
    {
451
        write(Pout);
452
453
454
455
456
457
458
459
    }
}


Foam::uniformSet::uniformSet
(
    const word& name,
    const polyMesh& mesh,
460
    const meshSearch& searchEngine,
461
462
463
464
465
466
467
468
469
470
471
472
    const dictionary& dict
)
:
    sampledSet(name, mesh, searchEngine, dict),
    start_(dict.lookup("start")),
    end_(dict.lookup("end")),
    nPoints_(readLabel(dict.lookup("nPoints")))
{
    genSamples();

    if (debug)
    {
473
        write(Pout);
474
475
476
477
478
479
480
481
482
483
484
    }
}


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

Foam::uniformSet::~uniformSet()
{}


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