triSurfaceSearch.C 11.3 KB
Newer Older
1
2
3
4
/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
laurence's avatar
laurence committed
5
    \\  /    A nd           | Copyright (C) 2011-2013 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

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

#include "triSurfaceSearch.H"
#include "triSurface.H"
28
29
#include "PatchTools.H"
#include "volumeType.H"
30

31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //

bool Foam::triSurfaceSearch::checkUniqueHit
(
    const pointIndexHit& currHit,
    const DynamicList<pointIndexHit, 1, 1>& hits,
    const vector& lineVec
) const
{
    // Classify the hit
    label nearType = -1;
    label nearLabel = -1;

    const labelledTri& f = surface()[currHit.index()];

    f.nearestPointClassify
    (
        currHit.hitPoint(),
        surface().points(),
        nearType,
        nearLabel
    );

    if (nearType == 1)
    {
        // near point
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79

        const label nearPointI = f[nearLabel];

        const labelList& pointFaces =
            surface().pointFaces()[surface().meshPointMap()[nearPointI]];

        forAll(pointFaces, pI)
        {
            const label pointFaceI = pointFaces[pI];

            if (pointFaceI != currHit.index())
            {
                forAll(hits, hI)
                {
                    const pointIndexHit& hit = hits[hI];

                    if (hit.index() == pointFaceI)
                    {
                        return false;
                    }
                }
            }
        }
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
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
129
130
    }
    else if (nearType == 2)
    {
        // near edge
        // check if the other face of the edge is already hit

        const labelList& fEdges = surface().faceEdges()[currHit.index()];

        const label edgeI = fEdges[nearLabel];

        const labelList& edgeFaces = surface().edgeFaces()[edgeI];

        forAll(edgeFaces, fI)
        {
            const label edgeFaceI = edgeFaces[fI];

            if (edgeFaceI != currHit.index())
            {
                forAll(hits, hI)
                {
                    const pointIndexHit& hit = hits[hI];

                    if (hit.index() == edgeFaceI)
                    {
                        // Check normals
                        const vector currHitNormal =
                            surface().faceNormals()[currHit.index()];

                        const vector existingHitNormal =
                            surface().faceNormals()[edgeFaceI];

                        const label signCurrHit =
                            pos(currHitNormal & lineVec);

                        const label signExistingHit =
                            pos(existingHitNormal & lineVec);

                        if (signCurrHit == signExistingHit)
                        {
                            return false;
                        }
                    }
                }
            }
        }
    }

    return true;
}


131
// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
132

133
134
135
136
137
138
139
Foam::triSurfaceSearch::triSurfaceSearch(const triSurface& surface)
:
    surface_(surface),
    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
    maxTreeDepth_(10),
    treePtr_(NULL)
{}
140
141


142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
Foam::triSurfaceSearch::triSurfaceSearch
(
    const triSurface& surface,
    const dictionary& dict
)
:
    surface_(surface),
    tolerance_(indexedOctree<treeDataTriSurface>::perturbTol()),
    maxTreeDepth_(10),
    treePtr_(NULL)
{
    // Have optional non-standard search tolerance for gappy surfaces.
    if (dict.readIfPresent("tolerance", tolerance_) && tolerance_ > 0)
    {
        Info<< "    using intersection tolerance " << tolerance_ << endl;
    }
158

159
160
161
162
163
164
165
166
167
168
169
170
171
172
    // Have optional non-standard tree-depth to limit storage.
    if (dict.readIfPresent("maxTreeDepth", maxTreeDepth_) && maxTreeDepth_ > 0)
    {
        Info<< "    using maximum tree depth " << maxTreeDepth_ << endl;
    }
}


Foam::triSurfaceSearch::triSurfaceSearch
(
    const triSurface& surface,
    const scalar tolerance,
    const label maxTreeDepth
)
173
174
:
    surface_(surface),
175
176
    tolerance_(tolerance),
    maxTreeDepth_(maxTreeDepth),
177
    treePtr_(NULL)
178
179
180
181
182
183
{}


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

Foam::triSurfaceSearch::~triSurfaceSearch()
184
{
185
186
187
188
189
190
191
    clearOut();
}


void Foam::triSurfaceSearch::clearOut()
{
    treePtr_.clear();
192
193
194
195
196
}


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

197
198
199
200
201
202
const Foam::indexedOctree<Foam::treeDataTriSurface>&
Foam::triSurfaceSearch::tree() const
{
    if (treePtr_.empty())
    {
        // Calculate bb without constructing local point numbering.
203
        treeBoundBox bb(vector::zero, vector::zero);
204

205
        if (surface().size())
206
        {
207
208
            label nPoints;
            PatchTools::calcBounds(surface(), bb, nPoints);
209

210
211
212
213
214
215
216
217
218
            if (nPoints != surface().points().size())
            {
                WarningIn("triSurfaceSearch::tree() const")
                    << "Surface does not have compact point numbering."
                    << " Of " << surface().points().size()
                    << " only " << nPoints
                    << " are used. This might give problems in some routines."
                    << endl;
            }
219

220
221
222
223
224
225
226
227
228
            // Random number generator. Bit dodgy since not exactly random ;-)
            Random rndGen(65431);

            // Slightly extended bb. Slightly off-centred just so on symmetric
            // geometry there are less face/edge aligned items.
            bb = bb.extend(rndGen, 1e-4);
            bb.min() -= point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
            bb.max() += point(ROOTVSMALL, ROOTVSMALL, ROOTVSMALL);
        }
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251

        scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
        indexedOctree<treeDataTriSurface>::perturbTol() = tolerance_;

        treePtr_.reset
        (
            new indexedOctree<treeDataTriSurface>
            (
                treeDataTriSurface(true, surface_, tolerance_),
                bb,
                maxTreeDepth_,  // maxLevel
                10,             // leafsize
                3.0             // duplicity
            )
        );

        indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
    }

    return treePtr_();
}


252
// Determine inside/outside for samples
253
254
255
256
Foam::boolList Foam::triSurfaceSearch::calcInside
(
    const pointField& samples
) const
257
258
259
260
261
262
263
{
    boolList inside(samples.size());

    forAll(samples, sampleI)
    {
        const point& sample = samples[sampleI];

mattijs's avatar
mattijs committed
264
        if (!tree().bb().contains(sample))
265
266
267
        {
            inside[sampleI] = false;
        }
268
        else if (tree().getVolumeType(sample) == volumeType::INSIDE)
269
270
271
272
273
274
275
276
277
278
279
280
        {
            inside[sampleI] = true;
        }
        else
        {
            inside[sampleI] = false;
        }
    }
    return inside;
}


281
void Foam::triSurfaceSearch::findNearest
282
283
(
    const pointField& samples,
284
285
    const scalarField& nearestDistSqr,
    List<pointIndexHit>& info
286
287
) const
{
288
289
    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
290

291
    const indexedOctree<treeDataTriSurface>& octree = tree();
292

293
    info.setSize(samples.size());
294

295
    forAll(samples, i)
296
    {
297
298
299
300
301
302
        static_cast<pointIndexHit&>(info[i]) = octree.findNearest
        (
            samples[i],
            nearestDistSqr[i],
            treeDataTriSurface::findNearestOp(octree)
        );
303
304
    }

305
    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
306
307
308
}


309
Foam::pointIndexHit Foam::triSurfaceSearch::nearest
310
(
311
    const point& pt,
312
    const vector& span
313
314
)
const
315
{
mattijs's avatar
mattijs committed
316
    const scalar nearestDistSqr = 0.25*magSqr(span);
317

318
319
    return tree().findNearest(pt, nearestDistSqr);
}
320
321


322
323
324
325
326
327
328
329
void Foam::triSurfaceSearch::findLine
(
    const pointField& start,
    const pointField& end,
    List<pointIndexHit>& info
) const
{
    const indexedOctree<treeDataTriSurface>& octree = tree();
330

331
332
333
334
335
336
337
338
339
340
341
342
    info.setSize(start.size());

    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();

    forAll(start, i)
    {
        static_cast<pointIndexHit&>(info[i]) = octree.findLine
        (
            start[i],
            end[i]
        );
343
344
    }

345
    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
346
347
348
}


349
void Foam::triSurfaceSearch::findLineAny
350
(
351
352
353
354
    const pointField& start,
    const pointField& end,
    List<pointIndexHit>& info
) const
355
{
356
    const indexedOctree<treeDataTriSurface>& octree = tree();
357

358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
    info.setSize(start.size());

    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();

    forAll(start, i)
    {
        static_cast<pointIndexHit&>(info[i]) = octree.findLineAny
        (
            start[i],
            end[i]
        );
    }

    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
373
374
375
}


376
377
void Foam::triSurfaceSearch::findLineAll
(
378
379
380
381
    const pointField& start,
    const pointField& end,
    List<List<pointIndexHit> >& info
) const
382
{
383
    const indexedOctree<treeDataTriSurface>& octree = tree();
384

385
    info.setSize(start.size());
386

387
388
    scalar oldTol = indexedOctree<treeDataTriSurface>::perturbTol();
    indexedOctree<treeDataTriSurface>::perturbTol() = tolerance();
389

390
391
    // Work array
    DynamicList<pointIndexHit, 1, 1> hits;
392

393
    DynamicList<label> shapeMask;
394

395
    treeDataTriSurface::findAllIntersectOp allIntersectOp(octree, shapeMask);
396

397
398
399
400
    forAll(start, pointI)
    {
        hits.clear();
        shapeMask.clear();
401

402
403
        while (true)
        {
404
            // See if any intersection between pt and end
405
406
407
408
409
410
            pointIndexHit inter = octree.findLine
            (
                start[pointI],
                end[pointI],
                allIntersectOp
            );
411

412
            if (inter.hit())
413
            {
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
                vector lineVec = end[pointI] - start[pointI];
                lineVec /= mag(lineVec) + VSMALL;

                if
                (
                    checkUniqueHit
                    (
                        inter,
                        hits,
                        lineVec
                    )
                )
                {
                    hits.append(inter);
                }
429

430
                shapeMask.append(inter.index());
431
432
433
            }
            else
            {
434
                break;
435
436
            }
        }
437
438

        info[pointI].transfer(hits);
439
    }
440
441

    indexedOctree<treeDataTriSurface>::perturbTol() = oldTol;
442
443
444
}


445
// ************************************************************************* //