polyLineSet.C 10.9 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
#include "polyLineSet.H"
27
28
29
30
31
32
33
34
35
36
#include "meshSearch.H"
#include "DynamicList.H"
#include "polyMesh.H"

#include "addToRunTimeSelectionTable.H"

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

namespace Foam
{
37
38
    defineTypeNameAndDebug(polyLineSet, 0);
    addToRunTimeSelectionTable(sampledSet, polyLineSet, word);
39
40

    const scalar polyLineSet::tol = 1e-6;
41
42
43
44
45
}


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

46
bool Foam::polyLineSet::trackToBoundary
47
(
48
    passiveParticleCloud& particleCloud,
49
    passiveParticle& singleParticle,
50
51
52
53
    label& sampleI,
    DynamicList<point>& samplingPts,
    DynamicList<label>& samplingCells,
    DynamicList<label>& samplingFaces,
54
    DynamicList<scalar>& samplingCurveDist
55
56
) const
{
57
    particle::TrackingData<passiveParticleCloud> trackData(particleCloud);
58

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

62
    while (true)
63
64
65
66
67
68
69
70
71
72
    {
        // Local geometry info
        const vector offset = sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
        const scalar smallDist = mag(tol*offset);

        point oldPos = trackPt;
        label facei = -1;
        do
        {
            singleParticle.stepFraction() = 0;
73
            singleParticle.track(sampleCoords_[sampleI+1], trackData);
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
        }
        while
        (
            !singleParticle.onBoundary()
         && (mag(trackPt - oldPos) < smallDist)
        );

        if (singleParticle.onBoundary())
        {
            //Info<< "trackToBoundary : reached boundary"
            //    << "  trackPt:" << trackPt << endl;
            if
            (
                mag(trackPt - sampleCoords_[sampleI+1])
              < smallDist
            )
            {
                // Reached samplePt on boundary
                //Info<< "trackToBoundary : boundary. also sampling."
                //    << "  trackPt:" << trackPt << " sampleI+1:" << sampleI+1
                //    << endl;
                samplingPts.append(trackPt);
                samplingCells.append(singleParticle.cell());
                samplingFaces.append(facei);

                // trackPt is at sampleI+1
100
                samplingCurveDist.append(1.0*(sampleI+1));
101
102
103
104
105
106
107
108
109
110
111
112
113
            }
            return true;
        }

        // Reached samplePt in cell
        samplingPts.append(trackPt);
        samplingCells.append(singleParticle.cell());
        samplingFaces.append(-1);

        // Convert trackPt to fraction inbetween sampleI and sampleI+1
        scalar dist =
            mag(trackPt - sampleCoords_[sampleI])
          / mag(sampleCoords_[sampleI+1] - sampleCoords_[sampleI]);
114
        samplingCurveDist.append(sampleI + dist);
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129

        // go to next samplePt
        sampleI++;

        if (sampleI == sampleCoords_.size() - 1)
        {
            // no more samples.
            //Info<< "trackToBoundary : Reached end : sampleI now:" << sampleI
            //    << endl;
            return false;
        }
    }
}


130
void Foam::polyLineSet::calcSamples
131
132
133
134
135
(
    DynamicList<point>& samplingPts,
    DynamicList<label>& samplingCells,
    DynamicList<label>& samplingFaces,
    DynamicList<label>& samplingSegments,
136
    DynamicList<scalar>& samplingCurveDist
137
138
139
140
141
) const
{
    // Check sampling points
    if (sampleCoords_.size() < 2)
    {
142
        FatalErrorInFunction
143
144
145
146
            << "Incorrect sample specification. Too few points:"
            << sampleCoords_ << exit(FatalError);
    }
    point oldPoint = sampleCoords_[0];
147
    for (label sampleI = 1; sampleI < sampleCoords_.size(); sampleI++)
148
149
150
    {
        if (mag(sampleCoords_[sampleI] - oldPoint) < SMALL)
        {
151
            FatalErrorInFunction
152
153
154
155
156
157
158
159
160
161
                << "Incorrect sample specification."
                << " Point " << sampleCoords_[sampleI-1]
                << " at position " << sampleI-1
                << " and point " << sampleCoords_[sampleI]
                << " at position " << sampleI
                << " are too close" << exit(FatalError);
        }
        oldPoint = sampleCoords_[sampleI];
    }

162
163
164
    // Force calculation of cloud addressing on all processors
    const bool oldMoving = const_cast<polyMesh&>(mesh()).moving(false);
    passiveParticleCloud particleCloud(mesh());
165

166
167
168
169
170
171
172
173
174
    // current segment number
    label segmentI = 0;

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

    label sampleI = 0;

    point lastSample(GREAT, GREAT, GREAT);
175
    while (true)
176
177
178
179
180
181
182
183
184
185
186
187
188
189
    {
        // Get boundary intersection
        point trackPt;
        label trackCellI = -1;
        label trackFaceI = -1;

        do
        {
            const vector offset =
                sampleCoords_[sampleI+1] - sampleCoords_[sampleI];
            const scalar smallDist = mag(tol*offset);


            // Get all boundary intersections
190
191
192
193
194
            List<pointIndexHit> bHits = searchEngine().intersections
            (
                sampleCoords_[sampleI],
                sampleCoords_[sampleI+1]
            );
195
196
197
198

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

199
            if (bHits.size())
200
201
202
203
204
205
206
207
208
209
210
211
212
            {
                bPoint = bHits[0].hitPoint();
                bFaceI = bHits[0].index();
            }

            // Get tracking point

            bool isSample =
                getTrackingPoint
                (
                    sampleCoords_[sampleI],
                    bPoint,
                    bFaceI,
213
                    smallDist,
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233

                    trackPt,
                    trackCellI,
                    trackFaceI
                );

            if (isSample && (mag(lastSample - trackPt) > smallDist))
            {
                //Info<< "calcSamples : getTrackingPoint returned valid sample "
                //    << "  trackPt:" << trackPt
                //    << "  trackFaceI:" << trackFaceI
                //    << "  trackCellI:" << trackCellI
                //    << "  sampleI:" << sampleI
                //    << "  dist:" << dist
                //    << endl;

                samplingPts.append(trackPt);
                samplingCells.append(trackCellI);
                samplingFaces.append(trackFaceI);

234
                // Convert sampling position to unique curve parameter. Get
235
236
237
238
                // fraction of distance between sampleI and sampleI+1.
                scalar dist =
                    mag(trackPt - sampleCoords_[sampleI])
                  / mag(sampleCoords_[sampleI+1] - sampleCoords_[sampleI]);
239
                samplingCurveDist.append(sampleI + dist);
240
241
242
243
244
245
246
247
248

                lastSample = trackPt;
            }

            if (trackCellI == -1)
            {
                // No intersection found. Go to next point
                sampleI++;
            }
249
        } while ((trackCellI == -1) && (sampleI < sampleCoords_.size() - 1));
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265

        if (sampleI == sampleCoords_.size() - 1)
        {
            //Info<< "calcSamples : Reached end of samples: "
            //    << "  sampleI now:" << sampleI
            //    << endl;
            break;
        }

        //
        // Segment sampleI .. sampleI+1 intersected by domain
        //

        // Initialize tracking starting from sampleI
        passiveParticle singleParticle
        (
266
            mesh(),
267
268
269
270
271
272
            trackPt,
            trackCellI
        );

        bool bReached = trackToBoundary
        (
273
            particleCloud,
274
            singleParticle,
275
276
277
278
            sampleI,
            samplingPts,
            samplingCells,
            samplingFaces,
279
            samplingCurveDist
280
281
282
        );

        // fill sampleSegments
283
        for (label i = samplingPts.size() - 1; i >= startSegmentI; --i)
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
        {
            samplingSegments.append(segmentI);
        }

        if (!bReached)
        {
            //Info<< "calcSamples : Reached end of samples: "
            //    << "  sampleI now:" << sampleI
            //    << endl;
            break;
        }
        lastSample = singleParticle.position();


        // Find next boundary.
        sampleI++;

        if (sampleI == sampleCoords_.size() - 1)
        {
            //Info<< "calcSamples : Reached end of samples: "
            //    << "  sampleI now:" << sampleI
            //    << endl;
            break;
        }

        segmentI++;

        startSegmentI = samplingPts.size();
    }
313
314

    const_cast<polyMesh&>(mesh()).moving(oldMoving);
315
316
317
}


318
void Foam::polyLineSet::genSamples()
319
320
321
322
323
324
{
    // Storage for sample points
    DynamicList<point> samplingPts;
    DynamicList<label> samplingCells;
    DynamicList<label> samplingFaces;
    DynamicList<label> samplingSegments;
325
    DynamicList<scalar> samplingCurveDist;
326
327
328
329
330
331
332

    calcSamples
    (
        samplingPts,
        samplingCells,
        samplingFaces,
        samplingSegments,
333
        samplingCurveDist
334
335
336
337
338
339
    );

    samplingPts.shrink();
    samplingCells.shrink();
    samplingFaces.shrink();
    samplingSegments.shrink();
340
    samplingCurveDist.shrink();
341
342
343
344
345
346
347

    setSamples
    (
        samplingPts,
        samplingCells,
        samplingFaces,
        samplingSegments,
348
        samplingCurveDist
349
350
351
352
353
354
    );
}


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

355
Foam::polyLineSet::polyLineSet
356
357
358
(
    const word& name,
    const polyMesh& mesh,
359
    const meshSearch& searchEngine,
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
    const word& axis,
    const List<point>& sampleCoords
)
:
    sampledSet(name, mesh, searchEngine, axis),
    sampleCoords_(sampleCoords)
{
    genSamples();

    if (debug)
    {
        write(Info);
    }
}


376
Foam::polyLineSet::polyLineSet
377
378
379
(
    const word& name,
    const polyMesh& mesh,
380
    const meshSearch& searchEngine,
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
    const dictionary& dict
)
:
    sampledSet(name, mesh, searchEngine, dict),
    sampleCoords_(dict.lookup("points"))
{
    genSamples();

    if (debug)
    {
        write(Info);
    }
}


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

398
Foam::polyLineSet::~polyLineSet()
399
400
401
402
{}


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